text:
# On Solar Radius Variation with magnetic Activity

## 1 Introduction

By analyzing the f-mode frequency data from MDI for the period 1996–2001, Dziembowski et al. (2001) claimed that the Sun shrinks in radius from the solar minimum to the solar maximum, at the rate of about 1.5 km per year. If correct, then this will certainly be a most dramatic discovery in solar physics. However, this claim has not yet been supported by other researchers, who often present conflicting results. Antia (2003) provides a review of the subject and concludes that there is no compelling evidence for solar radius change. It may be noted that some earlier studies (Nöel, 1997; Emilio et al., 2000) reported much larger radius variation (up to 700 km) based on other (non-helioseismic) methods. It is clear that more work is needed before a definitive conclusion can be reached. The aim of this paper is to discuss whether on theoretical grounds we expect a solar radius variation with solar cycle. Other than some early work by Spruit (1982, 1994), we have not been able to find any other theoretical analysis of this problem. In the light of recent developments of solar physics, this theoretical question certainly needs revisiting.

There are two ways in which a magnetic field present in the solar interior may lead to an increase of solar radius.

1. The pressure of the magnetic field embedded in the gas may cause the gas to expand, thereby inflating the size of the Sun.
2. The inhibition of convection by the magnetic field may cause a piling up of heat below the magnetic field. This would result in heating of the surrounding regions, leading to an expansion.

From both of these arguments, it appears that the radius of the Sun should increase when more magnetic fields are present. Since we normally think of the solar maximum as being the time when the Sun has more magnetic fields in the interior, we would expect the Sun to be larger in size during the solar maximum. Shrinking of the Sun with increasing magnetic activity, therefore, would seem improbable at the first sight.

We shall argue that the recent solar dynamo model of Nandy and Choudhuri (2002), in principle, allows for the possibility of shrinking of the Sun during the solar maximum. The important theoretical question is whether the theoretically expected radius decrease agrees with the claim of Dziembowski et al. (2001). Spruit (1982, 1994) studied the effect of heat blocking by magnetic fields and concluded that the total radius variation cannot be more than 0.14 km. We mainly look at the other effect of magnetic pressure and present order-of-magnitude estimates, suggesting that the radius variation, for reasonable assumptions, cannot be as large as what Dziembowski et al. (2001). We have also made a mathematical formulation of the problem, to go beyond order-of-magnitude estimates. To our surprise, we found that the mathematical problem is much more difficult to solve than what we anticipated. Although our mathematical analysis so far has not yielded very definite results, we present a brief account of it to help future researchers who may wish to analyze this problem in greater depth.

Section 2 addresses the question whether we expect the radius to increase or decrease during the solar maximum. An order-of-magnitude estimate is presented in Section 3, showing that a much smaller variation of radius is expected than what is reported by Dziembowski et al. (2001). Our incomplete mathematical analysis is given in Section 4. Finally, our conclusions are summarized in Section 5.

## 2 Increase or decrease of radius with magnetic activity?

Both the effects listed above—expansion due to magnetic pressure and heat blocking by inhibition of convection—would lead to an increase of the solar radius when there is more magnetic flux. So it appears that the Sun can shrink with increasing magnetic activity only if there is less magnetic flux inside the Sun at the time of solar maximum. We know that the Sun has much more magnetic flux near the surface in the form of sunspots during the solar maximum. It used to be tacitly assumed that the magnetic flux in the interior of the Sun should also have a peak during the solar maximum. Some recent dynamo models, however, raise doubts about this assumption.

It has been realized for some time that the meridional circulation plays an important role in the solar dynamo (Choudhuri et al., 1995). This meridional circulation is poleward near the surface and must be equatorward at the bottom of the solar convection zone (SCZ). With the mapping of angular velocity distribution in the solar interior by helioseismology, it is clear that the strongest gradient of angular velocity occurs at high latitudes within the tachocline and that is the most likely region for the creation of toroidal field. To explain why active regions appear at low latitudes, Nandy and Choudhuri (2002) have developed a solar dynamo model in which the meridional flow penetrates a little bit below the tachocline and advects the strong toroidal field (created at high latitudes inside the tachocline) to low latitudes through convectively stable layers below the tachocline. Only when the meridional circulation rises in the low latitudes, the toroidal fields are brought within the SCZ and rise quickly due to magnetic buoyancy to form sunspots. In the model of Nandy and Choudhuri (2002), therefore, the toroidal flux which erupts at the solar surface was created a few years earlier at high latitudes within the tachocline, since it took time for the flux to be transported from there. If the meridional flow speed at the base of the SCZ is taken to be of order 1.2 m s-1 (Hathaway et al., 2003), then the time taken for flux transport should be in the range 5-10 yr. Hence the flux should have been created a few years before the solar maximum. In this scenario, the solar maximum is the time when the Sun merely gets rid of the flux stored below the SCZ. In contrast to the traditional dynamo models in which during the solar maximum the Sun generates the strong toroidal field and the interior flux has a peak value, the model of Nandy and Choudhuri (2002) suggests the solar maximum to be the time when the magnetic flux in the Sun is reduced after being created a few years earlier. Although the solar maximum is the period when the magnetic flux near the surface in the form of sunspots peaks, it may also be a time when the magnetic flux at the bottom of the SCZ actually decreases. So, to settle the question whether a magnetic effect would lead to increase of decrease of the solar radius, we need to ascertain whether the magnetic effect is caused by the magnetic flux near the surface or the magnetic flux at the bottom of the SCZ. We point out that several simulations of active region formation (Choudhuri and Gilman, 1987; Choudhuri, 1989; D’Silva and Choudhuri, 1993; Fan et al., 1993; D’Silva and Howard, 1993; Caligari et al., 1995) have pinned down the value of the magnetic field at the bottom of SCZ to be around $`10^5`$ G.

We first make a few remarks on Spruit’s analysis of the heat blocking problem (Spruit, 1982, 1994). If the magnetic field is sufficiently strong to suppress convection, the effects of heat blocking become independent of the magnetic field (i.e. it does not matter if the magnetic is just strong enough to suppress the convection or much stronger than that). On the other hand, results of radius variation due to magnetic pressure depend critically on the value of the magnetic field, as we shall see below. Spruit (1982) mainly considered heat blocking by vertical flux tubes near the solar surface. Although we believe in the existence of strong toroidal flux tubes with $`10^5`$ G magnetic field at the bottom of SCZ, heat blocking by these flux tubes will not have much effect for two reasons. Firstly, due to the high concentration of magnetic field inside these flux tubes, the flux tubes would occupy a smaller volume compared to the volumes of flux tubes below the sunspots. Hence convection would be suppressed in much smaller volumes. Secondly, the bottom of the SCZ is the region where convective heat transport just begins to win over the radiative transport. A suppression of convection there is not expected to have as much effect as it would have in the interior of SCZ. We thus conclude that we can still trust Spruit’s results and his conclusion that the radius variation can at most be of order $`2\times 10^7R_{}`$, which is only $`0.14`$ km. It is not only considerably less than what is reported by Dziembowski et al. (2001), it also has the wrong sign. More heat blocking below sunspots takes place during the solar maximum and the radius should increase at that time.

If magnetic flux decreases at the bottom of SCZ during the solar maximum, as suggested by the model of Nandy and Choudhuri (2002), then the magnetic pressure is reduced and it is possible that the Sun shrinks due to this at the time of the solar maximum. We now turn to making an estimate of how much this radius decrease is expected to be.

## 3 Basic considerations and rough estimates

We now consider radius variations due to magnetic pressure of the concentrated flux tubes at the bottom of SCZ. When the magnetic flux is created, it disturbs the hydrostatic equilibrium. It is well known that the time scale for restoration of hydrostatic equilibrium is very short (see, for example, Schwarzschild, 1958, p. 32). As Schwarzschild (1958) comments: “There is no doubt that punishment would follow swiftly for any disobedience of the hydrostatic law.” We, therefore, assume that the Sun would fairly quickly adjust itself to a new hydrostatic equilibrium when the magnetic field is created. So we shall basically compare hydrostatic equilibria with and without magnetic fields embedded inside, without bothering about the adjustment time scales. Our approach thus is going to be quite different from the approach of Spruit (1982) while studying the heat blocking problem, in which two time scales (the diffusion time and the Kelvin-Helmholtz time) play important roles.

To calculate the expansion of a flux tube due to magnetic pressure, we need to specify the thermodynamic conditions. If the temperature in the interior of the flux tube is the same as the outside temperature, then certainly the flux tube expands due to magnetic pressure (see, for example, Choudhuri, 1998, §14.7). However, such a flux tube also tends to be buoyant and to rise against gravity. The flux tubes are believed to be stored in the mildly sub-adiabatic region immediately below the bottom of the SCZ. If a buoyant flux tube rises in a region of sub-adiabatic temperature gradient, it can be easily shown that it becomes less buoyant (Parker, 1979, §8.8; Moreno-Insertis, 1983). One expects that the flux tube would rise until it becomes neutrally buoyant, i.e. until the interior temperature decreases sufficiently to make the density inside and outside exactly equal. If the flux tube is neutrally buoyant, then its density is no different from the outside density and there is no increase in the overall size of the system. Thus, if all the magnetic flux inside the Sun remained neutrally buoyant, then there would be no change in the solar radius at all. In view of our lack of understanding of the physics of flux tubes at the bottom of SCZ, let us allow for the possibility that this may not be the case. The penetrating meridional circulation proposed in the model of Nandy and Choudhuri (2002) would drag the flux tubes downward, making them hotter and not allowing them to come up by magnetic buoyancy. Considering a favourable circumstance that the flux tubes are not colder than the surroundings, let us see if we can have sufficient radius increase of the Sun.

We now make a rough estimate by considering a one-dimensional atmosphere in rectangular geometry. Suppose a horizontal magnetic field $`B`$ is created within this atmosphere in a layer between $`z`$ and $`z+L`$. Because of the magnetic pressure inside this layer, the gas pressure would be less than the pressure $`p`$ that would exist there if the magnetic field were not present. To find the corresponding density decrease, we assume the region with the magnetic field to have the same temperature as the surrounding region. Then we can give the same arguments which we give for flux tubes in thermal equilibrium with surroundings (see, for example, Choudhuri, 1998, §14.7) and conclude that the density decrease is given by

$$\frac{\delta \rho }{\rho }\frac{B^2}{8\pi p}.$$
$`(1)`$
In a one-dimensional atmosphere, such a density decrease implies that the magnetic layer would be inflated by an amount $`\delta L`$ compared to thickness which the gas in this layer would have in the absence of the magnetic field. Hence, when the magnetic field is created, the overlying atmosphere will be pushed up by this amount $`\delta L`$, which is given by

$$\frac{\delta L}{L}\frac{\delta \rho }{\rho }.$$

Note that we are not bothering about the signs of $`\delta \rho `$ and $`\delta L`$. Then from (1),

$$\frac{\delta L}{L}\frac{B^2}{8\pi p}.$$
$`(2)`$
Since

$$F=BL$$

is the magnetic flux through the magnetic layer (per unit length of the layer measured in a direction perpendicular to $`B`$ in the layer), we can write (2) as

$$\delta L\frac{F^2}{8\pi pL}.$$
$`(3)`$
It is clear from (3) that the rise of the atmosphere is more if the same flux $`F`$ is concentrated in a narrower layer (i.e. if the magnetic field is stronger).

Taking $`p=5\times 10^{13}`$ dyn cm-2 at the bottom of the SCZ, (2) tells us that for magnetic field created there we must have

$$\delta L7.7\times 10^6\left(\frac{B}{10^5}\right)^2L.$$
$`(4)`$
For magnetic field of order $`10^5`$ G, we expect $`\delta L/L`$ to be only of order $`10^5`$.

We now estimate the radius shrinkage of the Sun due to the release of magnetic flux from the interior during the solar maximum. Assume that there are about $`10^4`$ eruptions during the maximum and each eruption brings out flux in the range from $`10^{20}`$ Mx to $`10^{21}`$ Mx. Then the total flux coming out of the interior of the Sun is about $`10^{24}`$$`10^{25}`$ Mx. This flux comes from the bottom of SCZ, where the magnetic field is $`10^5`$ G. This gives a cross-sectional area of $`10^{19}`$$`10^{20}`$ cm2. To compare with the one-dimensional model, we distribute this flux in a uniform shell at the bottom of SCZ. Taking $`L`$ to be the thickness of this shell, the cross-sectional area through which the magnetic field passes should be equal to $`\pi R_{\mathrm{bot}}L`$, where $`R_{\mathrm{bot}}`$ is the radial distance of the bottom of SCZ. Equating this cross-sectional area to $`10^{19}`$$`10^{20}`$ cm2, we get $`L`$ in the range from about $`10^8`$ cm to $`10^9`$ cm (i.e. from about 1000 km to 10,000 km). It follows from (4) that the expected maximum shrinkage of the Sun can only be in the range from 0.1 km to 1 km. If this radius shrinkage takes place over 5 years, then even the most favourable value 1 km of radius shrinkage will give a shrinkage rate of only 0.2 km per year. According to (3), the magnetic field at the bottom of SCZ has to be squeezed in a layer narrower by one order of magnitude if the radius shrinkage rate had to be as large as what Dziembowski et al. (2001) find. This would imply the horizontal magnetic field at the bottom of SCZ should be of order $`10^6`$ G. Such a strong field would create problems with the theoretical explanation of Joy’s law (D’Silva and Choudhuri 1993).

It should be kept in mind that the radius shrinkage rate estimated above is an over-estimate. If the magnetic flux tends to become neutrally buoyant in some regions, then certainly the radius variation will be less. Additionally, Choudhuri (2003) has argued that magnetic flux would exist in concentrated form only in limited regions at the bottom of SCZ, the field being more diffuse in other regions. Since diffuse magnetic field would not cause much radius variation, in accordance with (3), the actual radius variation would be smaller if concentrated magnetic field regions are intermittent.

## 4 A hydrostatic model with a magnetic layer

After the order-of-magnitude estimate, we now try to construct hydrostatic models of the SCZ, with and without magnetic fields in the interior, to determine the distance by which the outer surface shifts between these two models. After considering a hydrostatic model without magnetic fields, we shall assume that the magnetic field is created in a shell $`r_b<r<r_t`$ and calculate the new equilibrium. To our great surprise, we found this problem to be much more complicated than what we expected. Although we have not been able obtain very definite results, we discuss the formulation of the problem to provide guidance for future researchers who may wish to tackle the problem with more realistic models.

Let us first consider the case without the magnetic field and construct a hydrostatic model from the bottom of the overshoot layer at $`R_{BCZ}=0.67R_{}`$ to the surface at $`R_{}`$. Let $`p_0(r),\rho _0(r),T_0(r)`$ be the pressure, density and temperature structures in this case. The equation for hydrostatic equilibrium is given by

$$\frac{dp_0(r)}{dr}=\frac{\rho _0(r)GM_{}}{r^2}.$$
(5)
The temperature gradient is slightly sub-adiabatic in the overshoot layer and slightly super-adiabatic in the SCZ. We assume an adiabatic stratification throughout and use the equation of state

$$p_0(r)=\kappa \rho _0(r)^\gamma .$$
(6)
Substituting (6) in (5) and solving the resulting ODE with the boundary condition $`\rho _0(R_{})=0`$, we have the following expressions for $`\rho _0(r),p_0(r)`$, and $`T_0(r)`$. To compute the temperature we have used the ideal gas law $`p=\rho R_{gas}T`$.

$`\left(\begin{array}{c}\rho _0(r)\\ p_0(r)\\ T_0(r)\end{array}\right)=\left(\begin{array}{c}[\frac{\gamma 1}{\gamma }\frac{GM}{\kappa }\{\frac{1}{r}\frac{1}{R_{}}\}]^{\frac{1}{\gamma 1}}\\ \kappa [\frac{\gamma 1}{\gamma }\frac{GM}{\kappa }\{\frac{1}{r}\frac{1}{R_{}}\}]^{\frac{\gamma }{\gamma 1}}\\ \frac{\gamma 1}{\gamma }GM\{\frac{1}{r}\frac{1}{R_{}}\}\end{array}\right)`$ (7)
We now consider that a magnetic field $`B(r)`$ is created in the shell $`r_b<r<r_t`$, i.e. the magnetic field remains zero in the layers below and above. The hydrostatic equation is satisfied piecewise in all the three layers. The subscript ’$`1`$’ indicates quantities below the magnetic shell, the subscript ’$`2`$’ denotes quantities inside the magnetic layer, whereas ’$`3`$’ indicates those above. In the field-free regions ’1’ and ’3’, pressure and density should satisfy the hydrostatic equation (5) and the adiabatic equation (6) (with the subscripts ’1’ and ’3’ replacing the subscript ’0’). Then the solutions also will be similar to (7). We can assume the layer below to remain unchanged after the introduction of the magnetic field. So $`\rho _1(r)`$, $`p_1(r)`$ and $`T_1(r)`$ are equal to $`\rho _0(r)`$, $`p_0(r)`$ and $`T_0(r)`$ below $`r=r_b`$. In the case of the upper layer, we assume that outer radius is now $`R`$ instead of $`R_{}`$. So $`\rho _3(r)`$, $`p_3(r)`$ and $`T_3(r)`$ can simply be obtained from (7) by replacing $`R_{}`$ by $`R`$ (which essentially means changing the constant of integration), i.e.

$$\rho _3(r)=\left[\frac{(\gamma 1)GM_{}}{\gamma \kappa }\left(\frac{1}{r}\frac{1}{R}\right)\right]^{\frac{1}{\gamma 1}}$$
(8)
with similar corresponding expressions for $`p_3(r)`$ and $`T_3(r)`$.

Let us now focus our attention in the magnetic layer $`r_b<r<r_t`$. The hydrostatic equation there is

$$\frac{d}{dr}\left[p_2(r)+\frac{B^2(r)}{8\pi }\right]=\frac{\rho _2(r)GM_{}}{r^2}.$$
(9)
Actually, there should also be a magnetic tension term which would break the spherical symmetry. Apart from breaking the spherical symmetry, the other problem which such a term introduces is that it cannot be balanced by a combination of gradient forces, like the pressure gradient and gravity (which comes from the gradient of gravitational potential). For flux tubes at the bottom of SCZ, the magnetic tension term is supposed to be negligible compared to other term (the gradient of magnetic pressure becomes much larger near the boundaries of flux tubes). Hence we have neglected this troublesome tension term. We further assume that the pressure and density inside the magnetic layer satisfy an adiabatic relation

$$p_2(r)=\kappa _2\rho _2(r)^\gamma $$
(10)
with a constant $`\kappa _2`$ which can be different from $`\kappa `$. Since we do not expect vigorous convection inside the magnetized region, the adiabatic assumption there is somewhat questionable even though it simplifies our life.

To find the structure of the magnetic layer, we need to solve (9) and (10) subject to certain conditions. Since we want the pressure to be continuous across the top and bottom of the magnetic layer, we should have

$$p_1(r_b)p_2(r_b)=\frac{B^2(r_b)}{8\pi },$$
(11)
$$p_3(r_t)p_2(r_t)=\frac{B^2(r_t)}{8\pi }.$$
(12)
Additionally, total mass of the system does not change while the magnetic field is being created, implying

$`{\displaystyle _{R_{BCZ}}^R_{}}\rho _0(r)r^2𝑑r`$ $`=`$ $`{\displaystyle _{R_{BCZ}}^{r_b}}\rho _1(r)r^2𝑑r+{\displaystyle _{r_b}^{r_t}}\rho _2(r)r^2𝑑r+{\displaystyle _{r_t}^R}\rho _3(r)r^2𝑑r`$ (13)
$`=`$ $`{\displaystyle \frac{0.02M_{}}{4\pi }},`$
where we have assumed that only a total $`0.02M_{}`$ amount of mass is contained in the outer regions of the Sun we are considering, justifying our approximating the gravitational field by an inverse square law in (5). It easily follows from (7) and (13) that

$$\kappa =\frac{(\gamma 1)GM_{}}{\gamma }\left(\frac{4\pi }{0.02M_{}}\right)^{\gamma 1}R_{}^{3\gamma 4}\left(_{1/a}^1\xi ^{\frac{2\gamma 3}{\gamma 1}}(1\xi )^{\frac{1}{\gamma 1}}𝑑\xi \right)^{\gamma 1}$$
(14)
where $`a=R_{}/R_{BCZ}`$ and $`\xi `$ is a dummy non-dimensional variable. We would also like to specify a thermal condition of the magnetic layer, i.e. to specify that its temperature is not very different from the temperatures of the surroundings. However, for the cases we considered, we do not seem to have this freedom.

To solve the problem, we need to specify the magnetic field in the flux layer. The simplest possibility is to assume that we have a magnetic field of constant $`B`$ inside this layer. Then (9) reduces to the form (5) and can easily be integrated in combination with (10) to give

$$\rho _2(r)=\left[\frac{(\gamma 1)GM_{}}{\gamma \kappa _2}\left(\frac{1}{r}\frac{1}{R_2}\right)\right]^{\frac{1}{\gamma 1}},$$
(15)
where $`1/R_2`$ is essentially the integration constant. The expressions for pressure and temperature can easily be written down. We now have to satisfy the conditions (11), (12) and (13). If $`B`$ is assumed to be given, then $`R`$, $`R_2`$ and $`\kappa _2`$ have to be found such these three conditions are satisfied. This is a mathematically well-posed problem, since we have to find three quantities from three equations. It is straightforward to solve this problem numerically. It may be noted that we can use the incomplete $`\beta `$ function defined as follows

$$\beta (x,\frac{3\gamma 4}{\gamma 1},\frac{\gamma }{\gamma 1})=\frac{1}{\beta (\frac{3\gamma 4}{\gamma 1},\frac{\gamma }{\gamma 1})}_0^x\xi ^{\frac{2\gamma 3}{\gamma 1}}(1\xi )^{\frac{1}{\gamma 1}}𝑑\xi $$

to write down (13) in the following form

$`\kappa ^{\frac{1}{\gamma 1}}R_{}^{\frac{3\gamma 4}{\gamma 1}}\left\{1\beta ({\displaystyle \frac{r_o}{R_{}}},{\displaystyle \frac{3\gamma 4}{\gamma 1}},{\displaystyle \frac{\gamma }{\gamma 1}})\right\}`$
$`=\kappa _2^{\frac{1}{\gamma 1}}R_2^{\frac{3\gamma 4}{\gamma 1}}\left\{\beta ({\displaystyle \frac{r_m}{R_2}},{\displaystyle \frac{3\gamma 4}{\gamma 1}},{\displaystyle \frac{\gamma }{\gamma 1}})\beta ({\displaystyle \frac{r_o}{R_2}},{\displaystyle \frac{3\gamma 4}{\gamma 1}},{\displaystyle \frac{\gamma }{\gamma 1}})\right\}`$ (16)
$`+\kappa ^{\frac{1}{\gamma 1}}R^{\frac{3\gamma 4}{\gamma 1}}\left\{1\beta ({\displaystyle \frac{r_m}{R}},{\displaystyle \frac{3\gamma 4}{\gamma 1}},{\displaystyle \frac{\gamma }{\gamma 1}})\right\},`$
which facilitates the numerical computation. On solving the problem numerically, we found that the radius increase was even much smaller than what we found in our order-of-magnitude estimate, and what is more puzzling is that the radius variation is less for a stronger magnetic field—a clearly unphysical result! Figure 1a shows $`\mathrm{\Delta }R=RR_{}`$ for various values of $`B`$ (the magnetic flux has been kept the same, so a larger $`B`$ implies a thinner magnetic layer). This unphysical result obviously is due to the fact that the formulation used by us left us no freedom to ensure that the temperature of the magnetic layer remained close to the temperature of the surroundings. In fact, $`\kappa _2`$ was found to be smaller than $`\kappa `$, making the density in some regions of the magnetic layer more than the density in the surrounding region as shown in Figure 1b (with temperatures reduced inside the magnetic layer). It is clear that Figure 1a cannot be taken to present very realistic values of radius variation with magnetic activity.

We next try to construct a solution by assuming the plasma beta $`\beta =8\pi p_2/B^2`$ to be constant inside the magnetic shell. The hydrostatic equation (9) inside the magnetic shell now takes the form

$$\frac{p_2}{r}=\frac{GM_{}\rho _2}{(1+\beta ^1)r^2}$$
(17)
Consequently $`\rho _2(r)`$ has the expression

$$\rho _2(r)=\left[\frac{(\gamma 1)GM_{}}{\gamma \kappa _2(1+\beta ^1)}\left(\frac{1}{r}\frac{1}{R_2}\right)\right]^{\frac{1}{\gamma 1}}$$
(18)
When we demanded that the conditions (11), (12) and (13) are satisfied, we were able to find only one solution of the problem, which corresponds to $`\beta ^1=0`$ (i.e. no magnetic field) for which there is no change in radius. This is certainly a consistent solution, but not the solution we have been looking for.

Although we have not succeeded in constructing a satisfactory hydrostatic model, we decided to present this brief discussion as a warning to future researchers that the simplest and most obvious things do not work here. There is no doubt that (9), (11), (12) and (13) should hold. Perhaps (10) should be replaced by a more realistic thermodynamic equation and we should use the physics of magnetic field creation in more detail to decide how $`B(r)`$ may vary within the magnetic layer. For these more complicated cases, we shall not be able to write down an analytical solution of density as we do in (15) or (18), and consequently the solution of the problem will be much more complicated. We felt that it will be pre-mature to attempt such a solution now. Given the uncertainties in the observational data, the order-of-magnitude estimate should more than serve our purpose at present.

## 5 Conclusion

Within the currently understood framework of solar MHD, we have not been able to provide a theoretical explanation of the controversial claim made by Dziembowski et al. (2001) that the Sun shrinks with increasing magnetic activity at the rate of 1.5 km yr-1. Spruit (1982, 1994) estimated the radius change due to heat blocking by sunspots near the solar surface. Not only the radius variation was found to be much smaller, the radius is expected to increase with the solar maximum. We considered the possibility if the excess magnetic pressure of flux tubes near the base of the SCZ can give us a different result. Even though traditional dynamo models regard the solar maximum as the time when the magnetic flux in the solar interior peaks and would suggest a radius increase at the time of the solar maximum, we pointed out the model of Nandy and Choudhuri (2002) provides a different scenario. Although this model would predict a decrease of the solar radius with increasing magnetic activity (i.e. with increasing loss of magnetic flux from the storage region), various reasonable assumptions give a radius decrease rate about one order of magnitude smaller than what is claimed by Dziembowski et al. (2001). It is true that our estimate is based on a very rough calculation. However, we believe that it is still an over-estimate rather than an under-estimate. Only if the magnetic field at the base of the SCZ was concentrated to values as high as $`10^6`$ G, the radius shrinkage would have been as large as reported by Dziembowski et al. (2001). There is no theoretical reason to expect a magnetic field as strong as $`10^6`$ G at the base of the SCZ. In fact, it is not easy to generate even a $`10^5`$ G magnetic field there (Choudhuri, 2003). We conclude that either the claim of Dziembowski et al. (2001) is incorrect, or else we do not understand some basic physics of the solar cycle.

## Acknowledgements

We wish to thank S. C. Tripathy for valuable discussions.


format: text

text:
# Axisymmetric dynamical models for SAURON and OASIS observations of NGC~3377

## 1 Introduction

NGC 3377 is a prototypical ‘disky’ E5-6 galaxy with ‘boxy’ outer isophotes in the Leo I group (e.g. Peletier et al. 1990) at an assumed distance of $`D=9.9`$ Mpc. It has a power-law central luminosity profile (Faber et al. 1997; Rest et al. 2001), and its total absolute magnitude of $`19`$ ($`B`$) is intermediate between that of the classical ‘boxy’ giant ellipticals and the ‘disky’ lower-luminosity objects (e.g. Kormendy & Bender 1996). Previous dynamical models of this galaxy suggest the presence of a central massive black hole (BH) of $`10^8`$ $`M_{\mathrm{}}`$ (Kormendy et al. 1998, hereafter K+98; Gebhardt et al. 2003, hereafter G+03), while the $`M_{}\sigma `$ relation (Tremaine et al. 2002, and references therein) predict a BH of $`\mathrm{4\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$.

In this paper, we present a unique combined set of nested integral-field spectroscopic observations of NGC 3377, from SAURON/William Herschel Telescope (WHT) and OASIS/Canada-France-Hawaii Telescope (CFHT), on which we will base our dynamical modeling.

Bacon et al. (2001a) presented the large-scale two-dimensional stellar kinematics of NGC 3377, obtained with the panoramic integral-field spectrograph SAURON as part of a representative survey of nearby early-type galaxies (de Zeeuw et al. 2002). The resulting kinematic maps cover $`32\mathrm{}\times 43\mathrm{}`$, with an effective spatial resolution of $`2\mathrm{}`$ FWHM, and reveal modest but significant deviations from axisymmetry, not only in the stellar motions, but also in the morphology and kinematics of the emission-line gas. This is unexpected, as lower-luminosity steep-cusped systems such as NGC 3377 were assumed to be axisymmetric (e.g. Davies et al. 1983; Gebhardt et al. 1996; Valluri & Merritt 1998, but see Holley-Bockelmann et al. 2002).

Here we also present for the first time observations of the inner $`5\mathrm{}\times 6\mathrm{}`$ of NGC 3377 with a spatial resolution of $`0\stackrel{}{.}6`$ FWHM, obtained at the CFHT with OASIS in its adaptive-optics-assisted mode. The resulting nested set of high-quality integral-field maps make NGC 3377 a nearly ideal case for detailed dynamical modeling, aimed at determining the stellar mass-to-light ratio $`M/L`$, the internal orbital structure, and also the BH mass $`M_{}`$.

In this paper, we construct axisymmetric dynamical models which incorporate all the integral-field kinematic data. Our approach is similar to that followed by Verolme et al. (2002) for M32, but we use an independent version of the modeling software, include all the available integral-field datasets and study the full line-of-sight velocity distributions (LOSVDs). By construction, we have to ignore any signature of non-axisymmetry, but we will include them in triaxial models in a follow-up work, in the way done by Verolme et al. (2003) for similar nested observations of NGC~4365. The comparison of the axisymmetric and triaxial models will help clarify which (if any) of the parameters based on the axisymmetric modeling are robust. This is relevant because many of the nearly twenty galaxies for which BH masses have been determined by means of (edge-on) axisymmetric models (Tremaine et al. 2002; Gebhardt et al. 2003) show similar signs of non-axisymmetry.

This paper is organized as follows. Sect. 2 presents the photometric data and the mass model adopted for NGC 3377, and Sect. 3 describes the kinematic data. The dynamical models are detailed in Sect. 4, and discussed in Sect. 5. Sect. 6 presents the conclusions.

## 2 Photometry and mass model

In this Sect., we describe the construction of an accurate mass model for NGC 3377, using the Multi-Gaussian Expansion (MGE) method (Monnet et al. 1992; Emsellem et al. 1994a) applied to both wide-field and high spatial resolution photometry. This allows us to probe the central regions as well as to cover the full extent of the galaxy. The detailed procedure is described in Emsellem et al. (1999, see also ).

### 2.1 Photometric data

The large-scale image of NGC 3377 shown in Fig. 1 (left panel) was kindly provided by R. Michard. It was obtained in 1995 with the 1.2 m telescope of the Observatoire de Haute-Provence (OHP) in the $`I`$-band, with a spatial resolution of $`2\stackrel{}{.}1`$ FWHM sampled at 0$`\stackrel{}{.}`$84, and a field-of-view (FoV) of $`4\mathrm{}\times 7\mathrm{}`$. This image was reduced in the usual way (bias, dark, flat-field, cosmic rays and cosmetics), and the sky contribution was estimated from the outer part of the frame and then subtracted. The flux normalization in $`L_{\mathrm{}}`$/pc2 was carried out using published aperture photometry (Poulain 1988; Goudfrooij et al. 1994), accessible via the HYPERCAT service111http://www-obs.univ-lyon1.fr/hypercat/.

For the high spatial resolution photometric data (see Fig. 1, right panel), we use $`I`$-band (F814W) HST/WFPC2 images, retrieved through the ST/ECF archives (PI Faber, ID 5512). The 5 individual exposures ($`2\times 80`$ s and $`3\times 350`$ s) were reduced in the standard way, and normalized in flux (in $`L_{\mathrm{}}`$/pc2) using the most recent PHOTFLAM conversion factors. The exposures were combined (with cosmics removal) after verifying they were properly centered.

The comparison between the two ground- and space-based flux calibrated exposures shows an offset of $`0.1`$ mag, mostly due to the difference in the zero points. We decided to use the HST/WFPC2 F814W image as a reference and renormalized the OHP image accordingly. The agreement between the two exposures is then excellent.

### 2.2 The MGE surface brightness model

We have used the Multi-Gaussian Expansion (MGE) method to model the surface brightness distribution for NGC 3377 using the photometric data described above. The routine provides an analytical model of the observed luminosity distribution, taking into account the convolution effect of the point spread function (PSF). The method assumes that both the PSF and the unconvolved surface brightness distribution of the galaxy can be described by a sum of two-dimensional Gaussians, whose best-fitting parameters are determined using an iterative approach.

Each two-dimensional Gaussian $`G_j^{}`$ is described by its maximum intensity $`I_j^{}`$, its dispersion $`\sigma _j^{}`$ along its major axis, its flattening $`q_j^{}`$, its center given by its coordinates $`X_j^{},Y_j^{}`$, and its position angle PAj. In the general case, all components are free to have different PA and centers, an $`N`$-Gaussian model thus depending on $`6N`$ free parameters.

The MGE method has several advantages:

* the proper account of the individual PSFs allows the simultaneous use of complementary datasets with different spatial resolution (space- and ground-based observations);
* its flexibility makes it well suited for complex, multi-component galaxies such as NGC 3377, which exhibits both diskyness and boxyness (e.g. Peletier et al. 1990; Scorza & Bender 1995);
* the fitting procedure allows linear and non-linear constraints on the parameters. In the present case, since we require the deprojected model to be axisymmetric, we force the two-dimensional model to be bi-symmetric by imposing components to share the same PA and center, resulting in 3 free parameters per Gaussian.

The PSF of the OHP exposure was approximated with an MGE fit (two concentric circular Gaussians) of the extracted image of a star located several arc-minutes from the galaxy, where its background is negligible. Using this model PSF, the large-scale surface brightness distribution of NGC 3377 is well described by a sum of 6 Gaussians. Since the nuclear regions are strongly affected by seeing ($`2\stackrel{}{.}1`$ FWHM), we only keep the three outer components of this model (with $`14\stackrel{}{.}5\sigma _j^{}72\stackrel{}{.}4`$, Gaussians $`G_{11}^{}`$ to $`G_{13}^{}`$ in the final model, see Table 1), and proceed by then fitting the central part using the HST/WFPC2 exposure.

The PSF of the HST/WFPC2-PC/F814W exposure was derived with the dedicated PSF simulator TinyTim222http://www.stsci.edu/software/tinytim/tinytim.html (v4.4), and consequently adjusted by a sum of three concentric circular Gaussians. The presence of dust lanes in the very central parts of NGC 3377, already noted by K+98, makes the determination of the inner Gaussians more sensitive, and hence the convergence process is slower. These dust lanes correspond to a maximal absorption of $`10\%`$, and the areas affected were discarded from the MGE fit. After subtraction of the properly aligned three-Gaussian ‘ground-based model’, the residual image was fitted by a sum of 10 Gaussians, describing the inner parts of NGC 3377 ($`0\stackrel{}{.}04\sigma _j^{}6\stackrel{}{.}04`$, Gaussians $`G_1^{}`$ to $`G_{10}^{}`$ in the final model, see Table 1). Gaussians $`G_5^{}`$ and $`G_8^{}`$ are significantly flatter than the others, betraying the presence of a disk.

The final MGE model of NGC 3377 consists of the resulting set of 13 Gaussians (see Table 1). Fig. 1 displays the contour maps of the $`I`$-band images with superimposed contours of the appropriately convolved MGE model, as well as cuts along the major axis. The resulting fit to the observed photometry is excellent. There is no sign of departures from axisymmetry in the central regions, and the PA of the major axis is constant with radius (e.g. Peletier et al. 1990). The total luminosity of the MGE model is $`L_I=\mathrm{1.1\hspace{0.17em}10}^{10}`$ $`L_{\mathrm{}}`$, yielding $`M_I=21.0`$. Taking a total $`BI=1.85`$ (Idiart et al. 2002), we get $`M_B=19.15`$, in good agreement with the value of $`M_B=19.2`$ available in the LEDA database333http://leda.univ-lyon1.fr (Paturel et al. 1989).

### 2.3 Deprojection

The MGE modeling allows the spatial luminosity density to be analytically derived from the surface brightness model, assuming an inclination angle $`i`$, and that the luminosity density associated with each individual three-dimensional Gaussian is stratified on concentric ellipsoids. Each two-dimensional Gaussian $`j`$ then uniquely deprojects into a three-dimensional Gaussian, whose parameters $`I_j`$, $`\sigma _j`$ and $`q_j`$ can be derived from $`I_j^{}`$, $`\sigma _j^{}`$, $`q_j^{}`$ and $`i`$ following relations provided in Monnet et al. (1992). The associated gravitational potential $`\mathrm{\Phi }`$ can be derived easily assuming a given (constant) mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ (see Emsellem et al. 1994a, for details).

The inclination $`i`$ of the best-fitting MGE model (see Table 1) is (formally) constrained by that of the flattest two-dimensional Gaussian to be $`90\mathrm{°}i76.6\mathrm{°}=\mathrm{arccos}(\mathrm{min}\{q_j^{}\}=0.232)`$. This constraint should however be treated with some caution since it is model dependent. However, no good MGE fit could be found with the additional requirement $`\mathrm{min}\{q_j^{}\}0.35`$. This sets a robust minimum inclination of $`i_{\mathrm{min}}70\mathrm{°}`$. In the following, we will only consider the edge-on model, *i.e.,* $`i=90\mathrm{°}`$, since the inclination hypothesis is likely not to be the most restrictive one (see discussion in Sect. 5).

## 3 Kinematical data

In this Sect., we summarize the stellar kinematic measurements obtained with the integral-field spectrographs OASIS and SAURON.

### 3.1 OASIS observations

NGC 3377 was observed on April 1 and 2 1998 with OASIS mounted on the adaptive optics (AO) bonnette PUEO of the Canada-France-Hawaii Telescope (see Bacon et al. 1998, 2001b, for a full description of OASIS). In order to take advantage of the better performance of adaptive optics in the red, we used the MR3 configuration covering the $`{}_{}{}^{}\mathrm{Ca}_{}^{}`$ triplet ($`8500`$ Å) region. We selected a spatial sampling of 0$`\stackrel{}{.}`$16 per hexagonal-shape lens, which provides a FoV of $`6\stackrel{}{.}2\times 5\stackrel{}{.}0`$. The instrumental setup is given in Table 2.

Seven 30 mn exposures centered on the nucleus (but slightly dithered to avoid systematics) were acquired, the AO being locked on the central cusp of the galaxy. The atmospheric conditions were photometric, and observations were carried out at low airmass ($`<1.16`$). Natural seeing conditions were mediocre with FWHM between 1″ and 1$`\stackrel{}{.}`$5, providing an AO corrected PSF with FWHM $`0\stackrel{}{.}6`$ (see Sect. 3.3.2). Neon arc lamp exposures were obtained before and after each object integration. Other configurations exposures (bias, dome flat-field, micro-pupil) were usually acquired at the beginning or the end of the nights. Twilight sky exposures were obtained at dawn or sunrise.

### 3.2 SAURON observations

NGC 3377 was observed on February 17 1999 with SAURON during its first scientific run on the William Herschel Telescope (Bacon et al. 2001a). The instrumental setup is summarized in Table 3. Four centered but slightly dithered 30 mn exposures were acquired, providing a total FoV of $`32\mathrm{}\times 43\mathrm{}`$ with an original sampling of 0$`\stackrel{}{.}`$94. The spectral coverage of SAURON includes the $`{}_{}{}^{}\mathrm{Mg}_{}^{}`$ absorption triplet as well as $`{}_{}{}^{}\mathrm{Fe}_{}^{}`$ and $`{}_{}{}^{}\mathrm{Ca}_{}^{}`$ absorption lines, and the \[O iii\] and $`{}_{}{}^{}\mathrm{H}_{}^{}`$$`\beta `$ emission lines when present.

### 3.3 IFS data reduction

#### 3.3.1 Spectra extraction and calibration

The integral-field spectroscopic data were reduced according to the usual procedure described in Bacon et al. (2001b) for OASIS and Bacon et al. (2001a) for SAURON. The standard reduction procedures include CCD preprocessing, optimal extraction of the spectra, wavelength calibration, spectro-spatial flat-fielding and cosmic-ray removal.

Given the small FoV of OASIS and the high surface brightness of the nucleus of NGC 3377 ($`\mu _V=17.7`$ at 3$`\stackrel{}{.}`$3), no sky subtraction was required for this instrument. Furthermore, no flux calibration was performed for OASIS, since this is not required for measurement of the stellar kinematics.

On the SAURON exposures, after spectral resolution rectification to a common value of 1.91 Å $`108`$ $`\mathrm{km}\mathrm{s}^1`$(see Bacon et al. 2001a, for details), the night-sky spectrum was estimated from the dedicated sky lenslets and subtracted from the object spectra.

#### 3.3.2 Spatial PSF estimates and datacubes merging

The knowledge of the spatial PSF is of particular importance in the modeling process, as we aim at using datasets with very different spatial resolutions. The spatial PSFs presented below have been estimated from the comparison between the images reconstructed from our IFS data and higher resolution reference frames (see details of the method in Bacon et al. 2001b). Once the individual frames were properly aligned and renormalized, the merging of multiple exposures was carried out according to the prescriptions detailed in Bacon et al. (2001a). Note that the data reduction software simultaneously provides the variance along each spectrum, which will used to estimate the local signal-to-noise in the datacube.

##### OASIS.

The PSF of the individual OASIS exposures is well approximated by a sum of two concentric circular Gaussians, and we estimated it by comparison with the HST/WFPC2 F814W image described in Sect. 2.1. We ranked the 7 OASIS PSFs, labeled from i to vii, according to their FWHM (see Fig. 2): 0$`\stackrel{}{.}`$56 for exposure v, $`0\stackrel{}{.}62`$ for exposures i,iii,iv and vi, and $`0\stackrel{}{.}73`$ for exposures ii and vii.

Single exposures have a signal-to-noise ratio $`S/N`$ which is insufficient for individual use in deriving the stellar kinematics. We have therefore constructed two merged datacubes using two different sets of exposures (trading-off between $`S/N`$ and resolution):

includes all the individual exposures except cubes ii and vii, in order to optimize the spatial resolution;

includes all seven individual exposures, in order to maximize $`S/N`$.

We approximate the PSF of each merged datacube by a sum of three Gaussians (see Table 4 and Fig. 3). It appears that there is not much difference in terms of effective spatial resolution between cubes ‘A’ and ‘B’, while the global $`S/N`$ of datacube ‘B’ is slightly higher. Hence we consider only cube ‘B’, which consists of 637 spectra with a final sampling of 0$`\stackrel{}{.}`$25.

##### SAURON.

All four individual SAURON exposures have similar spatial resolution, with a FWHM ranging from 2$`\stackrel{}{.}`$0 to 2$`\stackrel{}{.}`$5 (estimated from direct images obtained before and after the SAURON observations). All of them were therefore merged in a cube of 2957 spectra over-sampled at 0$`\stackrel{}{.}`$68 to allow a proper analysis of the spatial PSF.

We approximate the spatial PSF of the SAURON merged datacube by a sum of two Gaussians and estimated their parameters by comparison with the HST/WFPC2 F555W images obtained and reduced as described in Sect. 2.1 (see Table 4 and Fig. 3). The SAURON PSF is significantly non-gaussian: the two-Gaussian approximation has a global FWHM of 2$`\stackrel{}{.}`$15, while the best fit with a single Gaussian has a FWHM of 2$`\stackrel{}{.}`$62.

In order to reduce significantly the amount of data to handle during the modeling, while retaining all the initial information, the four individual SAURON data-cubes were then merged again in a final cube of 1534 spectra sampled at 0$`\stackrel{}{.}`$94, corresponding to a strict Nyquist sampling of the spatial PSF.

#### 3.3.3 Spatial binning

In order to increase the $`S/N`$, and reduce the number of independent apertures for the dynamical modeling, we applied an adaptive quadtree spatial binning to the SAURON datacube (see Appendix A). We did not apply this technique to the OASIS datacube, as we wished to retain the highest spatial resolution available. The SAURON binned datacube finally includes 475 independent spectra (see Fig. 14).

### 3.4 Stellar kinematics

We obtained reference stellar templates from dedicated exposures of HD 132737 for OASIS, and HD 85990 for SAURON, both of spectral type K0III. The stellar spectra were optimally summed over a spatial aperture of $`R=2\mathrm{}`$ to maximize the $`S/N`$ ($`1300`$ for OASIS and $`500`$ for SAURON). The resulting individual stellar spectra, as well as the final OASIS and SAURON NGC 3377 datacubes, were then continuum divided and rebinned in $`\mathrm{ln}\lambda `$. We used our own version of the Fourier Correlation Quotient (FCQ) method (Bender 1990; Bender et al. 1994) to derive the non-parametric LOSVD at every point of the final datacubes. When needed, the LOSVDs were parametrized using a simple Gaussian and complemented using 3rd and 4th order Gauss-Hermite moments $`h_3`$ and $`h_4`$ (van der Marel & Franx 1993; Gerhard 1993). We checked that similar results were obtained with an updated version of the cross-correlation method (Statler 1995), and that these did not significantly depend on the details of the continuum subtraction procedure.

Fig. 4 and 5 present the kinematic maps (mean velocity $`V`$, velocity dispersion $`\sigma `$, Gauss-Hermite moments $`h_3`$ and $`h_4`$) for OASIS and SAURON respectively, in the setup used for the dynamical modeling (see Sect. 4).

Estimates of the errors on each velocity bin of the non-parametric LOSVDs are obtained by means of a Monte-Carlo approach (see Copin 2000): for each galaxy position, a noise-free galaxy spectrum is built by convolving the template spectrum with the corresponding parametrized LOSVD. We then add noise realizations (consistent with the derived noise spectra of the datacube), and extract the LOSVD via FCQ. The error is defined as the variance of the distribution at each velocity bin after 100 realizations. If needed, errors on the measured kinematical parameters are derived from the same procedure, which produces reasonably realistic error-bars (see e.g. Fig. 6 and 7; see also de Bruyne et al. 2003).

#### 3.4.1 Comparisons with different datasets

We checked that the OASIS and SAURON datasets were consistent with each other, not only regarding their surface brightness distribution from which the PSFs were estimated, but also by comparing $`V`$ and $`\sigma `$ maps. The OASIS FoV is too small to allow proper convolution to the SAURON resolution. We therefore built a simple analytical model which approximates the ‘intrinsic’ $`V`$ and $`\sigma `$ fields, and used the MGE photometry model to convolve this model with the adequate PSFs. The agreement between the kinematics obtained with OASIS and SAURON is excellent, as shown in Fig. 6.

We also compared the OASIS and SAURON stellar kinematics with measurements based on long-slit data published by Bender et al. (1994, Calar Alto, slit width of 2$`\stackrel{}{.}`$1, spectral resolution of 46 $`\mathrm{km}\mathrm{s}^1`$, hereafter BSG94), K+98 (spectrograph SIS/CFHT, slit width $`0\stackrel{}{.}3`$–0$`\stackrel{}{.}`$5, spectral resolution $`40`$–60 $`\mathrm{km}\mathrm{s}^1`$) and Simien & Prugniel (2002, spectrograph CARELEC/OHP, slit width of 1$`\stackrel{}{.}`$5, spectral resolution of 25 $`\mathrm{km}\mathrm{s}^1`$, hereafter SP02). A proper comparison between these different datasets is difficult because of the very different spatial resolutions and instrumental setups. We therefore restricted these comparisons to the datasets which share *similar* spatial resolutions, *i.e.,* OASIS with K+98 and SAURON with BSG94 and SP02. The IFS data were thus binned according to the characteristics of the long-slit data and compared as shown in Fig. 7: the datasets are in good agreement, particularly considering the difficulty mentioned above. The (marginal) 0$`\stackrel{}{.}`$15-offset reported by K+98 in the mean velocity profile is not seen in our OASIS data.

We chose not to use the FOS observations of the nucleus of NGC 3377 (G+03, square aperture of 0$`\stackrel{}{.}`$21, spectral resolution of $`100`$ $`\mathrm{km}\mathrm{s}^1`$, no aperture illumination corrections applied) for three reasons: (a) the lack of extensive spatial information makes the comparison with such high-spatial resolution data hazardous, (b) the proper determination of aperture positioning is known to be a critical procedure (van der Marel et al. 1997), potentially leading to critical consequences for the dynamical modeling, (c) by excluding these measurements, our dynamical models are *fully* independant to the ones developed by G+03.

#### 3.4.2 First results from the SAURON and OASIS maps

As mentioned earlier, the OASIS and SAURON kinematics presented in Fig. 4 and 5 respectively are fully coherent and in excellent agreement with the long-slit datasets previously published, up to the highest spatial resolutions. The steep velocity gradient, unresolved at the SAURON resolution, reaches a maximum of $`110`$ $`\mathrm{km}\mathrm{s}^1`$ at $`r1\stackrel{}{.}5`$, and then decreases to a plateau of $`90`$ $`\mathrm{km}\mathrm{s}^1`$. The velocity dispersion peak is barely spatially resolved in the OASIS observations, and reaches $`170`$ $`\mathrm{km}\mathrm{s}^1`$. The $`h_3`$ coefficient is anti-correlated with the mean velocity as expected from the superposition of a bright dynamically cold component on a hotter spheroid.

Integral-field spectroscopy enables us to analyse the ‘morphology’ of the kinematics, something not possible using long-slit spectra. Both OASIS and SAURON mean velocity fields reveal a twist of the zero velocity curve of $`10^{}`$ with respect to the photometric minor-axis up to $`r10\mathrm{}`$, while the photometry does not show any variation of position angle (Peletier et al. 1990). As a consequence, the velocities along the minor-axis are not null, a fact not reported so far444Long-slit (unpublished) data showed a weak rotation along the minor-axis of NGC 3377 ., and which is clearly visible with an amplitude of $`7`$ $`\mathrm{km}\mathrm{s}^1`$ for OASIS and $`12`$ $`\mathrm{km}\mathrm{s}^1`$ for SAURON.

The SAURON velocity dispersion-map also presents interesting morphological features: while the $`\sigma `$-peak is elongated along the photometric minor-axis – consistent with the diskyness of the light distribution –, it displays a twist of about $`10^{}`$ in the direction *opposite* to the kinematic misalignment. These modest but significant departures from axisymmetry are probable signatures of a triaxial intrinsic shape. This is also apparent in the ionized gas distribution and motions, which exhibit a noticeable spiral-like morphology and strong departures from circular motions (Bacon et al. 2001a). These issues will be discussed in more detail in a forthcoming paper .

## 4 Dynamical models

### 4.1 Schwarzschild models

We construct dynamical models based on the orbit superposition technique of Schwarzschild (1979, 1982). Similar orbit-based models were constructed to measure the masses of central BHs (van der Marel et al. 1998; Cretton & van den Bosch 1999; Verolme et al. 2002; Cappellari et al. 2002; Gebhardt et al. 2003), and dark halo parameters (Rix et al. 1997; Gerhard et al. 1998; Saglia et al. 2000; Kronawitter et al. 2000). Our implementation is described in detail in Cretton et al. (1999, hereafter C99). We summarize it briefly here.

We sample the stellar orbits using a grid in integral space, which includes: the energy $`E`$, the vertical component of the angular momentum $`L_z`$ and an effective third integral $`I_3`$. We use 20 values of $`E=1/2R_c\mathrm{\Phi }/R+\mathrm{\Phi }(R_c,0)`$, sampled through the radius of the circular orbit $`R_c`$. We take a logarithmic sampling of $`R_c`$ in $`[0\stackrel{}{.}01,300\stackrel{}{.}0]`$, since more than 99% of the total mass of our MGE model lies inside this range. 14 values of $`L_z`$ per $`E`$ in $`[L_{z,\mathrm{max}},+L_{z,\mathrm{max}}]`$ and 7 values of $`I_3`$ per ($`E,L_z`$) were adopted (see C99).

The orbit library is constructed by numerical integration of each trajectory for an adopted amount of time (200 periods of the circular orbit at that $`E`$) using a Runge-Kutta scheme. During integration, we store the fractional time spent by each orbit in a Cartesian ‘data-cube’ $`(x^{},y^{},v_{\mathrm{los}})`$, where $`(x^{},y^{})`$ are the projected coordinates on the sky and $`v_{\mathrm{los}}`$, the line-of-sight velocity. We use an $`E`$-dithering scheme (see C99) to make each orbit smoother in phase space. The data-cube of each individual orbit is further convolved with the PSF and eventually yields the *orbital* LOSVD $`_{(x^{},y^{})}(v_{\mathrm{los}})`$ at each position $`(x^{},y^{})`$ on the sky.

In C99, we adopted a parametrized form for the (orbital and observed) LOSVDs using the Gauss-Hermite series expansion. However, as emphasized in Cretton & van den Bosch (1999), some problems may arise in the modeling of dynamically ‘cold’ systems (*i.e.,* with a high $`V/\sigma `$), where the use of the linear expansion can lead to spurious counter-rotation. To avoid this potential shortcoming, we choose to constrain the dynamical models using directly the full non-parametrized LOSVDs. As a consequence, the maps presented in Fig. 4 and 5 are shown for illustration purposes, but were not used in the modeling technique.

Orbital occupation times are also stored in logarithmic polar grids in the meridional plane and in the $`(x^{},y^{})`$ plane, to make sure the final orbit model reproduces the MGE mass model. The mass on each orbit is computed with the NNLS algorithm (Lawson & Hanson 1974), such that the non-negative superposition of all orbital LOSVDs aims at reproducing the observed LOSVDs within the errors. In addition, the model has to fit the intrinsic and projected MGE mass profiles. Smoothness in integral space can be enforced through a regularization technique (see C99). Models with different values of BH mass $`M_{}`$ and mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ are constructed and compared to the data. As mentioned in Sect. 2.3, the mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ is assumed to be constant over the whole extent of the galaxy. The quality of the fit is assessed through a $`\chi ^2`$-scheme and we use the $`\delta \chi ^2=\chi ^2\chi _{\mathrm{min}}^2`$ statistic to assign confidence values to the iso-$`\chi ^2`$ contours (e.g. Cretton et al. 2000; Gebhardt et al. 1996).

The main differences with respect to models previously published are (a) the use of the full non-parametrized LOSVDs (but see Gebhardt et al. 2000; Bower et al. 2001; Gebhardt et al. 2003), (b) the two-dimensional coverage of the data (but see Cappellari et al. 2002; Verolme et al. 2002). Most previous studies made use of long-slit 1D-data at a few position angles (e.g. van der Marel et al. 1998), but the advance of IFU spectrographs (e.g. OASIS/CFHT, SAURON/WHT, FLAMES and VIMOS at the Very Large Telescope, etc.) will deliver two-dimensional data for many objects. A more detailed study of the effect of using such two-dimensional constraints will be presented in a companion paper (Cretton & Emsellem, in preparation).

### 4.2 SAURON constraints with original errors

We compute a grid of orbit libraries with $`M_{}`$ ranging from $`0.0`$ to $`\mathrm{1.82\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$ and $`\mathrm{{\rm Y}}_I`$ between $`2.2`$ and $`2.6`$555A Jeans model based on the previous MGE model (Sect. 2.2) gives $`\mathrm{{\rm Y}}_I2.5`$ (see e.g. Emsellem et al. 1999).. First, we constrain these models only with the SAURON dataset and $`\delta \chi ^2`$-contours are showed in Fig. 8 (left panel). At 99.7% confidence level, the allowed $`M_{}`$ is *smaller* than $`\mathrm{2\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$. This appears to be surprisingly small, given the low spatial resolution ($`2\stackrel{}{.}1`$ FWHM) of the SAURON dataset. Indeed, in the case of NGC 3377 with a characteristic central velocity dispersion of $`\sigma 135`$ km/s and an assumed distance of $`D=9.9`$ Mpc, this spatial resolution corresponds roughly to the radius of influence of a $`\mathrm{1.6\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$ BH (de Zeeuw 2001).

In order to check that the tightness of the $`\delta \chi ^2`$-contours is not an artifact of our choice of technical implementation, we have constructed several orbit libraries, changing one parameter of the library at a time: (a) the length of each orbit integration has been increased up to 1000 radial periods, (b) the $`E`$-sampling of the library has been refined, (c) the individual orbital time-steps have been decreased, and (d) the size of the orbit library has been increased by eight-fold, $`(N_E,N_{L_z},N_{I_3})=(40,28,14)`$ instead of $`(20,14,7)`$. For the latter model, the resulting $`\delta \chi ^2`$-contours are shown in the right panel of Fig. 8: although the orbit library was expanded by a factor of 8, the upper limit for $`M_{}`$ is still surprisingly small. In fact, none of the above modifications affected the tightness of the $`\delta \chi ^2`$-contours in a significant way. We conclude that this peculiar result is not related to specific technical details of our modeling method.

### 4.3 Departures from axisymmetry

To test further the origin of the tightness of the $`\delta \chi ^2`$-contours for the SAURON dataset, we constructed fake constraints drawn from an isotropic distribution function $`f(E,L_z)`$, using the Hunter & Qian (1993) method (hereafter HQ). The even part of such a DF is uniquely specified by the mass profile and we choose the odd part such as to mimic the true observed SAURON kinematics. A $`10^8`$ $`M_{\mathrm{}}`$ BH has been included into the HQ model. HQ-LOSVDs are then computed in the SAURON setup: they are convolved with the SAURON PSF and binned into the SAURON spatial elements. This artificial dataset looks very similar to the SAURON data, but a dynamical fit gives very different $`\delta \chi ^2`$-contours, allowing BH masses up to $`\mathrm{1.5\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$ at 99.7% confidence level. For this fit, we have used ad-hoc error-bars equal to 5% of the largest LOSVD value in each pixel. In each SAURON pixel, these errors are therefore independent of $`v_{\mathrm{los}}`$. With such a choice of errors, we are conservative in the sense that real observed errors are (on average) three times higher and would therefore induce even wider $`\delta \chi ^2`$-contours.

By construction, the HQ data correspond to a perfectly axisymmetric galaxy, and are therefore fully \[anti-\]bi-symmetric, *i.e.,* $`_{(x^{},y^{})}(v_{\mathrm{los}})=_{(x^{},y^{})}(v_{\mathrm{los}})=_{(x^{},y^{})}(v_{\mathrm{los}})`$. In that sense, as mentioned earlier, the observed SAURON data show noticeable departure from axisymmetry, and the fit of an axisymmetric model to them will significantly increase both $`\chi ^2`$ and $`\delta \chi ^2`$.

If one still wishes to model non-axisymmetric data with an axisymmetric code, there are two possibilities: (a) symmetrize the data, while keeping the statistical error-bars, (b) increase the error-bars to encompass the systematic errors (if any) between the 4 quadrants. While the first solution was adopted by G+03, we chose the second option to keep as much as possible the original spatial extent of the integral-field data. For the central parts of the SAURON dataset (roughly $`|x^{}|12\mathrm{}`$ and $`|y^{}|6\mathrm{}`$), for which measures from the 4 quadrants are accessible, the error-bar including systematics is taken as $`\sigma _{+\mathrm{syst}}^2=\sigma ^2_4+\mathrm{Var}_4(\sigma ^2)`$, where $`\mathrm{Var}_4`$ is the variance between the four LOSVDs of the four quadrants and $`\sigma ^2_4`$ is the mean statistical error-bar of the four LOSVDs. In the outer parts, the spatial binning and the incomplete coverage of the SAURON data make the computation more difficult: we estimate the systematic errors by taking the immediate neighbor spatial element (after folding all 4 quadrants on one) if there is no symmetric correspondent in the other quadrants.

As expected, we obtain the same contours using either original or symmetrized SAURON data when the error-bars including systematics are used.

### 4.4 IFU constraints with errors including systematics

##### SAURON dataset.

Fig. 9 (left panel) shows the $`\delta \chi ^2`$-contours using only the SAURON data with the errors including systematics. While the constraint on the upper limit of $`M_{}`$ is noticeably relaxed (the largest allowed BH mass is now $`10^8`$ $`M_{\mathrm{}}`$), $`\mathrm{{\rm Y}}_I`$ is strongly restricted to the range $`[2.05,2.45]`$ (99.7% confidence level), due to the large SAURON FoV.

##### OASIS dataset.

The OASIS data also show signs of non-axisymmetry (see Fig. 4), so we applied the procedure described in the previous Sect. to compute systematic error-bars. Fig. 9 (middle panel) shows the $`\delta \chi ^2`$-contours constrained only with the OASIS data (restricted to the central $`4\mathrm{}\times 2\mathrm{}`$ part, see Fig. 10). The 99.7% confidence level is much larger than in the SAURON case: as could be expected (e.g. G+03) $`\mathrm{{\rm Y}}_I`$ is very weakly constrained due to the small extension of the data. However, a model *without* BH is excluded by the OASIS data at a 99.7% confidence level.

##### Combined datasets.

Fig. 9 (right panel) shows the results of a combined fit on both datasets. Since this fit includes all available IFU data, it presumably provides the best estimate of the BH mass and the stellar $`\mathrm{{\rm Y}}_I`$. To improve such estimates based on a grid interpolation, we compute additional models corresponding to new values of $`\mathrm{{\rm Y}}_IM_{}`$. The best fitting model has $`M_{}=\mathrm{7.0\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ and $`\mathrm{{\rm Y}}_I=2.05`$ (see the best-fit maps in Fig. 10 and 11, to be compared with the observed maps in Fig. 4 and 5). The minimum and maximum BH mass allowed by these data (99.7% confidence level) are $`\mathrm{2.5\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ and $`\mathrm{1.25\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$ respectively.

### 4.5 Regularized models

In previous Sects., we have explored models with no requirements on the smoothness of the solution. As a result, the orbital weights are jagged and can strongly vary from one orbit to the next. One can include regularization constraints in order to obtain smoothly varying solutions (see e.g. C99). From such regularized solutions smooth functions characterizing the internal dynamical structure can be computed (see Fig. 13). But there is a much more important reason to derive regularized models. Recently, Valluri et al. (2002) have shown that (unregularized) orbit-based models do not provide reliable BH mass estimates, because the $`\chi ^2`$-contours on which they are based depend strongly on the number of orbits in the library: if one increases the number of orbits, the contours get wider and the indeterminacy on the BH grows. In a forthcoming paper, we have explored this issue with component-based $`f(E,L_z)`$ models and have also observed such a significant widening of the contours when the orbit library is expanded. Our tests also show that regularization may provide a way to ‘stabilize’ the $`\chi ^2`$-contours, while increasing the number of orbits . Although we do not have yet a clear physical justification for the regularization scheme, we decided to apply such a procedure to construct regularized models of NGC 3377. We refer the reader to the companion paper for a discussion of these issues, outside the scope of this paper.

Fig. 12 shows the combined fit, including regularization constraints. The level of regularization has been adjusted from a comparison with the HQ models according to the prescription of C99 (their Fig. 8). The contours are then very similar to the case without regularization, and the best fitting values are almost unchanged with respect to the unregularized model with the exception of the slightly restricted maximum allowed BH mass, which is now $`\mathrm{1.1\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$. The values we adopt are thus $`M_{}`$ in the range 2–$`\mathrm{11\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ (99.7% confidence level), with a best fit value of $`\mathrm{7\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ for a mass-to-light ratio of $`\mathrm{{\rm Y}}_I=2.1\pm 0.2`$, and an assumed distance of 9.9 Mpc. Note that at the 68% confidence level, our model is well constrained with $`M_{}`$ between $`6`$ and $`\mathrm{8\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$.

### 4.6 Internal dynamical structure

To describe the internal dynamical structure of NGC 3377, we compute the velocity dispersion profiles of the best-fit orbit model: $`\sigma _r`$, $`\sigma _\theta `$ and $`\sigma _\varphi `$, where $`(r,\theta ,\varphi )`$ are the usual spherical coordinates (see Fig. 13). We use the regularized models constrained on both SAURON and OASIS datasets with errors including systematics (see previous Sect.). In Fig. 13, we also plot the ratios $`\sigma _r/\sigma _\theta `$, and $`\sigma _r/\sigma _\varphi `$ to better estimate the anisotropy.

The meridional plane $`(R,z)`$ is divided into a polar grid with seven angular sectors for $`\theta [0,\pi /2]`$ (see Sect. 4.1); the first sector is close to the symmetry axis, while the last one is near the equatorial plane. We decide to discard the sector closest to the symmetry axis (*i.e.,* the first one), because only few orbits with very small $`L_z`$ can reach it. To reduce the noise, we average sectors two by two: #2 with #3, #4 with #5 and finally sectors #6 with #7, in the left, middle and right panels respectively (see Fig. 13).

In each panel, the vertical dashed lines indicate the extent of the kinematic data: it thus delineates the region where the models are directly constrained. Going from the symmetry axis toward the equatorial plane, the best-fit model (with $`M_{}=\mathrm{7\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$) shows an increasing radial anisotropy. It is however relatively close to isotropy, with dispersion ratios $`\sigma _r/\sigma _\theta `$ and $`\sigma _r/\sigma _\varphi `$ of about 1.2. Except for the peak of $`\sigma _r/\sigma _\varphi `$ in the central part of the equatorial plane profile, the model shows a rather constant anistropy along each (averaged) sector.

Although our mass model contains two flat Gaussian components ($`G_5^{}`$ and $`G_8^{}`$, see Table 1) with $`q_j^{}<0.24`$, the internal structure of our best-fit model of NGC 3377 does not exhibit a strong dynamical signature of a disk like in NGC~4342 (see Fig. 13 and §7.3 of Cretton & van den Bosch 1999). Indeed, in the case of NGC 3377, these two flat components never dominate the local mass density, whereas in NGC 4342 the outer disk starts to dominate outside of 4″, inducing a radial anisotropy.

G+03 observed a ‘peak’ in $`\sigma _r/\sigma _\theta `$ between $`r=0\stackrel{}{.}1`$ and 1″ in the equatorial plane (see their Fig. 13). This ratio reaches 1.5 in their best-fit model, whereas we observed values below 1.2 in the same radial range. The differences are certainly due to the significantly different spatial coverage of the datasets used to constrain the models.

### 4.7 Comparison with previous models

The dynamics of NGC 3377, and the presence of a central massive BH, have been first studied by K+98 with isotropic models, and in greater details by G+03. The modeling technique (3$`I`$-Schwarzschild) and initial assumptions (axisymmetry, edge-on, constant mass-to-light ratio) used by G+03 are similar to the ones we applied on our NGC 3377 integral field data. A few important differences should be emphasized though. Firstly, G+03 constrained the spatial luminosity to be constant on homothetic ellipsoids, which may cause significant differences in the resulting dynamical model (and internal structure). However, considering that isophotes of NGC 3377 are reasonably well approximated by ellipses (with constant ellipticity), this is probably not critical in the present case. Secondly, G+03 included the FOS apertures in their list of observable constraints. As noted by G+03, the inclusion of HST measurements improves the significance of the BH detection, but does not alter significantly the upper limit of the mass estimates. The uncertainty on the actual location of the FOS apertures and its corresponding kinematics (see Sect. 3.4.1) furthermore motivated our choice of excluding this dataset.

G+03 found a BH mass of $`M_{}=1.0_{0.1}^{+0.9}\mathrm{\hspace{0.17em}10}^8`$ $`M_{\mathrm{}}`$ ($`1\sigma `$ range, one degree of freedom, for an assumed distance of 11.2 Mpc). This translates into a mass range for the $`M_{}`$ of about $`8\mathrm{17\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ at a distance of 9.9 Mpc. However, as explained by the authors, the lower limit on $`M_{}`$ is solely constrained by the 2 nuclear FOS LOSVDs: the fit on their ground-based data – stellar $`V`$ and $`\sigma `$ measurements from K+98 (see Sect. 3.4.1) – is consistent with the absence of central BH (their Fig. 7). Our best fit value of $`M_{}=\mathrm{7\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ is just at the edge of the 95% confidence contours derived by G+03, but outside their 68% confidence mass range (see their Fig. 3).

Regarding the mass-to-light ratio, the best-fit value of $`\mathrm{{\rm Y}}_V=2.9_{0.6}^{+0.1}`$ at $`D=11.2`$ Mpc found by G+03 (99.7% level) corresponds to $`\mathrm{{\rm Y}}_I\mathbf{2.16}_{0.45}^{+0.08}`$ at 9.9 Mpc (for a mean color index $`(VI)=1.14`$, Goudfrooij et al. 1994666Note that although the colour profiles presented in Goudfrooij et al. (1994) are right, the corresponding $`\overline{VI}`$ value quoted in Table 5 is incorrect (Goudfrooij, private communication)., and a standard Cousins $`(VI)_{\mathrm{}}=0.685`$, Bessell 1979), consistent with our best fit value of $`\mathrm{{\rm Y}}_I=2.1`$. At the 68% level, we obtain a narrow range of allowed $`\mathrm{{\rm Y}}_I`$, between 2.02 and 2.12, a significant improvement compared to the corresponding range of 1.71–2.24 found by G+03. The main reason for our better mass-to-light ratio constraint certainly lies in the spatial coverage of our kinematical datasets: G+03 only disposed of major-axis (long-slit) profiles (with the addition of the FOS apertures). As shown in Fig. 9, the SAURON dataset with its two-dimensional coverage on a relative large field strongly constrains $`\mathrm{{\rm Y}}`$ (assumed to be constant).

Using the combined SAURON and OASIS datasets, we therefore significantly improve the constraints on both $`M_{}`$ and $`\mathrm{{\rm Y}}`$. While our ground-based data spatial resolution is coarser than the HST spatial-grade resolution, we believe, given our aforementioned concerns regarding the uncertainties on critical kinematic measurements derived from FOS apertures, that our BH mass estimate is overall more robust than G+03’s. In addition, we think that our mass-to-light ratio estimate is better because of the extended spatial coverage of our kinematical datasets.

## 5 Discussion

The kinematic data used to constrain the dynamical modeling presented in this paper were exclusively obtained from two integral-field spectroscopic observations of NGC 3377. The full spatial coverage allows a proper estimate of the respective PSF (see Sect. 3.3.2), and a detailed comparison of the datasets, insuring their internal consistency (see Sect. 3.4.1). Furthermore, the spatial complementarity of the datasets – SAURON covers a relatively large FoV while OASIS probes the central parts with sharper spatial resolution – plays a key role in the modeling process (see Sect. 4.4).

Our implementation of the Schwarzschild technique has been extended to support 2D-coverage and non-parametrized LOSVD fit, and includes regularization constraints minimizing the ‘contours widening’ effect described by Valluri et al. (2002). However, our models still depend on three major assumptions, namely a given inclination, axisymmetry, and a constant mass-to-light ratio.

##### Inclination.

The galaxy NGC 3377 is known *not to be* strictly edge-on, as strongly hinted from the ionized gas distribution (Bacon et al. 2001a, if the gas is confined in the equatorial plane). However, NGC 3377 is certainly close to edge-on ($`i70^{}`$, see Sect. 2.3) given its apparent flattening (E5-6). Running a full grid of dynamical models over the allowed range of inclinations (as in, e.g. Gebhardt et al. 2000; Verolme et al. 2002), might not be relevant if, as suspected, the assumption of axisymmetry is the most restrictive one.

##### Axisymmetry.

As already mentioned, a $`10^{}`$ kinematic misalignment is observed both in the OASIS and SAURON datasets, in contradiction with an intrinsic axisymmetric morphology (while the photometry does not show any position angle twist). This non-axisymmetry could be due to the presence of an inner bar . While the observed non-axisymmetry weakens the results of our axisymmetric models, it is not clear whether it is critical for the global understanding of the intrinsic dynamics of this ‘nearly oblate’ galaxy (Davies et al. 1983) – contrary to more extreme cases (e.g. NGC 4365, Davies et al. 2001; Verolme et al. 2003) – and wether it prevents our $`(M_{},\mathrm{{\rm Y}})`$ best-fit parameters to be taken into consideration.

##### Mass-to-light ratio.

We can compare the dynamical estimate of our $`\mathrm{{\rm Y}}`$ (obtained via Schwarzschild modeling) to the one independently inferred from observed colours and line strengths. Idiart et al. (2003) estimated the metallicity of NGC 3377 using broad-band colours, and derived a central and global $`[Fe/H]`$ of $`+0.0`$ and $`0.2`$ respectively. Using Vazdekis et al. (1996) stellar population synthesis models with observed colors $`BV0.86`$ and $`VI1.14`$, this corresponds to a $`\mathrm{{\rm Y}}_V`$ in the (loose) range 2.7–3.7. Our best-fit model has $`\mathrm{{\rm Y}}_I=2.1`$, or $`\mathrm{{\rm Y}}_V=3.2`$, fully consistent with these estimates (but see discussion in Maraston 1998; Thomas et al. 2003). We also checked that the SAURON line strength values ($`{}_{}{}^{}\mathrm{H}_{}^{}`$$`\beta `$, $`{}_{}{}^{}\mathrm{Mg}_{}^{}`$b, $`{}_{}{}^{}\mathrm{Fe}_{}^{}`$5015, $`{}_{}{}^{}\mathrm{Fe}_{}^{}`$5270) observed for NGC 3377 are within the ranges predicted by the Vazdekis et al. models. This independent check is important as it links the stellar populations to the dynamical observables.

Regarding our central BH mass measurement, the only independent estimate we can provide, beside comparison to previously published dynamical models (see Sect. 4.7), is by using the $`M_{}\sigma `$ relation: Merritt & Ferrarese (2001) and Tremaine et al. (2002) used stellar dispersion (central or effective) values of 131 and 145 $`\mathrm{km}\mathrm{s}^1`$ respectively, providing $`M_{}=1.47(\pm 0.72)\mathrm{\hspace{0.17em}10}^7`$ and $`M_{}3.7(\pm 1.7)\mathrm{\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$, respectively. The range of allowed $`M_{}`$ derived from our Schwarzschild models overlap with both estimates, although our best-fit value is clearly outside these two ranges.

Studying the importance of the effect of triaxiality is outside the scope of this paper and will be done in the future using generalized modelling techniques . In that respect, it is noticeable that only a two-dimensional coverage can properly trace the non-axisymmetry. This has consequences not only for the specific case of NGC 3377, the non-axisymmetry of which was ignored before, but probably for many other galaxies and related studies.

## 6 Conclusions

We have presented a unique set of nested integral-field spectroscopic observations of the E5-6 ‘disky/boxy’ cuspy galaxy NGC 3377, obtained with the OASIS and SAURON instruments. Both sets display regular stellar kinematics, with indications for an inner disk, coherent with the diskyness of the light distribution up to $`30\mathrm{}`$. However, the IFS observations also reveal a significant twist of the line of zero velocity of $`10^{}`$ up to $`r10\mathrm{}`$ – with no corresponding isophotal twist – indicative of a moderately triaxial intrinsic shape or of a central bar.

Based on this complementary and intrinsically coherent dataset, we have constructed general (three-integral) axisymmetric dynamical models for this galaxy using the orbit superposition technique of Schwarzschild, non-parametrically fitting the full LOSVDs delivered by the 3D spectrographs. In order to properly take into account the observed non-axisymmetry of the kinematic maps, we have added a systematic component to the Monte-Carlo-computed statistical errors on the LOSVDs used in the fit. The best-fit model has a BH mass of $`M_{}=7_5^{+4}\mathrm{\hspace{0.17em}10}^7`$ $`M_{\mathrm{}}`$ and a mass-to-light ratio $`\mathrm{{\rm Y}}_I=2.1\pm 0.2`$ (99.7% confidence level) for an assumed distance of $`D=9.9`$ Mpc.

The use of fully coherent IFS measurements on the observation side, and of proper errors and regularization constraints on the modeling side allowed us to derive this significantly improved constraints on both $`\mathrm{{\rm Y}}`$ and $`M_{}`$. However, some assumptions used in the modeling – the most important presumably being the axisymmetry – are in noticeable contradiction with our IFS observations. Therefore, BH mass estimates for NGC 3377 based on axisymmetric models should be considered with caution.

###### Acknowledgements.

YC’s research was supported through a European Community Marie Curie Fellowship at Leiden Observatory. NC thanks H.-W. Rix for many insightful discussions. EE wishes to warmly thank R. Michard for making available the ground-based photometry of NGC 3377, and the visitor program at ESO, during which part of this work was done. We also would like to thank Michele Cappellari, Roger Davies, Harald Kuntschner and Tim de Zeeuw for a critical reading of the manuscript.

## Appendix A Adaptive quadtree two-dimensional binning

In order to increase the $`S/N`$ in the outer parts of the SAURON exposure, we used an adaptive spatial binning algorithm. A detailed account of binning techniques is described in Cappellari & Copin (2003).

While filtering (*i.e.,* ‘smoothing’) spectrographic data is equally easy for long-slit and integral-field spectrographs, data binning is much trickier for two-dimensional data, due to the topological constraints (*i.e.,* the proper tiling of the plane). Still, proper binning of two-dimensional data is necessary when we wish to compare them with theoretical or numerical models.

The first step is to estimate the flux and mean $`S/N`$ at each spatial element (hereafter *spaxel*). In the present case, each spaxel $`i`$ of the input cube is associated to the spectrum $`S_i(\lambda )`$, and one can compute the total flux $`𝒮_i=S_i(\lambda )d\lambda `$ by integrating along the spectral dimension, as well as the mean signal-to-noise ratio $`\eta _i=S_i(\lambda )/\sigma _i(\lambda )_\lambda `$, where $`\sigma _i(\lambda )`$ is the chromatic variance as estimated during the data-reduction process. Furthermore, the real $`(𝒮_i,\eta _i)`$-spaxels are embedded in a $`2^n\times 2^n`$-grid, otherwise completed with $`(0,0)`$-spaxels.

The algorithm is then as follows:

All the spaxels of the input grid are associated in a single bin, with flux $`_1=_i𝒮_i`$ and mean signal-to-noise $`\zeta _1=_1/_i(𝒮_i^2/\eta _i^2)`$. Hopefully, this 1-bin has a $`S/N`$ high enough to be divided further and initialize the process.

If $`m`$-bin $`j`$ is indeed a bin – *i.e.,* contains at least 4 spaxels – and has high enough a $`S/N`$ with respect to the target $`S/N`$, it can be divided further into 4 subsequent $`(m+1)`$-bins. Since the initial grid had a power-of-two side length, all the child-bins also have a power-of-two side length.

All the bins which are not completely covered by real data (in our case, 80%) or whose $`S/N`$ is lower than a given threshold are then discarded.

This method leads naturally to an *Unbalanced Quadtree spatial binning* of the original datacube (e.g. Bern & Eppstein 1995). Contrarily to the optimal binning presented in Cappellari & Copin (2003), this quadtree binning does not produce a very homogeneous $`S/N`$ distribution over the field. However, it provides simple shapes – squares – for the final bins, and therefore simplify the comparison with models.

For the specific case of the SAURON NGC 3377 data, we have used a target $`S/N`$ of 50 and a minimum $`S/N`$ of 30. This resulted into 475 independent binned spectra as shown in Fig. 14.