# The Three-Dimensional Mass Distribution in NGC 17001footnote 11footnote 1Observations reported in this paper were obtained at the Multiple Mirror Telescope Observatory, a joint facility of the University of Arizona and the Smithsonian Institution.

## 1 Introduction

Evidence is accumulating for a dichotomy between two classes of elliptical galaxies: high-luminosity, slowly rotating systems with shallow central brightness cusps and a tendency for boxy isophotes, and lower luminosity, rotationally supported systems with steeper central cusps and a tendency for diskiness (Lauer et al. 1995, Kormendy & Bender 1996, Faber et al. 1997, and references therein). Photometric shape differences between high- and low-luminosity ellipticals (Tremblay & Merritt 1996) are roughly consistent with the general expectation that rapidly rotating systems should have oblate axisymmetry while slow rotators are more likely triaxial. There are exciting suggestions that central cusps may drive evolution towards axisymmetry through orbital stochasticity (Merritt & Fridman 1996, Merritt & Valluri 1996). But where the causal arrows actually point among these various properties is still far from clear.

There is consequently a growing need for intrinsic shape measurements for a subsample of well studied ellipticals. This paper continues our efforts to obtain such measurements. Previously we have derived approximate shapes for the “standard” elliptical NGC 3379 (Statler 1994c; hereafter S94c) and NGC 4589 (Statler 1994b; S94b) using photometry and kinematic data from the literature. In this paper we take up the case of the high surface brightness elliptical NGC 1700. The availability of high quality multi-position-angle kinematics allows an intrinsic shape determination of unprecedented precision, as well as interesting constraints on the presence of dark matter in the outer parts of the system.

NGC 1700 is an elliptical galaxy of fairly average luminosity ($`M_B=22.3`$ for $`H_0=50\mathrm{km}\mathrm{s}^1`$ $`\mathrm{Mpc}^1`$), but comparatively high central dispersion ($`230\mathrm{km}\mathrm{s}^1`$) and small ($`14\mathrm{}`$) effective radius. This puts it roughly $`2\sigma `$ off the Fundamental Plane in the sense of having uncharacteristically low $`M/L`$ for objects of this mass. NGC 1700 is notable for being the most luminous of the steep-cusped “power law” ellipticals in the sample observed by HST (Lauer et al. 1995; most recent compilation in Faber et al. 1997). The high surface brightness at $`r_e`$ makes it particularly attractive for a study of intrinsic shape and $`M/L`$ gradient, since long-slit spectra with good signal-to-noise can be obtained out to several optical scale lengths.

Statler, Smecker-Hane, & Cecil (1996, hereafter SSC) present kinematic profiles for NGC 1700 out to $`4r_e`$ on four position angles. This is the first data set to meet the requirements for accurate shape fitting spelled out in S94b; accuracy in the mean velocity field (VF) is better than $`10\mathrm{km}\mathrm{s}^1`$ inside $`2r_e`$ and better than 15% of the maximum inside $`3r_e`$. In addition to the usual profiles of $`v`$, $`\sigma `$, $`h_3`$, and $`h_4`$, SSC derive profiles of the velocity field asymmetry parameters $`(V_1,V_2,V_3)`$ (Statler & Fry 1994), and show all three to be consistent with zero over a factor $`10`$ in radius. From the VF’s high degree of symmetry, SSC surmise that the triaxiality of NGC 1700 is probably small; i.e., that the shape is close to oblate. They also argue qualitatively that the shallow slope of the projected RMS velocity profile is consistent with a standard dark halo. Our objective in this paper is to make both of these statements quantitative.

We proceed as follows: in Section 2 we apply the dynamical shape fitting method of S94b to SSC’s kinematic data and photometry from Franx et al. (1989). We find that NGC 1700 is very nearly oblate. This result lets us exploit existing tools for modeling axisymmetric systems to assess the mass-to-light ratio. In § 3 we demonstrate that two-integral models with constant $`M/L`$ are strongly inconsistent with the data, while a model with $`M/L`$ increasing outward by a factor of $`2`$ over the range of the observations is viable. We then construct anisotropic three-integral models, and show that even anisotropy cannot completely obviate the need for an $`M/L`$ gradient. We summarize and discuss the broader implications of these results in § 4.

## 2 Intrinsic Shape

### 2.1 Method

The shape fitting method is described in detail in S94b and S94c. The underlying dynamical models (Statler 1994a, hereafter S94a) rest on the following assumptions: (1) mean rotation in the galaxy arises from internal streaming in a potential with negligible figure rotation; (2) the mean-motion streamlines of short-axis and long-axis tube orbits can be represented by coordinate lines in a confocal coordinate system;<sup>2</sup><sup>2</sup>2Anderson & Statler (1998) demonstrate the validity of this assumption for realistic potentials and show that the triaxiality parameter recovered from the streamlines closely agrees with the triaxiality of the mass distribution. (3) the luminosity density $`\rho _L`$ is stratified on similar ellipsoids, i.e., $`\rho _L(r,\theta ,\varphi )=\overline{\rho }_L(r)\rho _L^{}(\theta ,\varphi )`$; and (4) the velocity field obeys a “similar flow” ansatz outside the tangent point for a given line of sight, $`𝐯(r,\theta ,\varphi )=\overline{v}(r)𝐯^{}(\theta ,\varphi )`$.

Assumptions (3) and (4) are used only for projecting the models. As long as $`\overline{\rho }_L(r)`$ and $`\overline{\rho }_L(r)\overline{v}(r)`$ decrease faster than $`r^2`$, nearly all of the contribution to the projection integrals will come from radii near the tangent point, and the result will be insensitive to the structure at larger $`r`$. We refer to this approximation as quasi-local projection. The requirement on $`\overline{\rho }_L(r)\overline{v}(r)`$ generally limits the validity of the models to regions where the rotation curve is not steeply rising, which for NGC 1700 means $`rr_e`$.

A single model is described by the parameters $`(T,c_L,\mathrm{\Omega },𝐝)`$, where $`T`$ is the triaxiality of the total mass distribution, $`c_L`$ is the short-to-long axis ratio of the luminosity distribution, $`\mathrm{\Omega }(\theta _E,\varphi _E)`$ is the orientation of the galaxy relative to the line of sight, and the vector $`𝐝`$ represents the remaining dynamical parameters which are defined in S94b. The projected model predicts the ellipticity and the radial velocities on each sampled position angle, from which we calculate the likelihood of the observations. Repeating for $`10^7`$ models covering the full parameter space, we determine the multidimensional likelihood $`L(T,c_L,\mathrm{\Omega },𝐝)`$. The Bayesian estimate of the galaxy’s shape is then the probability distribution $`P(T,c_L)`$, obtained by multiplying the likelihood by a model for the parent distribution from which Nature drew the galaxy and integrating over all other parameters:

$$P(T,c_L)=𝑑\mathrm{\Omega }𝑑𝐝F_p^{}(T,c_L,\mathrm{\Omega },𝐝)L(T,c_L,\mathrm{\Omega },𝐝).$$
(1)   
In Bayesian terminology, the integrand is called the posterior density and $`F_p^{}`$ the prior. The marginal posterior density $`P`$ is the integral of the posterior over the nuisance parameters. We use $`P`$ to stand for any marginal density; for instance, the joint distribution in shape and orientation is given by

$$P(T,c_L,\mathrm{\Omega })=𝑑𝐝F_p^{}(T,c_L,\mathrm{\Omega },𝐝)L(T,c_L,\mathrm{\Omega },𝐝).$$
(2)   
We can also factor the parent distribution according to

$$F_p^{}(T,c_L,\mathrm{\Omega },𝐝)=\frac{1}{4\pi }F_p(T,c_L)F_p^d(𝐝).$$
(3)   
The first factor describes a population with random orientations,<sup>3</sup><sup>3</sup>3Actually, the population of ellipticals may have a preferred orientation if, for instance, a significant number of E’s are misclassified S0’s seen far from edge-on. and the dependence on $`𝐝`$ can be factored out without loss of generality because the $`𝐝`$ themselves can implicitly be functions of $`T`$ and $`c_L`$.

In modeling NGC 3379 and NGC 4589 it was necessary to average the kinematics over a single large radial bin on each PA, yielding only a single average shape for each system. For NGC 1700 we can extend the method because we have excellent data over nearly a factor of 4 in radius in the region consistent with the assumptions of the models. SSC’s radial binning gives 5 independent measurements along each PA, with bin centers at $`12\stackrel{}{\mathrm{.}}6`$, $`15\stackrel{}{\mathrm{.}}6`$, $`21\stackrel{}{\mathrm{.}}0`$, $`30\stackrel{}{\mathrm{.}}6`$, and $`46\stackrel{}{\mathrm{.}}8`$. We first model each annulus independently in the quasi-local approximation. Then, assuming that the intrinsic principal axes remain aligned through the galaxy, we can combine the posterior probabilities requiring that the five annuli match their local kinematics and ellipticities in the same orientation. If we further assume that the triaxialities of the mass and light distributions are equal, we can add the constraint that the triaxiality profile must reproduce the observed isophotal twist. This procedure is described more fully in Appendix A.

We argued previously (SSC) that the parallelogram-shaped isophotal distortions seen beyond $`R60\mathrm{}`$ arise from an incompletely phase-mixed, differentially precessing ring or disk, presumably acquired by accretion. If so, then the photometric and kinematic twists that start near $`R40\mathrm{}`$ may be the result of intrinsic twisting of the isodensity surfaces rather than a shape gradient in an intrinsically aligned system. In fact, we strongly suspect that this is the case; nonetheless, in this paper we proceed under the assumption that the isodensity surfaces are aligned and derive constraints on the shape gradient in that context. If the outermost annulus is actually misaligned, our results for the inner four annuli still hold and are virtually unaffected, as we have verified by recomputing the posterior densities excluding the last annulus.

### 2.2 Data

#### 2.2.1 Photometry

Ellipticity and position angle profiles are taken from Franx et al. (1989) and averaged over the approximate ranges of SSC’s radial bins. At smaller radii the intervals over which the profiles are averaged overlap slightly in order to better reflect the ellipticity gradient inside $`25\mathrm{}`$; this has a negligible effect on the adopted profile but (intentionally) increases the error bars where the gradient is steeper. The results are given in columns 2 and 3 of Table 1. Data are weighted according to the inverse square of Franx et al.’s error bars. The tabulated uncertainty is the larger of the mean observational error per data point and the variance over the bin.

The logarithmic slope of the surface brightness profile wavers non-monotonically between about $`1.7`$ and $`2.3`$ for $`10\mathrm{}r55\mathrm{}`$. For simplicity in projection we adopt a power-law radial profile for the luminosity density, $`\overline{\rho }r^k`$, compute models for $`k=(\frac{11}{4},3,\frac{13}{4})`$, and average the results. It turns out that the results are insensitive to the choice of $`k`$.

#### 2.2.2 Kinematics

Mean velocities, corrected for the non-Gaussian shape of the LOSVD, are taken directly from columns 10 and 11 of SSC’s Table 2 (as plotted in their Figure 10), folded about the origin, and averaged. Results are given in columns 4 – 7 of Table 1 for PAs 0 (near minor axis), 225, 270, and 315.

The kinematic sampling is at a different position angle with respect to the photometric axes in each annulus because the major axis PA is not strictly constant in radius. However, the PA twist for $`10\mathrm{}<r<35\mathrm{}`$ is so small that we assume the sampling in the inner four annuli is the same for the purpose of calculating the projected model velocities; this saves a factor of 3 in computing time.

Over the region being modeled, the observed mean rotation amplitude is constant or slowly increasing with radius, though on some PAs the rotation curve does fall at large $`R`$. Again for simplicity in projection, we take a power-law form for the velocity scaling law $`\overline{v}(r)r^l`$ and compute models for $`l=(0,\pm \frac{1}{2})`$. As for the density index $`k`$, the results are not sensitive to the value of $`l`$, and we simply average the derived probability densities.

#### 2.2.3 Forbidden Velocities on PA 225?

SSC reported a striking reversal in the rotation profile along PA 225, at $`r40\mathrm{}`$ to the northeast, and argued that this feature, most prominent in the outermost data point, is real. Two of us (TSS and TS-H) subsequently obtained additional spectra along this PA with the MMT Red Channel on 1995 February 3–4 UT. The instrument was rotated $`180\mathrm{°}`$ relative to the earlier setup, and the galaxy was displaced $`60\mathrm{}`$ to the southwest, putting the region of interest near the center of the slit. Reductions were performed as described by SSC.

The resulting rotation curve shows no velocity reversal. We now believe the earlier result to have been spurious, caused by a distortion of the cross-correlation peak by some unnoticed or inadequately corrected systematic effect (possibly scattered light) in the weakly exposed outer part of the image.

The data were reduced, unfortunately, after the present models had all been computed, but the results for intrinsic shape are not compromised. The radius at which the velocity had appeared to change sign is outside the region we are modeling. The revised rotation profile alters only our outermost annulus, through one data point which is averaged with its counterpart on the opposite side of the galaxy. Were we to include the revision, the value of $`v_{225}`$ in the last row of Table 1 would change from $`53.1\pm 26.9\mathrm{km}\mathrm{s}^1`$ to approximately $`84\pm 21\mathrm{km}\mathrm{s}^1`$. This is a change of just over $`1\sigma `$ in one datum (with the second-largest error bar) out of 20 in the velocity field. We have rerun models over part of the parameter space with this change included and find only small differences; we will point out below where the effect is noticeable.

### 2.3 Results

Models are computed over essentially the same grid of dynamical parameters used previously (S94b, S94c). At each of the 9 pairs of $`(k,l)`$ values there are 8 choices of the “contrast function” (giving the relative streaming amplitude in short-axis and long-axis tubes), 4 choices of the function $`v^{}(t)`$ (giving the angular dependence of the mean velocity away from the symmetry planes), and two treatments of radial mean motions, with streamlines assumed to lie on either spherical or ellipsoidal shells. The roles of the contrast and $`v^{}(t)`$ functions are summarized in Appendix B. We do not calculate models in which the mass and luminosity distributions have different triaxialities since our earlier tests showed that the triaxiality of the latter has little bearing on that of the former. When we include the alignment and isophotal twist constraints we implicitly assume the two triaxialities are equal; this is justifiable a posteriori since nearly axisymmetric potentials turn out strongly preferred. By Monte Carlo sampling over $`(T,c,\mathrm{\Omega })`$,<sup>4</sup><sup>4</sup>4From here on we omit the $`L`$ subscript on $`c`$, but it is still the axis ratio of the luminosity distribution. the likelihood $`L(T,c,\mathrm{\Omega },𝐝)`$ for each annulus is computed for each set of dynamical parameters $`𝐝`$, and the individual likelihoods are combined by the methods described in Section 2.1.

#### 2.3.1 The “Maximal Ignorance” Estimate

The marginal posterior densities $`P_{12345}(T_i,c_i)`$, including the assumption of intrinsic alignment and the isophotal twist constraint (cf. equation 19), are shown in Figure 1 for all models. This should be considered the “maximal ignorance” estimate of intrinsic shape, since as few restrictions as possible are placed on the dynamics. In each panel, oblate spheroids are at the right margin and prolate spheroids at the left. We rather conservatively allow values of $`c`$ down to $`0.35`$, but since no models with $`c>0.8`$ are consistent, at even the $`4\sigma `$ level, with the observed $`ϵ>0.26`$, the upper part of the parameter space is truncated. White contours enclosing $`68\%`$ and $`95\%`$ of the total probability are the Bayesian equivalent of $`1\sigma `$ and $`2\sigma `$ error bars.

The features to note in Figure 1 are: (1) the strong similarity of the distributions for the inner 4 annuli, with prominent peaks at $`T=0`$; (2) the overall trend toward smaller $`c`$ with increasing radius, consistent with the gradual radial increase in $`ϵ`$; and (3) the significantly broader distribution in the outermost annulus, indicating an increase in $`T`$. The diagonal ridges in the distribution and the scalloped contours toward smaller $`c`$ are numerical artifacts arising from undersampling of the sharply peaked densities $`P_i(T_i,c_i,\mathrm{\Omega }_i)`$. Changing the outermost velocity on PA 225 (Sec. 2.2.3) causes a minuscule shift of the distributions for the inner four annuli toward smaller $`T`$ and larger $`c`$, and in the fifth annuls flattens the secondary maximum near $`T=0.2`$, $`c=0.45`$ and truncates the weak appendage extending toward $`T=1`$, $`c=0.7`$.

Figure 1 suggests that NGC 1700 is likely to be close to axisymmetric, at least in the inner $`35\mathrm{}`$ or so. Considering its central cuspiness, it is important to derive a quantitative constraint on the triaxiality $`T`$. By integrating the distributions in Figure 1 over $`c`$, we obtain the marginal densities for $`T`$ shown in the first two panels of Figure 2. Figure 2$`a`$ shows the differential distributions; the corresponding cumulative distributions are plotted in Figure 2$`b`$. The distributions for the inner 4 annuli are so similar that one is naturally led to suspect that $`T`$ may be constant over this range of radii. Indeed, plotting the integrated distributions tends somewhat to mask a strong preference in the joint distribution for constant $`T`$. This is demonstrated in Figure 2$`c`$, in which we plot the differential probabilities for the change in triaxiality between annuli 1 and 2, 2 and 3, and 3 and 4. All are sharply peaked at $`\mathrm{\Delta }T=0`$. The strong similarity of these distributions and the small scatter in mean values are consistent with $`\mathrm{\Delta }T`$ being drawn from a parent distribution that is a delta function at $`\mathrm{\Delta }T=0`$, indicating constant $`T`$ in the inner 4 annuli.

If $`T`$ is assumed to be the same in annuli 1 – 4, then this assumption can be included as “prior information” (essentially by inserting factors of the form $`\delta (T_iT_{i+1})\delta (T_{i+1}T_{i+2})\mathrm{}`$ in equation ). The resulting cumulative distribution for the single value of $`T`$ is shown in figure 2$`d`$. Upper limits on $`T`$ at the $`1\sigma `$ and $`2\sigma `$ confidence levels are drawn as dotted lines. We obtain

$$T<0.16\mathrm{for}2.2h^1\mathrm{kpc}<r<6.7h^1\mathrm{kpc}(95\%\mathrm{confidence}).$$
(4)   
A more intuitive but less precise statement is that the middle-to-long axis ratio $`b/a`$ is greater than $`0.9`$ at $`>95\%`$ confidence.

#### 2.3.2 Statistical Effects of Orientation

Because $`P(T,c)`$ is a marginal density integrated over orientation $`\mathrm{\Omega }`$, its value at a given $`(T,c)`$ is proportional to the amount of solid angle over which a model of that shape will be seen to reproduce the observations (averaged over the other parameters). Consequently, there is a complementarity between constraints on shape and orientation. Shapes able to fit the data in many orientations are statistically favored while those requiring special orientations are not; orientations that are required by many shapes earn high marginal probabilities while those useful to only a few score low. This leads to the somewhat counterintuitive result that constraining the models to be at the same orientation has the effect of narrowing the marginal distribution in orientation while broadening that in intrinsic shape. This effect can be rigorously demonstrated for analogous one-dimensional problems. For NGC 1700’s highly symmetric velocity field, triaxial models requiring lines of sight in or near the $`xz`$ plane to suppress apparent minor axis rotation are made relatively more likely by the alignment constraint than axisymmetric models which can be more freely oriented. The resulting marginal densities $`P(\mathrm{\Omega })`$ thus show a strong preference for such orientations.

One should remember, however, that this estimate is made in isolation, nearly independent of knowledge of the rest of the elliptical galaxy population. That the data may imply a special orientation for a single object is hardly alarming. It is only in conjunction with other objects that the statistics are affected by the requirement that the observed sample must be consistent with a random drawing from an isotropic parent distribution. Since most observed ellipticals are known not to have strongly triaxial-like kinematics, this added constraint will push down the peak in the orientation distribution, eliminating triaxial models that require that special orientation and narrowing the marginal distribution for $`T`$. Methods for including such constraints have been developed (Bak & Statler 1997) and will be discussed fully and applied in a forthcoming paper (Bak & Statler 1998, in preparation). For now it suffices to acknowledge that including such constraints in the modeling of NGC 1700 would tighten the limits on $`T`$, so that the limit in equation (4) can be regarded as conservative.

#### 2.3.3 The Dynamical Prior

The maximal ignorance estimate intentionally places as few restrictions as possible on the orbit populations, because we have yet to understand how Nature chooses to populate orbits. Such knowledge may be forthcoming, however, either from observations or numerical simulations. Prior restrictions on the dynamics will alter the shape estimate.

The infrequent occurrence of minor axis rotation in ellipticals seems to indicate that the “intrinsic misalignment” between the angular momentum vector and the short axis tends to be small (Franx, Illingworth, & de Zeeuw 1991). It may be the case that net streaming motions are present only in the short-axis tube orbits. If all models with long-axis tube streaming are excluded from the fits, the result is the set of posterior densities shown in Figure 3. The galaxy comes out more triaxial, since, with long-axis tube streaming disallowed, triaxial systems are more likely to be seen in orientations where the apparent minor axis rotation is zero. In Figure 4$`a`$ we show the marginal densities for $`T`$ alone. Those for the inner four annuli are not quite the same, but the differences, compared with Figure 2$`a`$, are consistent with the larger amount of numerical noise in the restricted set of models. Adding the prior constraint that $`T`$ is constant over this range of radii gives the cumulative distribution in Figure 4$`b`$, and a 95%-confidence upper limit of $`T<0.22`$. Note that this is a change of only $`0.06`$ from the maximal ignorance limit.

By far the most important dynamical parameter is not the relative loading of long- and short-axis tubes, but the form of the $`v^{}(t)`$ function, which gives the latitude dependence of the mean streaming velocity across the $`xz`$ plane. Restricting the models to $`v^{}(t)`$ types 1 and 4 or to types 2 and 3 gives the marginal $`T`$ probabilities shown in figure 5. (No restriction is placed on long-axis tube streaming.) The former case shows a marked preference for triaxial shapes over axisymmetric ones, the latter just the opposite.

The reason for this difference is the way in which the models fit the position-angle dependence of the observed velocity field. SSC find that velocities on the sampled PAs are consistent with an almost exactly sinusoidal dependence at a given radius. To understand how this affects the models, we note first that the results for types 1 and 4 separately are nearly the same, as are those for types 2 and 3; this similarity indicates that the model velocity fields are dominated by short-axis tubes. In triaxial systems, the sequence of short-axis tubes terminates in the $`xz`$ plane at a latitude $`\mathrm{cos}^1T^{1/2}`$ up from the $`xy`$ plane. The streamlines diverge to either side of this constriction (see fig. 2b of S94a). The mean flow in the vicinity of the short axis will be seriously diluted, either by this divergence alone, or by additional geometrical factors if the line of sight is not close to the $`xz`$ plane — unless either the flow velocity through the constriction is sufficiently large near the terminating latitude, or $`T`$ is sufficiently small that the terminating latitude is itself close to the $`z`$ axis. The type 2 and 3 models, where the flow velocity through the constriction falls off smoothly to zero with increasing latitude, rotate too slowly on the diagonal position angles relative to the major axis in projection, unless they are so close to oblate that the short-axis tubes extend nearly all the way to the $`z`$ axis. Conversely, the type 1 and 4 models, where the velocity is assumed constant right up to the terminating latitude, rotate too fast on the diagonal PAs unless the terminating latitude is rather low, meaning $`T`$ is significantly above zero.

Knowing the actual rate of occurrence of triaxiality in ellipticals has major consequences for the abundance of steep central cusps and/or massive black holes, and is an important test of the idea that central mass concentrations can destroy triaxiality. The above result implies that triaxiality constraints could be considerably strengthened if we could restrict, on physical grounds, the form of $`v^{}(t)`$. It is, therefore, extremely important to obtain the $`v^{}(t)`$ function from simulated ellipticals, e.g., from $`N`$-body merger experiments (Dutta, Statler, & Weil 1998). Furthermore, understanding how this function is determined physically by different mechanisms could make it possible, with a sufficiently large sample, to use present-day kinematics to constrain galaxy formation physics.

## 3 Mass-to-Light Ratio and Anisotropy

Given the upper limit on $`T`$ derived in the previous section, we can safely use axisymmetric models to estimate the mass-to-light ratio in NGC 1700. A wealth of techniques have been developed over the past several years to fit axisymmetric models to observations, of which two suffice for our purposes.

### 3.1 Moment Models

#### 3.1.1 Method

For axisymmetric systems in which the phase space distribution function $`f`$ depends only on energy and the $`z`$ component of angular momentum, the solution of the Jeans equations for the low order moments of $`f`$ can be broken down into one-dimensional integrals along lines parallel to the symmetry axis and derivatives of these integrals (Satoh 1980; Binney, Davies, & Illingworth 1990, hereafter BDI). Efficient methods for constructing such models were developed by BDI and refined by van der Marel, Binney, & Davies (1990; hereafter MBD) and Magorrian & Binney (1994). Our implementation essentially parallels that of BDI and MBD, with minor modifications.

The $`R`$-band surface photometry is first interpolated onto an elliptic grid, spaced uniformly in position angle $`\psi `$ and logarithmically in elliptic radius $`m_p=R(\mathrm{cos}^2\psi +\mathrm{sin}^2\psi /c_p^2)^{1/2}`$; here $`R`$ is projected radius on the sky and we take $`c_p=0.7`$. Since we are assuming axisymmetry, PA twists are ignored. The photometry is extrapolated to small and large radii using separate $`r^{1/4}`$-law fits to the inner and outer parts of the radial profile.

The deprojection for a given inclination $`i`$ starts with a parametric least-squares fit of a luminosity density of the form

$$\rho _L(m)=\frac{\rho _0a^3}{m^2(m^2+a^2)^{1/2}}.$$
(5)   
In equation (5), $`a`$ is a scale length and $`m`$ is the spheroidal radius $`r(\mathrm{cos}^2\theta +\mathrm{sin}^2\theta /c^2)^{1/2}`$, where $`c=[(c_p^2\mathrm{cos}^2i)/\mathrm{sin}i]^{1/2}`$ is the intrinsic axis ratio that projects to the apparent ratio $`c_p`$ at inclination $`i`$. This model projects to a surface brightness

$$I(m_p)=\frac{I_0a}{m_p}\mathrm{tan}^1\frac{a}{m_p},$$
(6)   
which is a better fit to the NGC 1700 data than the Hernquist and Jaffe models used by MBD. The model is evaluated on a spheroidal grid, and the fit is then improved by iterating via Lucy’s Method, giving a final luminosity model that is still axisymmetric but not spheroidal. Three iterations are sufficient to bring the photometric residuals down below $`0.05`$ magnitudes.

The mass density is found from

$$\rho (r,\theta )=\mathrm{{\rm Y}}\rho _L(r,\theta ),$$
(7)   
where we can take the mass-to-light ratio $`\mathrm{{\rm Y}}`$ to be either constant or a function of $`m`$. We obtain the gravitational potential $`\mathrm{\Phi }`$ by a spherical harmonic expansion, calculate the mean-square velocities using BDI’s equations (7) – (8), and project them onto the line of sight. We choose not to separate the ordered rotation from the random tangential velocities since this introduces more free parameters into the models. Instead, we add the observed dispersions and mean velocities in quadrature and compare the resulting RMS with the models.

#### 3.1.2 Results

The dotted lines in figure 6 show the RMS velocity profiles for the best constant-$`M/L`$ model. This model has $`i=90\mathrm{°}`$ (edge-on) and has been fitted by eye to the data in the inner $`10\mathrm{}`$. The model clearly fails to reproduce the high RMS velocities observed beyond $`r_e`$. Note also that for $`R<10\mathrm{}`$, the velocity on the major axis is slightly overestimated and that on the minor axis slightly underestimated. This mismatch only gets worse for flatter models at lower inclinations, since, to support the flattening, the tangential dispersion must increase at the expense of the radial and vertical dispersions (van der Marel 1991). We conclude that 2-integral models with constant $`M/L`$ cannot fit the data.

Two-integral models with variable $`M/L`$ fare better, though it is difficult to constrain the details of the mass profile. To estimate the magnitude of the required $`M/L`$ gradient, we adopt a simple law of the form

$$\mathrm{{\rm Y}}(m)=\mathrm{{\rm Y}}_0\left[1+(m/r_h)^b\right]^{1/b},$$
(8)   
in which $`\mathrm{{\rm Y}}`$ is stratified approximately on density surfaces. For large $`b`$, this describes an $`M/L`$ ratio that is constant at small radii, turns up abruptly at $`r_h`$, and then increases linearly with radius. Since the surface brightness is still falling roughly as $`R^2`$ where the photometry ends, $`\mathrm{{\rm Y}}r`$ corresponds to what one expects from a standard isothermal dark halo. Fixing the value of $`b=8`$ and minimizing $`\chi ^2`$ with respect to $`\mathrm{{\rm Y}}_0`$ and $`r_h`$ results in the model shown as the dashed lines in figure 6; the fit is statistically adequate, with a reduced $`\chi ^2`$ of $`1.16`$ (49 degrees of freedom). This model is inclined at $`90\mathrm{°}`$, since, as before, lower inclinations are unable to match the major to minor axis dispersion ratio. Making the halo spherical rather than spheroidal only exacerbates the problem, since then the galaxy cannot rely on the flattening of the potential to keep the luminosity density flat, and must become even more tangentially anisotropic.

The best fit local and cumulative $`M/L`$ profiles are shown in figure 7. We have converted to blue luminosity using a mean $`BR`$ of $`1.7`$ (FIH). The central value of $`5.5h`$ agrees very well with Bender, Burstein, & Faber’s (1992) result of $`2.8`$ for $`h=0.5`$, and the factor of $`2`$ cumulative increase out to a scale of $`10h^1\mathrm{kpc}`$ is consistent with standard halo models. Thus 2-integral models with dark matter appear able to fit the data.

This result, however, is less than entirely satisfactory for at least three reasons. First, the velocity field fitting of Section 2 implies flatter shapes and lower inclinations than obtained here. Second, the mere existence of a solution to the Jeans equations does not guarantee the existence of a non-negative distribution function $`f`$. And third, since mass can be traded off with anisotropy, 3-integral models are essential to determine whether an $`M/L`$ gradient is actually required by the data.

### 3.2 Quadratic Programming Models

#### 3.2.1 Method

Three-integral models constructed by quadratic programming (QP) score over moment models in that, since one solves explicitly for a positive definite phase space distribution function, there is no risk that an apparently successful model may require an unphysical negative density.

The numerical approach used here is essentially the one described by Dejonghe et al. (1996). In brief, the procedure is as follows: first, for an assumed inclination $`i`$, the surface photometry is deprojected as described above using Lucy’s method. The total mass distribution is assumed to have the same shape as the luminosity density, though not necessarily the same radial profile. This is accomplished by setting, for each term in the spherical harmonic expansion of the density,

$$\rho _{lm}(r)=\rho _{L,lm}(r)\mathrm{{\rm Y}}_0(1+Br^p);$$
(9)   
the case of constant $`M/L`$ then corresponds to $`B=0`$. The potential is computed from $`\rho `$ by harmonic expansion. A Stäckel fit is made to the potential, using the de Zeeuw & Lynden-Bell (1985) method as implemented by Dejonghe & de Zeeuw (1988). This fit is made only to define an approximate third integral $`I_3`$; wherever the potential is needed during the modeling, the original potential is used. The distribution function $`f`$ is built up as a sum of components of Fricke type:

$$f(E,L_z,I_3)=\underset{p,q,n}{}A_{pqn}^\pm F_{pqn}^\pm ,$$
(10)   
where $`n`$ is an integer but $`p`$ and $`q`$ may be real, and

$$F_{pqn}^\pm =\{\begin{array}{cc}E^p(EL_z^2/2)^q(EI_3)^n,\hfill & \pm L_z0,\hfill \\ 0\hfill & \pm L_z<0.\hfill \end{array}$$
(11)   
For components of this type, the velocity moments can be expressed in terms of algebraic functions of the coordinates and the potential. Models are constructed from a library of components with various combinations of $`p`$, $`q`$ and $`n`$, generated by the direct product of $`p=2\mathrm{}10`$, $`q=0\mathrm{}5`$ and $`n=0\mathrm{}4`$. The QP algorithm is used to compute the coefficients $`A_{pqn}^\pm `$ that optimize the fit to the photometric data and the mean velocity and dispersion profiles. This is done by minimizing a $`\chi ^2`$-like variable of the form

$$\chi ^2=\underset{i}{}w_i\left(\frac{\mu _{\mathrm{obs},i}\mu _{\mathrm{mod},i}}{\mu _{\mathrm{mod},i}}\right)^2$$
(12)   
where the $`\mu `$’s represent the photometric and kinematic data collectively and the sum is over all data points. While complete two-dimensional surface photometry is available, we choose for convenience to fit the photometry only at the points where there are kinematic measurements. This allows some of the models to come out a bit peanut-shaped, but the effect on the inferred mass profile is minimal. The weights $`w_i`$ are in principle arbitrary, but are here assigned the value 1.

#### 3.2.2 Results

The best constant-$`M/L`$ 3-integral model, at an inclination of $`90\mathrm{°}`$, is shown in figure 8. The global value of $`M/L_B=2.8h`$ has been chosen to optimize the overall fit, i.e., without assigning larger weight to any particular region, though the higher $`S/N`$ at smaller radii makes the model appear pinned to the central dispersion. The shallow slope of the RMS velocity profile cannot be reproduced; thus, allowing for anisotropy does not eliminate the need for an $`M/L`$ gradient. In fact, the fit is even slightly worse than for the 2-integral moment model of § 3.1.2. We ascribe this to the positivity requirement on $`f`$ in the QP models, which is not applied in the moment models. If the positivity constraint is lifted, the same QP components can produce a fit with a $`\chi ^2`$ a factor of 3 lower. We take this as an indication that the constant-$`M/L`$ moment model actually does not correpsond to a nonnegative distribution function.

Allowing for an $`M/L`$ gradient improves the fit substantially. The dashed lines in figure 9 show a model with $`i=90\mathrm{°}`$, $`p=1`$, $`B=1.27h^1\mathrm{kpc}^1`$, and $`\mathrm{{\rm Y}}_0=5.1h`$ ($`B`$ band, in solar units), compared with the observed rotation (squares) and dispersion (crosses) profiles. The fit to the kinematics alone, ignoring the photometry, is passable, with $`\chi ^2/\nu =2.03`$ for $`\nu =108`$ degrees of freedom. The discrepancy at large radii on PA 315 is not especially worrisome since this is where isophotal and kinematic twisting sets in, and axisymmetric models should not be expected to reproduce the detailed velocity field.

The local and cumulative $`M/L`$ profiles for this model are shown as the dashed lines in figure 10. The central values agree very well with those shown in Figure 7 for the 2-integral model, and appear to be insensitive to assumptions about the functional form (i.e., equation  vs. equation .) At large $`r`$ the local $`M/L`$ gradient is roughly a factor $`1.5`$ smaller than that obtained from the 2-integral model. This translates to a $`20\%`$ difference in the cumulative $`M/L`$ gradient because of the steepness of the mass profile.

Obtaining a strict lower limit on the $`M/L`$ gradient, while desirable, is not especially practical. A shallower gradient can be obtained at the expense of a more pathological distribution function at low binding energy, corresponding to distortions of the LOSVD at radii for which there are no data. As there is no objective criterion for accepting or rejecting such models, we have rather subjectively chosen the parameters $`B`$ and $`p`$ so as not to overly tax the QP algorithm; thus the distribution function is not on the verge of going negative, and the anisotropy is not extreme. The latter property is illustrated in the top two panels of figure 11, which show the ratios of the unprojected dispersions $`\sigma _R`$, $`\sigma _\varphi `$, and $`\sigma _z`$ in the meridional plane. The velocity ellipsoids are generally elongated in the tangential direction, and $`\sigma _z/\sigma _R`$ is within $`10\%`$ of unity nearly everywhere. The model becomes more tangential at larger $`R`$; but remember that the kinematic data end near $`70\mathrm{}`$, so the velocity distribution beyond this point is not well constrained.

Though we have not checked explicitly, it seems virtually certain that for any inclination $`i>50\mathrm{°}`$ (for which $`c>0.35`$) it is possible to find many values of $`B`$ and $`p`$ that produce an acceptable fit. For the shape fits of Section 2, the posterior densities in figures 1 and 3 give a most-probable intrinsic flattening $`c`$ between $`0.55`$ and $`0.65`$, which would correspond, for an axisymmetric system, to $`i60\mathrm{°}`$. A QP model at this inclination, using the same $`B`$ and $`p`$ as the $`i=90\mathrm{°}`$ model above, is shown as the dotted lines in figure 9. The formal best fit has been scaled up by $`3\%`$ in velocity, or $`6\%`$ in mass, to improve the fit to the central dispersion. The local and cumulative $`M/L`$ profiles are plotted as the dotted lines in figure 10, and the ratios of the principal dispersions are shown in the bottom two panels in figure 11. This rather flat model is tangential at all radii, and at large $`R`$ (past where the data end) becomes radially quite cold. Such a configuration is not at all impossible, especially if, as SSC suggest, the main body of the galaxy is surrounded by a differentially precessing disk acquired through a merger or other accretion event.

## 4 Discussion

By modeling its velocity field and surface photometry, we have found that NGC 1700 is nearly axisymmetric and oblate. We can say with better than $`95\%`$ confidence that if it were viewed down its short axis, it would appear as at most an E1 at radii $`36\mathrm{}`$ ($`6.7h^1\mathrm{kpc}`$). Beyond this point, if the density surfaces remain aligned, then photometric and kinematic twists indicate that it becomes significantly flatter and more triaxial.

The dominant source of systematic error in the determination of intrinsic shape is our lack of prior knowledge about the internal dynamical structure, specifically, how we should expect the mean streaming velocity to vary away from the symmetry planes. Axisymmetry is strongly favored if the velocity declines away from the $`xy`$ plane; triaxiality is preferred if it doesn’t. It is therefore important to understand how the latitude dependence of the mean streaming may be determined by physics. Preliminary results of an analysis of group merger simulations (Dutta, Statler, & Weil 1998), though limited by small-number statistics, suggest that the mean velocity does decline away from the $`xy`$ plane in objects produced by this mechanism. It remains to be seen whether different formation scenarios would produce different behavior.

NGC 1700 falls into the family of steep-cusped ellipticals delineated by WFPC photometry, standing out, in fact, as its most luminous member (so far). In general, the steep-cusped galaxies show characteristics one would expect of axisymmetric systems, such as disky isophotes and relatively rapid rotation. The near axisymmetry we have found is consistent with Merritt’s (1997) result that orbital stochasticity prevents the existence of self-consistent triaxial equilibria for the steep ($`r^2`$) cusped Jaffe profile unless $`T0.4`$, $`T0.9`$, or $`c0.8`$. It is tempting to interpret the change of shape beyond $`40\mathrm{}`$ as marking the limit out to which triaxiality has been eaten away by stochastic scattering of box orbits. But it is far from clear that the cusp in NGC 1700 is the cause of this structure. As discussed by SSC, there is reason to believe that the observed photometric and kinematic twists arise from an outer accreted stellar disk or ring which is incompletely phase-mixed and therefore intrinsically misaligned with the rest of the galaxy. In this case the dynamical time at the location of the twist, times a factor of order a few, indicates the time since the last significant merging event. The dynamical age of $`2`$ to $`4h^1\mathrm{Gyr}`$ (SSC) is corroborated by Whitmore et al. (1997), who have observed NGC 1700’s globular clusters and find colors and magnitudes consistent with a comparably old cluster population.

We have also found, using axisymmetric 2-integral moment models and 3-integral quadratic programming models, that NGC 1700 must have a radially increasing mass to light ratio. No constant-$`M/L`$ model can fit the data. The cumulative (enclosed) $`M/L_B`$ increases from a central value of $`5h`$ to $`10h`$ (in solar units) at $`r=12h^1\mathrm{kpc}`$. The case for dark matter in NGC 1700 is strong; even allowing for varying anisotropy, the best constant-$`M/L`$ model fits the RMS velocity profiles at no better than $`\chi ^2/\nu =5.74`$ (for $`\nu =50`$). In this sense the evidence of dark matter is comparable to that in NGC 2434 (Rix et al. 1997); in fact, NGC 1700 is arguably a more compelling case since only spherical models have been applied to NGC 2434. Other systems suspected on stellar dynamical grounds of harboring substantial dark halos, such as NGC 1399, 4374, 4472, 5813, 7626, 7796, and IC 4296 (Saglia, Bertin, & Stiavelli 1992, Saglia et al. 1993, Bertin et al. 1994) are statistically less convincing. NGC 1700 appears to represent the strongest stellar dynamical evidence to date for dark halos in ellipticals.

Bertin et al. (1994) point out that, of the ellipticals listed above with stellar dynamical indications of dark matter, most turn out to be known bright X-ray sources. Unfortunately, NGC 1700 seems never to have been observed by an X-ray telescope, though ROSAT programs are pending. NGC 1700 bears a more-than-passing resemblance to NGC 7626 (Saglia, Bertin, & Stiavelli 1992), having a very similar luminosity, central dispersion, effective radius, and distance. However, NGC 7626 is 2 orders of magnitude more powerful in the radio (Bender et al. 1989), and also resides in a rich cluster while NGC 1700 has only one nearby companion and lies among small groups. It will be very interesting to see whether these two apparently similar ellipticals, in very different environments, turn out to be X-ray twins.

###### Acknowledgements.

TSS acknowledges support of this work from NASA Astrophysical Theory Program Grant NAG5-3050 and NSF CAREER grant AST-9703036. We thank the director and staff of the Multiple Mirror Telescope Observatory for their generous assistance and allocations of time to this project, and Trevor Ponman and Joe Shields for timely information.    

## Appendix A: Combining Shape Estimates for Multiple Annuli

We begin by introducing the notation $`P_{ijk\mathrm{}}(T_i,c_i,\mathrm{\Omega }_i)`$ for the marginal probability density for the shape and orientation of annulus $`i`$, making use of the data from annuli $`i,j,k,\mathrm{}`$. Suppose that we have obtained $`P_1(T_1,c_1,\mathrm{\Omega }_1)`$ and $`P_2(T_2,c_2,\mathrm{\Omega }_2)`$ for annuli 1 and 2 individually (cf. equation ). (Here both $`T`$ and $`c`$ refer to the light distribution; we omit the $`L`$ subscripts for clarity.) We assume that the luminosity distribution is intrinsically aligned ($`\mathrm{\Omega }_1=\mathrm{\Omega }_2=\mathrm{\Omega }`$) and require that the observed position angle difference $`\delta _{12}\pm \sigma _{\delta _{12}}`$ be reproduced. Let $`L_{\delta _{12}}(T_1,T_2,\mathrm{\Omega })`$ be the likelihood of the observed twist given the pair of intrinsic shapes $`(T_1,c_1)`$ and $`(T_2,c_2)`$, and orientation $`\mathrm{\Omega }`$:

$$L_{\delta _{12}}(T_1,T_2,\mathrm{\Omega })=\mathrm{exp}\left\{\frac{[\delta _{12}\delta _{12}^{\mathrm{mod}}(T_1,T_2,\mathrm{\Omega })]^2}{2\sigma _{\delta _{12}}^2}\right\}.$$
(13)   
The axis ratios $`c_1`$ and $`c_2`$ do not appear in either the likelihood $`L_{\delta _{12}}`$ or the predicted twist $`\delta _{12}^{\mathrm{mod}}`$ because in the quasi-local approximation the major axis position angle in annulus $`i`$,

$$\mathrm{PA}_i=\frac{1}{2}\mathrm{tan}^1\frac{T_i\mathrm{cos}\theta _E\mathrm{sin}2\varphi _E}{|T_i(\mathrm{sin}^2\varphi _E\mathrm{cos}^2\varphi _E\mathrm{cos}^2\theta _E)\mathrm{sin}^2\varphi _E|}.$$
(14)   
is independent of $`c_i`$. The joint probability that the two shapes fit the local photometry and kinematics and the PA twist is given by

$$P(T_1,c_1,T_2,c_2)=𝑑\mathrm{\Omega }P_1(T_1,c_1,\mathrm{\Omega })P_2(T_2,c_2,\mathrm{\Omega })L_{\delta _{12}}(T_1,T_2,\mathrm{\Omega }),$$
(15)   
but this is a somewhat unwieldy quantity. Of more use is the marginal probability for, say, annulus 1 alone:

$$P_{12}(T_1,c_1)=𝑑T_2𝑑c_2P(T_1,c_1,T_2,c_2)=𝑑\mathrm{\Omega }P_1(T_1,c_1,\mathrm{\Omega })K_{12}(T_1,\mathrm{\Omega }),$$
(16)   
where the “twist kernel” $`K_{12}`$ is given by

$$K_{12}(T_1,\mathrm{\Omega })=𝑑T_2L_{\delta _{12}}(T_1,T_2,\mathrm{\Omega })P_2(T_2,\mathrm{\Omega }),$$
(17)   
and $`P_2(T_2,\mathrm{\Omega })`$ is the marginal distribution integrated over $`c_2`$, i.e.,

$$P_2(T_2,\mathrm{\Omega })𝑑c_2P_2(T_2,c_2,\mathrm{\Omega }).$$
(18)   
The posterior $`P_{12}(T_1,c_1)`$ is the shape estimate for annulus 1 considering all of the data in both annuli, but regardless of the actual shape of annulus 2.

In our case the $`P_i(T_i,c_i,\mathrm{\Omega }_i)`$ are obtained from the dynamical models, so the $`T_i`$ are mass triaxialities. If one is unwilling to assume a relation between the triaxialities of mass and luminosity, one can still combine the annular fits, assuming only that the principal axes are intrinsically aligned and foregoing the isophotal twist constraint, by letting $`\sigma _{\delta _{12}}\mathrm{}`$, which gives $`L_{\delta _{12}}=1`$ by equation (13).

The generalization to $`N`$ annuli uses the measured PA gradients $`\delta _{12},\delta _{23},\mathrm{},\delta _{N1,N}`$, as follows:

$$P_{12\mathrm{}N}(T_i,c_i,\mathrm{\Omega })=P_i(T_i,c_i,\mathrm{\Omega })K_{i,i1,\mathrm{},1}^{\mathrm{int}}(T_i,\mathrm{\Omega })K_{i,i+1,\mathrm{},N}^{\mathrm{ext}}(T_i,\mathrm{\Omega }),$$
(19)   
where the interior and exterior kernels are given by

$`K_{i,i1,\mathrm{},1}^{\mathrm{int}}(T_i,\mathrm{\Omega })`$ $`=`$ $`{\displaystyle 𝑑T_{i1}P_{i1}(T_{i1},\mathrm{\Omega })L_{\delta _{i1,i}}(T_{i1},T_i,\mathrm{\Omega })}`$   
$`{\displaystyle 𝑑T_{i2}P_{i2}(T_{i2},\mathrm{\Omega })L_{\delta _{i2,i1}}(T_{i2},T_{i1},\mathrm{\Omega })\mathrm{}}`$   
$`{\displaystyle 𝑑T_1P_1(T_1,\mathrm{\Omega })L_{\delta _{1,2}}(T_1,T_2,\mathrm{\Omega })};`$   
$`K_{i,i+1,\mathrm{},N}^{\mathrm{ext}}(T_i,\mathrm{\Omega })`$ $`=`$ $`{\displaystyle 𝑑T_{i+1}P_{i+1}(T_{i+1},\mathrm{\Omega })L_{\delta _{i,i+1}}(T_i,T_{i+1},\mathrm{\Omega })}`$   
$`{\displaystyle 𝑑T_{i+2}P_{i+2}(T_{i+2},\mathrm{\Omega })L_{\delta _{i+1,i+2}}(T_{i+1},T_{i+2},\mathrm{\Omega })\mathrm{}}`$   
$`{\displaystyle 𝑑T_NP_N(T_N,\mathrm{\Omega })L_{\delta _{N1,N}}(T_{N1},T_N,\mathrm{\Omega })}.`$   
As before, the shape estimate for annulus $`i`$, using all of the data but regardless of the shape at other radii, is the marginal distribution $`P_{12\mathrm{}N}(T_i,c_i)=𝑑\mathrm{\Omega }P_{12\mathrm{}N}(T_i,c_i,\mathrm{\Omega })`$. The construction of equations (19) and (20) guarantees that the marginal distribution in orientation, $`P_{12\mathrm{}N}(\mathrm{\Omega })=𝑑T_i𝑑c_iP_{12\mathrm{}N}(T_i,c_i,\mathrm{\Omega })`$, is independent of $`i`$, as it must be given the assumption of intrinsic alignment.

## Appendix B: The Contrast and $`v^{}(t)`$ Functions

In the models used in Section 2, properties of the phase space distribution function are subsumed into the “similar flow” ansatz, a scalar constant $`C`$, and a function of one variable $`v^{}(t)`$. The last two describe the mean velocity across the $`xz`$ plane on one fiducial spherical or ellipsoidal shell, which in turn determines the velocity field over the whole shell once $`T`$ and the luminosity density are specified.

The “contrast” $`C`$ is defined as the ratio of the $`y`$ component of the mean velocity on the $`x`$ axis to that on the $`z`$ axis, on the fiducial shell. One may think of this as the ratio of the maximum mean streaming velocities produced by short-axis and long-axis tubes. Models lacking short-axis or long-axis tube streaming have $`C=0`$ or $`C=\mathrm{}`$, respectively. Models are run with values of $`C`$ both independent of intrinsic shape and given by the four functions of $`(T,c)`$ shown in Figure 2 of S94c.

The function $`v^{}(t)`$ gives the angular dependence of the mean velocity across the $`xz`$ plane on the fiducial shell. The variable $`t`$ is a rescaled polar angle, defined so that $`t=0`$ corresponds to the $`x`$ axis, $`t=2`$ to the $`z`$ axis, and $`t=1`$ to the locus dividing the short-axis from the long-axis tubes. For spherical shells,

$$t=\{\begin{array}{cc}2\frac{\mathrm{sin}^2\theta }{T},\hfill & \theta <\mathrm{sin}^1\sqrt{T},\hfill \\ \frac{\mathrm{cos}^2\theta }{1T},\hfill & \theta >\mathrm{sin}^1\sqrt{T};\hfill \end{array}$$
(21)   
for ellipsoidal shells the relation between $`t`$ and $`\theta `$ is more complicated, and we refer the reader to Section 3.1 of S94c. By definition $`v^{}(0)=C`$ and $`v^{}(2)=1`$. In the intervals $`[0,1)`$ and $`(1,2]`$, $`v^{}(t)`$is assumed to be either piecewise-constant or piecewise linear, in the latter case falling to zero at $`t=1`$. This gives four possible forms for $`v^{}(t)`$: Type 1, constant in both intervals; Type 2, linear in both intervals; Type 3, linear in $`[0,1)`$ and constant in $`(1,2]`$; and Type 4, constant in $`[0,1)`$ and linear in $`(1,2]`$.