text:
# Gas Dynamics in the Milky Way: Second Pattern Speed and Large–Scale Morphology

## 1 Introduction

Observations of cold gas in the Milky Way \[MW\] have contributed substantially to our understanding of MW structure. No other tracer is observed in as large a part of the MW as are gas clouds. Longitude-velocity \[lv\] diagrams (Hartmann & Burton , Dame et al. ) show the distribution of gas velocities as a function of galactic longitude $`l`$, integrated over some range in latitude $`b`$. By observing the MW in different spectral lines, this gas can be traced at substantially different densities. The largest absolute velocity as a function of $`l`$ defines the terminal velocity curve (TVC). In an axisymmetric galaxy, the gas at these velocities is found at the “tangent point” where the line of sight is tangential to a circle around the Galactic centre. From this the rotation curve can be determined. However, due to the bar and spiral perturbations in the MW potential, the gas has substantial non-circular velocities, which are most evident in the central $`10^{}20^{}`$, due to the bar, but also as “bumps” in the TVC where spiral arm tangents perturb the gas flow by $`10`$-$`20\mathrm{km}\mathrm{s}^1`$. At sub-TVC velocities, crowding in both position and in velocity produces ridge-like structures in the ($`l,v`$) diagram. The Galactic spiral arms are visible as straight or curved such ridges.

A number of attempts have been made to model these observations. One group is formed by analytic models of spiral structure. The first exhaustive analytical formulation of a spiral arm theory was developed by Lin & Shu  and applied to the MW by Lin et al. . They proposed a two-armed model with pitch angle $`6^{}`$ and a pattern speed for the spiral structure $`\mathrm{\Omega }_{sp}13.5\mathrm{km}/\mathrm{s}/\mathrm{kpc}`$. Amaral & Lépine  fitted the rotation curve of the MW to an analytic mass model and found a self-consistent solution with a combined two- and four-armed spiral structure. In Lépine et al.  they extended this model to allow for a phase difference between the two-armed and the four-armed spiral pattern.

The second group, numerical simulations of the Galactic gas flow, also have a long tradition. A recent example is the smoothed particles hydrodynamics (SPH) models of Fux , who evolved a gas disk inside a self-consistent $`N`$-body model scaled to the COBE/DIRBE K-band map of the MW and the radial velocity dispersion of M giants in Baade’s window. The resulting gas flow was transient, but at specific times resembled closely a number of observed arms and clumps in the bar region. Weiner & Sellwood  compared predictions from fluid dynamical simulations for analytic mass densities with the observed outer velocity contours of the HI $`(l,v)`$ diagram to constrain the pattern speed and bar angle. Englmaier & Gerhard  \[hereafter, Paper I\] computed gas flows in the gravitational potential of the NIR luminosity distribution of Binney, Gerhard & Spergel , assuming constant NIR mass-to-light ratio ($`M/L`$). Their best SPH gas flow models reproduced quantitatively a number of observed gas flow features, including the positions of the five main spiral arm tangents at $`|l|60^{}`$ and much of the terminal velocity curve.

An important feature of all these gas flow models is the pattern speed of the non-axisymmetric component. Englmaier & Gerhard found a best pattern speed for the bar $`\mathrm{\Omega }_p60\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$. Weiner & Sellwood derived a bar pattern speed $`42\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$, whereas Fux determined $`50\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$ from his models. Dehnen  used resonant features in the Hipparcos stellar velocity distribution to argue that the Sun is located just outside the OLR of the exciting quadrupole perturbation, giving a pattern speed of $`51\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$ for solar constants $`R_0=8\mathrm{kpc}`$ and $`v_0=220\mathrm{km}\mathrm{s}^1`$. For their spiral structure model, Amaral & Lépine  and Lépine et al.  found $`\mathrm{\Omega }_{sp}20\mathrm{Gyr}^135\mathrm{Gyr}^1`$. Fernández et al.  obtained a somewhat higher $`\mathrm{\Omega }_{sp}30\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$ from Hipparcos data for OB stars and Cepheids. Debattista et al.  used the Tremaine-Weinberg method on a sample of intermediate age to 8 Gyr old OH/IR-stars in the inner Galactic disk. They found a pattern speed of $`59\pm 5\pm 10`$ (systematic) km/s/kpc, which may be driven by the bar in the center of the MW.

Thus the bar and spiral arms in the MW may not rotate with the same pattern speed. For a fast bar, a single $`\mathrm{\Omega }_p`$ would imply that the spiral arms are entirely outside their corotation radius. Observations of external galaxies (see, e.g., the Hubble Atlas of Galaxies, Sandage ) suggest that galaxies exist with dust lanes on the inner (concave) edges of their spiral arms. For a trailing spiral pattern, these arms would be inside their corotation radius. A lower pattern speed for the spiral structure than for the bar would remove this discrepancy. Indeed, Sellwood & Sparke  showed evidence for multiple pattern speeds in their $`N`$-body simulations. Rautiainen & Salo  analysed two-dimensional $`N`$-body simulations, some of them with a massless, dissipative gas component added. They confirmed the possibility of multiple pattern speeds in self-consistent N-body models of barred galaxies, and found a number of possible configurations. These included models with corotating bar and spirals, as well as models with different pattern speeds. Some of the models in the latter group show evidence for a non-linear mode-coupling (Tagger et al. ) between bar and spiral pattern, but others show no such evidence. In some of their models there exist separate inner spirals corotating with the bar, and outer spirals which rotate with their own, lower pattern speed. It is therefore tempting to analyse gas flow models of the MW with multiple pattern speeds.

What is the morphology of the MW spiral arms? Most authors infer four spiral arms from tracers which directly or indirectly measure the gas density, such as molecular clouds, HII regions, pulsars, and the galactic magnetic field (Georgelin & Georgelin , Sanders et al. , Caswell & Haynes , Grabelsky & et al. , Taylor & Cordes , Vallée ; however, Bash  infers a two-armed pattern from the same HII-data as used by Georgelin & Georgelin). The problem is that all spiral arm parameters other than the tangent point directions (e.g., Paper I) require distance information. It is also not clear whether all spiral arms seen in the MW gas are present in the old disk. Ortiz & Lépine  constructed a four-armed model which reproduces their star counts in the near infrared. Drimmel  preferred a two-armed structure from COBE/DIRBE K-band data, but a four-armed structure for the dust distribution seen in the $`240\mu \mathrm{m}`$ data. Drimmel & Spergel  project a luminosity model through the $`240\mu \mathrm{m}`$ dust model, to compare with the NIR J and K band COBE/DIRBE data. Their best model for the stellar distribution is four-armed, but dominated by two spiral arms. Drimmel & Spergel conclude that, if there are four arms in the K-band luminosity distribution, the Sag-Car arm is of reduced strength (by a factor of $`2.5`$).

In this paper we investigate the dynamical effects of the Galactic bar and spiral arms on the gas flow in the Milky Way. We investigate the possibility of different pattern speeds for bar and spiral arms, and the consequences this would have on the observed ($`l,v`$) diagrams. Our mass models for the inner Galaxy are based on the NIR luminosity density models of Bissantz & Gerhard  \[hereafter, Paper II\] which include spiral structure. We use SPH simulations to determine the gas flow in the MW, extending the work presented in Paper I, where eightfold symmetric mass models were used. This paper is organised as follows. In Section 2 we describe the luminosity models, methods, and observational data used in this work. In Section 3 we describe our best gas model for the observed 12CO ($`l,v`$) diagram. Then we compare models with different pattern speeds (Section 4), bar angles and spiral arm morphology (Section 5) to constrain these parameters, and finally give our conclusions in Section 6.

## 2 Gas dynamical model

### 2.1 Mass model of the Milky Way

Here we first describe the adopted model for the distribution of luminous mass in the MW. The model is based on the dust-corrected near-infrared maps of Spergel et al. , which they obtained from COBE/DIRBE data using a three-dimensional dust model derived from $`240\mu m`$ observations. From their L-band map we obtained a non-parametric luminosity distribution using the procedure described in Bissantz & Gerhard . Because the non-parametric model only covers the central $`5\mathrm{kpc}`$ part of the MW, we used a parametric best-fit model of the same L-band data to extend the model to larger radii. To convert the luminosity model to a model for the luminous mass density we assumed constant L-band mass-to-light ratio.

The model building procedure is described in detail in Paper II. Here we only summarise the most important points. Bissantz & Gerhard estimated the luminosity model iteratively from the L-band data by maximising a penalised likelihood function. The penalty terms encourage eightfold-symmetry with respect to the three main planes of the bar, as well as smoothness and a prescribed spiral structure model in the disk. The spiral structure term is based on an approximate model for the MW spiral arms from Ortiz & Lépine . The best models resulting from this approach reproduced the dust-corrected COBE/DIRBE maps with an RMS accuracy of 0.07 mag.

To define the eightfold symmetry penalty term we had to specify the position of the Sun in the MW. Bissantz & Gerhard set the distance of the Sun to the galactic centre to $`R_0=8\mathrm{kpc}`$, and the vertical distance from the galactic plane to $`z_0=14\mathrm{pc}`$. These values will be used throughout the present paper as well. The third parameter is the bar angle. Paper II compared photometric models for bar angles $`10^{}\phi _{\mathrm{bar}}44^{}`$, and concluded that the best models had bar angles $`20^{}\phi _{\mathrm{bar}}25^{}`$.

Paper II also showed that models which include spiral structure are better than models without spiral structure. In particular, because some spiral arm tangent points are evident in the L-band data, such models give a better description of the nearby disk. Also, the inclusion of spiral arms makes a model appear broader on the sky. For given data, this allows the bulge/bar to be more elongated in models with spiral structure. The larger bulge elongation in these models makes it possible to reproduce the asymmetries seen in the apparent magnitude distributions of clump giant stars in several bulge fields (Stanek et al. , ), for $`15^{}\mathrm{}<\phi _{\mathrm{bar}}\mathrm{}<30^{}`$. The shape of the bulge/bar in the model with $`\phi _{\mathrm{bar}}=20^{}`$ is about $`10:34:34`$, and its length is approximately $`3.5\mathrm{kpc}`$.

### 2.2 Gravitational potential

The gravitational force field and potential generated by the distribution of luminous mass were calculated using the multipole expansion method described in Paper I, modified to include phase terms in order to allow for the spiral arm components. The multipole expansion method was used because (i) it allows an independent treatment of disk, bar, and spiral arm components, (ii) it works for an arbitrary density distribution, and (iii) it does not require detailed boundary conditions as do other methods, since the integrations can easily be extended to infinite radii. The expansion was computed to maximum spherical harmonic orders $`l_{max}=8`$ and $`m_{max}=8`$. Odd orders were neglected, assuming point-symmetry with respect to the centre. Naively, one could interpret the $`m=2`$ component as due to the bar mode and all higher $`m`$ modes as due to spiral arm modes. However, some of the $`m=2`$ component beyond the bar’s corotation radius is due to an $`m=2`$ spiral arm mode in the luminosity model (see also Amaral & Lépine ).

Since the spatial resolution of the mass models obtained in Paper II is limited by that of the dust-corrected maps of Spergel et al., the central cusp in the MW’s density distribution and potential is incorrectly represented in our multipole expansion model. We attempted to correct for this by modifying the multipole coefficient functions in the centre as in Paper I, replacing the central mass distribution by a power law $`\rho ^{1.8}`$.

Optionally we add an analytical halo potential

$$\varphi _{\mathrm{Halo}}=\frac{1}{2}V_{\mathrm{inf}}^2\mathrm{ln}(r^2+a^2),$$

where $`V_{\mathrm{inf}}`$ is the circular velocity at infinity and $`a`$ is the core radius. In this paper we use $`V_{\mathrm{inf}}=220\mathrm{k}\mathrm{m}/\mathrm{s}`$ and determine the core radius such that the model in the potential with halo best fits the TVC near $`|l|=90^{}`$.

Figure 1 shows the rotation curve derived from the $`m`$=0–component of the reference luminosity density model of Paper II (for bar angle $`20^{}`$, with the parametric extension). For brevity, this mass density is called ”standard mass model” hereafter, and the associated gravitational potential ”standard potential”. Its rotation curve no longer shows a strong bump in the inner Galaxy’s rotation curve as for the models of Paper I (Fig. 4 there). This is because the new luminosity model reproduces light in the MW’s spiral arms significantly better, as compared to the density maxima $`3\mathrm{kpc}`$ down the minor axis of the bar in the earlier models. We note that the spiral arms in the luminosity models used here arbitrarily start at a galactocentric radius of $`3.5\mathrm{kpc}`$. This generates a slight distortion of the rotation curve around this radius.

Figure 2 shows the quadrupole and octopole terms of the potential, separated into the parts generated by the bar and the parts generated by the spiral arms. All multipoles are normalised by the value of the monopole term at the respective Galactic radius. Note that in this model the bar plays a major role for the non-axisymmetric forces, in particular the quadrupole moment, even beyond where it ends in the mass density.

The multipole representation of the potential is used in the orbital analysis and in the two-dimensional SPH hydrodynamics code, as described in Paper I. The bar and spiral components are both given a constant pattern speed. In most models described below, the bar and spiral arm patterns rotate with different pattern speeds, implying that the mass distribution and potential undergo periodic oscillations; in some models these pattern speeds are equal and the mass distribution is constant in the rotating frame.

Figure 3 shows a resonance diagram for the standard mass model with the same scaling as in Figure 1. (The scaling constant is different by $`\mathrm{}<1.5\%`$ for models 40 and 60, which have significantly different spiral pattern speeds from our standard model.) We will see below that the preferred bar pattern speed is $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$. Observe from Fig. 3 that in this case the corotation radius of the bar nearly coincides with the inner ultra-harmonic resonance of the spiral structure, if this has pattern speed $`\mathrm{\Omega }_{sp}=40\mathrm{G}\mathrm{y}\mathrm{r}^1`$, and with its inner Lindblad resonance, if $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$. Models with these values for the spiral arm pattern speed in conjunction with $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$ are particularly interesting to analyse (Tagger et al. ), and are discussed below.

### 2.3 Hydrodynamical method

We use the two-dimensional smoothed particles hydrodynamics \[SPH\] code described by Englmaier & Gerhard . The code solves Euler’s equation for an isothermal gas with effective sound speed $`c_s`$:

$`{\displaystyle \frac{𝐯}{t}}+(𝐯)𝐯=c_s^2{\displaystyle \frac{\rho }{\rho }}\mathrm{\Phi }.`$ (1)
This approach is based on a result of Cowie . He showed that an isothermal single fluid description crudely approximates the dynamics of the ISM. However, here the isothermal sound speed is not the thermal sound speed, but an effective sound speed representing the RMS random velocity of the cloud ensemble.

The SPH method has the advantage of allowing for a spatially adaptive resolution length. This is achieved by adjusting the smoothing length $`h`$ of a particle everywhere such that the number of particles that overlap a given particle is approximately constant. Fluid quantities are approximated by averaging over neighbouring particles. Furthermore, the SPH scheme includes an artificial viscosity to allow for shocks in the simulated gas flow. For further discussion of the method see Englmaier & Gerhard , Paper I, and Steinmetz & Müller .

Our SPH models contain $`5000060000`$ particles, except when indicated otherwise. The initial surface density of the models is taken to be constant inside $`8\mathrm{kpc}`$ galactocentric radius.

### 2.4 Comparison with observational data

Our main tools to compare our gas flow models with observations are the terminal velocity curve (TVC) and the ($`l,v`$) diagram, which shows the radial velocities of gas clouds as function of galactic longitude. Throughout this paper observed velocities are given with respect to the local standard of rest (LSR). For the distance of the Sun from the galactic centre we assume $`R_0=8\mathrm{kpc}`$. The LSR circular velocity is assumed to be $`v_0=220\mathrm{k}\mathrm{m}/\mathrm{s}`$, consistent with $`R_0=8\mathrm{kpc}`$ (Feast & Whitelock , Reid et al. , Backer & Sramek ). Note that a $`10\%`$ change in $`v_0`$ is not critical, amounting to radial velocity variations of only $`10\mathrm{km}\mathrm{s}^1`$ in the central $`l=45^{}`$ (cf. Paper I).

($`l,v`$) diagrams for our gas flow models are constructed as follows. The LSR observer is specified by the galactocentric radius $`R_0`$, the LSR circular velocity $`V_0`$, and by the angle $`\phi _{\mathrm{bar}}`$ relative to the bar. We first project all particle velocities $`(v_x,v_y)`$ onto the line-of-sight vector $`\stackrel{}{e}_p`$ from the LSR observer to the particle, subtracting the component of $`V_0`$ in the direction of the particle:

$$v_r=\stackrel{}{e}_p(v_x,v_y)v_0\mathrm{sin}(l).$$
(2)
We then construct a two-dimensional binned histogram of the particle distribution in the $`lv`$-plane, with bin size $`0.23^{^{}}\times 0.7\mathrm{km}/\mathrm{s}`$. Finally, we convert this histogram into a greyscale plot, using a lower surface density cutoff $`C_{\mathrm{lv}}=0.51\%`$ to enhance the contrast. $`C_{\mathrm{lv}}`$ varies between different models, and is selected so as to optimize the visibility of the spiral arm ridges and the terminal velocity envelope.

We compare a model with observations in a two step process. In the first step, we focus on the terminal velocities, comparing the model TVC with an observed TVC composed from the following data: HI velocities from Burton & Liszt ; HI velocities from Fich et al. , based on data from Westerhout ; unpublished $`140`$-ft single dish HI velocities, kindly provided by Dr. B. Burton; northern 12CO-velocities from Clemens , including error bars; and southern 12CO-velocities from Alvarez et al. . The Clemens  data were corrected for internal dispersion (Paper I). The observed velocities, corrected by the respective authors or Paper I to the pre-Hipparcos LSR frame, in which the Sun was assumed to move with approximately $`u_{}=10\mathrm{km}\mathrm{s}^1`$ inwards and $`v_{}=15\mathrm{km}\mathrm{s}^1`$ in the forward direction of Galactic rotation, are here transformed to the Hipparcos LSR frame ($`u_{}=10\mathrm{km}\mathrm{s}^1`$, $`v_{}=5\mathrm{km}\mathrm{s}^1`$, Dehnen & Binney ), by subtracting $`10\mathrm{sin}l\mathrm{km}\mathrm{s}^1`$.

A free parameter of our models is the mass-to-light ratio. From the comparison of the model TVC with the observed TVC we determine the best-fit scaling factor $`\xi `$ for the model velocities. This parameter $`\xi `$, which determines the mass scale of the model, varies by $`510\%`$ between all our models. We only use the TVC data for $`10^{}<|l|<50^{}`$ to determine $`\xi `$, because near the centre the resolution of our models is insufficient, and because for $`|l|\mathrm{}>50^{}`$ the halo contributes significantly (see below). In the fitting we take special care of the location of “bumps” in the TVC, because these indicate spiral arm tangents. This procedure assumes that the NIR disk and bulge are responsible for all of the observed velocities in the central parts of the MW.

In most models the gravitational potential is time dependent, because the pattern speeds of bar and spirals are different. In this case we select a “best” snapshot, corresponding to a specific phase and evolutionary age. The value of $`\xi `$ generally depends slightly on both the model and the evolutionary age of the snapshot.

In the second step we compare the model ($`l,v`$) diagram with the observed 12CO ($`l,v`$) diagram of Dame et al.  (Figure 4), using the scaling $`\xi `$ from the TVC. Important features in the observed ($`l,v`$) diagram are the ridges of emission which indicate the location of spiral arms. The white lines drawn in Fig. 4 reproduce approximately the locations of these spiral arm ridges. These lines are then transformed to the Hipparcos LSR \[Figure 4 was constructed assuming a solar motion of $`|𝒗_{}|=20\mathrm{km}\mathrm{s}^1`$ towards $`(l,b)=(56.2^{},22.8^{})`$), and are then overplotted on most model ($`l,v`$) diagrams. These observed ridges should be reproduced by a good model of the MW gas flow. The reverse need not always be true, however, because the visibility of a spiral arm ridge in the data may depend on the radial distribution of gas and the geometry of the arm with respect to the line-of-sight.

In the model ($`l,v`$) diagrams we also show molecular cloud observations of Dame et al.  and Bronfman  \[symbol “x” in the plots\], and HII region observations of Georgelin & Georgelin , Downes et al. , and Caswell & Haynes  \[symbol “+”\]. Most of our models do not contain a halo potential, and therefore underestimate the velocities for $`|l|\mathrm{}>4050^{}`$. Thus we omit the molecular cloud and HII observations in these models for $`|l|>40^{}`$. For the sake of clarity we have also left out clouds with less than $`10^{5.5}M_{}`$ from the Bronfman et al. data, as well as clouds in the smallest brightness bin from the Georgelin & Georgelin sample. For a more detailed discussion of these observations see Paper I.

## 3 Best-fit model for the Milky Way

We have investigated a number of gas flow models in the COBE potentials of Section 2, for different pattern speeds, bar angles, and stellar spiral arm morphologies. The analysis of these model sequences is deferred to Section 4. Here we begin with a description of our best model for the gas dynamics in the MW. This best-fit model (hereafter called ”standard model”) is based on the standard four-armed $`\phi _{\mathrm{bar}}=20^{}`$ luminosity model (Bissantz & Gerhard  – Paper II), as described in Section 2, it is point-symmetric, the bar pattern speed is $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$, implying corotation at $`R_{\mathrm{cr}}3.4\mathrm{kpc}`$, and the spiral arm pattern speed is $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$.

The gas distribution of the model is shown in Figure 5. In this plot the Sun is at $`(x,y)=(7.5,2.7)\mathrm{kpc}`$, with $`R_0=8\mathrm{kpc}`$. Both inside and outside corotation there exist four spiral arms which are connected in a complicated way through the corotation region of the bar. Outside corotation the spiral pattern in the gas response consists of a pair of strong arms and a pair of weaker arms.

To compare the model to the gas observations, we scale it to the observed terminal velocity curve (TVC) in the longitude range $`10^{}|l|50^{}`$, using an LSR circular velocity $`V_0=220\mathrm{km}\mathrm{s}^1`$. With this scaling the CO TVC is reproduced well by the model (Figure 6), including most distinct “bumps” in the observations except that at $`l50^{}`$. Vice-versa, the model TVC shows an extra bump at $`l15^{}`$ which is not seen in the data. Presumably the potential in this bar–disk transition region is not accurately modelled – there is little information in the NIR data on the mass distribution in this region. Overall, however, the fall-off with longitude and even the detailed form of the inner disk Galactic TVC are represented well by the model.

In the bulge region, for $`|l|\mathrm{}<10^{}`$, the model velocities are limited by the resolution of the hydrodynamic simulation. We thus do not expect to fit the observed large terminal velocities there, but several other effects may play a role as well. See §3.1 below and Englmaier & Gerhard  for a more detailed discussion. Because the gas particles in the numerical simulations flow inwards, however, the velocity structure in the inner $`11.5\mathrm{kpc}`$ is not of significant relevance for the gas flow and, particularly, for the shock structure in the main spiral arms well outside this region. Only the velocities of the pair of inner arms passing the minor axis of the bar laterally at $`1\mathrm{kpc}`$ could be somewhat affected. Since the main aim here is to investigate the large-scale gas flow and spiral arm morphology, we have therefore not attempted a detailed fit to the inner bulge terminal velocities.

The standard model does not contain a dark matter halo. For the assumed LSR velocity of $`V_0=220\mathrm{km}\mathrm{s}^1`$, its TVC falls below the observed TVC at $`|l|\mathrm{}>40^{}50^{}`$. At $`|l|=90^{}`$ the deficit in the model’s terminal velocity is $`25\mathrm{k}\mathrm{m}/\mathrm{s}`$. This can be corrected by a adding a halo potential that generates a rotation velocity $`v_{DM}120\mathrm{k}\mathrm{m}/\mathrm{s}`$ near the Sun. We thus constructed a new gas model in a potential which includes a suitable quasi-isothermal dark matter halo, which has circular velocity at infinity $`V_{\mathrm{inf}}=220\mathrm{k}\mathrm{m}/\mathrm{s}`$ and core radius $`a=10.7\mathrm{kpc}`$ (there is substantial freedom in these parameter values). The TVC of this model halo is also shown in Fig. 6; its best-fitting scaling factor $`\xi `$ and mass-to-light ratio is only slightly different from that of the standard model. Model halo fits the shape of the observed TVC out to $`R_0`$, but is not quite as good a match to the bumps in the TVC as the standard model. This suggests that the assumed halo model is oversimplified. Note the fact that the constant M/L standard model, a maximum disk model by construction, reproduces the observed terminal velocities inside galactocentric radius $`R5\mathrm{kpc}`$ and still accounts for most of the circular velocity near the Sun.

The ($`l,v`$) diagram for the standard model is shown in Figure 7. The dense ridges in this diagram show the locations of the spiral arms. Figure 11 of Paper II shows the correspondence between the locations of spiral arms in the Galactic plane and in the ($`l,v`$) diagram. The model spiral arm ridges in Fig. 7 generally coincide very well with the observed spiral arm ridges. Less good is the correspondence for the 3kpc-arm (again in the bar–disk transition region), which is at too small negative velocities in the model compared to the observations for $`l\mathrm{}>5^{}`$.

For model halo the lv-plot is shown in Figure 8, at enhanced contrast to emphasize the spiral arms (see §2.4). Compared to the standard model ($`l,v`$) diagram (without halo), significant differences are in the outer spiral tangents near $`l\pm 50^{}`$ which are relocated by a few degrees, and in the terminal velocities at $`|l|\mathrm{}>\pm 50^{}`$. This is a general result: For a number of models discussed later in this paper we have added a suitable halo potential and computed a new gas model, sometimes additionally changing the extent of the initial gas disk in the simulation from the standard $`8\mathrm{kpc}`$ to $`10\mathrm{kpc}`$. In all these cases the only significant change has been that the outer tangent points moved outwards in longitude by $`|\mathrm{\Delta }l|\mathrm{}<5^{}`$. Thus we do not include a halo potential in the remaining models discussed below.

Also shown in Fig. 8 are the observed spiral arm ridge lines, tangent points, and tracers from §2.4. The comparison with the grey scale plot for the model shows that model halo gives a very good description of the gas kinematics in the disk outside the bar. For many features deviations are less than $`10\mathrm{km}\mathrm{s}^1`$.

We end this section with showing in Figure 9 a map of the radial velocities of gas clouds with respect to the LSR for model halo. This map allows one to assess the likely errors made in determining kinematic distances from cloud radial velocities by assuming a circular orbit model.

### 3.1 Orbits and gas flow in the bulge region

In this section we consider in more detail the gas flow and TVC in the central $`|l|\mathrm{}<10^{}`$. The discussion of the large-scale morphology and pattern speed is continued in §4.

We begin with an analysis of closed orbits in the standard potential, including a central cusp in the mass model as described in Section 2.2, and assuming the same values for the pattern speeds as in the standard gas model. The closed orbits are found with a simple shooting algorithm. Despite the intrinsic time-dependence of the potential in the bar frame, the closed orbits in the inner kpc remain essentially unperturbed. Figure 10 shows closed x2\- and x1-orbits around the so-called ”cusped” orbit, the x1-orbit whose turning points on the bar’s major axis have a cusp shape. Closer to the galactic centre the more tightly bound x1-orbits become self-intersecting, and the x2-orbit family of stable orbits elongated perpendicular to the bar appears. Note that in this potential the outermost x2 orbits are (nearly) converging on the major axis of the bar, implying that gas clouds on these orbits would collide. This limits the radial extent of accessible x2-orbits for gas clouds in a hydrodynamic flow.

Between the cusped x1 orbit and the first non-intersecting x2 orbit there is a region in which no closed orbits suitable for gas flow exist (Fig. 10). Because of this, inflowing gas has to quickly pass this region: the mechanism described by Binney et al. , where gas moving in from the last $`x_1`$-orbits collides with gas on the outermost $`x_2`$-orbits, producing a spray that then forms an off-axis shock by hitting the far side of the $`x_1`$-orbits, cannot work as well in the potential here because of the lack of suitable outer $`x_2`$-orbits. Instead, the main place of dissipation of kinetic energy is likely to be the self-crossing loops of the $`x_1`$-orbits inside the cusped orbit, from where the gas moves inwards to hit the $`x_2`$-disk further in. This may explain why there is a mostly gas-free gap in the hydrodynamic simulations between the last non-intersecting x1 orbit and the first accessible x2 orbit, without any clear off-axis shocks like those found in other barred galaxy models (e.g. Athanassoula ).

How would gas clouds following the closed orbits in our standard model potential compare to the terminal velocity observations? Figure 11, plotted for bar angle $`20^{}`$, shows that the terminal velocities of the closed x1-orbits reproduce the observations surprisingly well. The cusped orbit in Fig. 11 is not only at the same longitude as the maximum in the observed TVC (at positive $`l`$ where we have data), but also appears to account for the decline of the observed terminal velocities at lower $`l`$. This is consistent with a gas-free gap between the cusped orbit and the first acceptable x2-orbit, due to which no strong leading shocks form like those observed in other barred galaxies. Episodic infall of gas clouds within the gap region may nonetheless form transient shocks similar to those observed by Hüttemeister et al. .

The closed orbits thus reproduce the observed terminal velocities, but why is there no gas at these velocities in the simulations? Possible explanations are as follows: First, the resolution length of our SPH code, which is dominated by the smoothing length of the SPH particles, may be too large to follow the strongly elongated x1-orbits near their cusped ends. Also, in the low density region further in, the relative velocities of neighbouring particles are large, so that the corresponding viscosity may lead to fast infall of SPH particles to the centre. This would depopulate the inner x1 orbits, where the largest terminal velocities are expected. Indeed, the terminal velocity curve attains larger peak velocities in high-resolution models with some $`10^5`$ particles (the standard gas model has $`610^4`$ particles), for example the maximum is at $`222\mathrm{k}\mathrm{m}/\mathrm{s}`$ on the $`l>0^{}`$ side, compared to $`210\mathrm{k}\mathrm{m}/\mathrm{s}`$ in the standard gas model.

Second, the detailed orbit shapes in the inner kpc of the MW depend on the fine-structure of the potential there, which is not accurately known, both because of the limited resolution of the underlying NIR data and luminosity model, and the difficulty of the deprojection in this region. The true x1 orbits could easily be somewhat less cuspy near their ends, or the $`x_2`$-orbits less converging, making them more easily populated with gas clouds.

Third, the hydrodynamic flow in the central MW could be genuinely slower than suggested by the closed orbits which it approximately follows. Increasing the model terminal velocities in the bulge region to match the observations would then require a somewhat larger bulge mass-to-light ratio than the value determined from the TVC fit at $`10^{}<|l|<50^{}`$. This might be plausible if the disk population is somewhat younger than that of the bulge. It is unlikely that the $`M/L`$ ratio is significantly larger in the entire central $`\mathrm{kpc}`$ of the MW because the peak velocitiy of gas on x2-orbits in our standard gas model is $`85\mathrm{k}\mathrm{m}/\mathrm{s}`$, which compares well to the observed $`\mathrm{}<80\mathrm{k}\mathrm{m}/\mathrm{s}`$ for CS-cloud cores at $`|l|\mathrm{}<0.7^{}`$, where the projected model x2-orbits are located. However, the potential near $`R1\mathrm{kpc}`$, which is most relevant for the x1 orbits, is more sensitive to the upper bulge component than that in the inner $`100\mathrm{pc}`$ where the x2 orbits are. The required variations of the bulge $`M/L`$ to give terminal velocities of $`\mathrm{}>250\mathrm{k}\mathrm{m}/\mathrm{s}`$ would seem consistent with the spread of K band mass-to-light ratios for other bulges. E.g., from NIR surface brightness photometry of early type spirals, Moriondo et al.  obtain a mean value and dispersion $`(0.6\pm 0.2)(M_{}/L_{})_K`$ \[the velocity scale of the standard gas model corresponds to $`(M/L)_K0.6(M_{}/L_{})_K`$\].

In summary, there are uncertainties in modelling the gas velocities in the inner few hundred parsec of the bulge. We will not pursue this further here, because as already discussed in the last subsection, this is not important for the gas flow in the main spiral arms further out. Rather, we now turn to the determination of the pattern speeds in the Milky Way.

## 4 The gas flow in models with separate bar and spiral arm pattern speeds

Figure 12 shows a time sequence for the gas distribution in our point-symmetric standard model with pattern speeds $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1,\mathrm{\Omega }_{\mathrm{sp}}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$. The inner arms emanating from the ends of the bar corotate with the bar pattern speed; they are clearly driven by the bar. These inner arms are connected to the outer spiral arms by a time-dependent transition region near bar corotation. Here a lateral arm from the distant end of the bar as seen from the Sun merges with the other inner arm from the nearer end in some frames (4-5), but not in others (1-2). In the latter case, the lateral arm continues into the outer spiral arm passing close to the Sun, in the former both arms join into a weaker outer arm staying well inside the Sun. The outer arms themselves move with respect to the bar frame. However, they do not rotate in the plot steadily around the model’s centre, as one might have expected if they were driven by the different pattern speed of the spiral arm potential, but apparently exhibit complicated back-and-forth oscillations, with respect to each other and with respect to the bar frame, and some arms merge and bifurcate at certain times; compare the vicinity of $`(x,y)=(5,3)`$ in the different panels.

The outer spiral structure evolves in such a complicated way because both the spiral arm potential and the bar simultaneously force the gas distribution with different pattern speeds. To investigate this further, we have computed a gas model similar to the standard model, but with spiral structure removed from the gravitational potential, so that the non-axisymmetric component of the potential is solely due to the bar. The resulting model has four spiral arms inside the bar corotation radius, and two symmetric outer spiral arms, all stationary in the bar frame. In the model with a driving spiral structure potential, the four outer arms can be regarded as a superposition of one component generated by the bar perturbation, which is at nearly constant position in the bar frame, and a second component driven by the spiral structure potential perturbation. The former is the stronger pair of arms in Fig. 12, one of which passes just inside the Sun symbol in the figure. The second component, visible as the weaker pair of arms between the bar-driven arms in Fig. 12, can be seen to fall behind the bar-driven arms through the sequence of frames in Fig. 12 ($`\mathrm{\Omega }_p>\mathrm{\Omega }_{sp}`$). At certain times the second component is seen to branch off the bar-driven arms (frame 1). Thereafter, the spiral-driven arm appears to fall behind, move inwards, until it finally collides with the opposite bar-driven arm (frames 2 through 6 $``$ frame 1 through frame 3). In the course of this evolution, the detailed morphology of the transition region around bar corotation changes. At most times, the spiral-driven arms are connected to the inner lateral arms (frames 2-5), at others the latter connect to the bar-driven arms (frame 1). At certain times, an arm may look fragmented, and its tangent point may split in longitude.

Thus when neither the bar nor the spiral structure component dominates the non-axisymmetric potential and gas flow in the disk, as in Fig. 12, it is difficult to deduce from a single snapshot of the arm morphology that the system supports a second, independent pattern speed. This may explain why it is difficult to say from the observed arm morphologies in barred spiral galaxies whether the spiral arms are driven by the bar or not (see Sellwood & Wilkinson ).

### 4.1 Pattern speeds in the Milky Way

We now investigate Galactic models with different combinations of the bar pattern speed $`\mathrm{\Omega }_p`$ and the spiral arm pattern speed $`\mathrm{\Omega }_{sp}`$. Particularly for the spiral arms the assumption of a constant pattern speed is probably still idealized, but represents a first step towards understanding realistic cases. In these models, the separation of bar and spiral arm components in the potential is based on the density distribution. We assume that the $`m2`$-multipoles of the density at galactocentric radii $`r<r_{\mathrm{cut}}=3.5\mathrm{kpc}`$ (bar) rotate with $`\mathrm{\Omega }_p`$, and outside of $`r_{\mathrm{cut}}`$ (spirals) with $`\mathrm{\Omega }_{sp}`$. The specific value of $`r_{\mathrm{cut}}`$ was chosen because it corresponds approximately to the end of the bar in the reference luminosity model of Paper II, and is also equal to the inner radius of the spiral pattern there. All models described in this Section are based on this luminosity model (cf. §2.1).

Multiple pattern speeds ($`\mathrm{\Omega }_p\mathrm{\Omega }_{sp}`$) imply a genuinely time-dependent potential in the frame corotating with the bar. For each model we have therefore investigated a sequence of snapshots at different evolutionary ages, separated by $`\mathrm{\Delta }\varphi 23^{}`$ in phase difference between the two components in the potential. This corresponds to steps in the evolutionary age of the model of $`0.010.02\mathrm{Gyr}`$, depending on the combination of pattern speeds. A finer analysis of our standard gas model in steps of $`\mathrm{\Delta }\varphi 9^{}`$ showed that no significant features in the models are missed with $`\mathrm{\Delta }\varphi 23^{}`$. Because the potentials used in this section are point-symmetric, we need to cover only a range of $`180^{}`$ in bar-spiral phase difference. All analysed snapshots have evolutionary ages at or around $`0.30\mathrm{Gyr}`$, the time after which the gas flow in the similar single pattern speed models of Paper I had become approximately quasi-stationary.

We begin with the three models 60, 40 and 20 with different spiral arm pattern speeds $`\mathrm{\Omega }_{sp}=61.4,40`$ and $`19.6\mathrm{Gyr}^1`$. The bar pattern speed is set to $`61.4\mathrm{Gyr}^1`$ ($`60\mathrm{G}\mathrm{y}\mathrm{r}^1`$ for model 20), based on the results of Paper I. Model 60 is the single pattern speed model; the two other values for $`\mathrm{\Omega }_{sp}`$ are motivated as follows. Firstly, Amaral & Lépine  preferred the lower value $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$, based on the positions and ages of open clusters, and from comparing their model with other data. Fernández et al.  obtained a somewhat higher $`\mathrm{\Omega }_{sp}30\mathrm{k}\mathrm{m}/\mathrm{s}/\mathrm{kpc}`$ from Hipparcos data for OB stars and Cepheids, supplemented with radial velocities or distances, respectively. Thus we selected $`\mathrm{\Omega }_{sp}=40\mathrm{G}\mathrm{y}\mathrm{r}^1`$ as an additional intermediate value between $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$ and the single pattern speed model. Secondly, for $`\mathrm{\Omega }_{sp}=60,20,40\mathrm{Gyr}^1`$, the bar corotation radius coincides approximately with the corotation, inner Lindblad, and inner ultra-harmonic (1:4) resonance of the spiral pattern, respectively (see Fig. 3).

What are the characteristic differences between these models? In Figure 13 we compare their spatial gas distributions. All three models have four inner spiral arms inside bar corotation (at $`3.4\mathrm{kpc}`$), and four outer spiral arms outside of a transition region beyond the end of the bar. Because we have always selected that snapshot of a model which best reproduces the observed ($`l,v`$) diagram (Fig. 4), the spiral arm tangent points are fixed relative to the position of the Sun. Thus in the figure the outer arms of all models appear approximately at the same positions.

However, near the bar corotation radius the models differ significantly. In the single pattern speed model 60, where the spiral structure corotation radius coincides with that of the bar, the gas in the corotation region moves with near-sonic velocities relative to the pattern. Consequently shocks are weak or nonexistent in this region, and the spiral arms in the gas response dissolve there. On the other hand, for model 40 with $`\mathrm{\Omega }_{sp}=40\mathrm{G}\mathrm{y}\mathrm{r}^1`$, the corotation radius of the spiral structure is at $`5\mathrm{kpc}`$. Again, the gas distribution shows gaps in the spiral arms at this location – one pair of arms nearly vanishes there, the other pair weakens – but there is no such gap in the spiral arms near bar corotation. Rather, the connection between the inner and outer spiral arms is dynamic, due to the different pattern speeds. Finally, in model 20 the spiral structure corotation radius occurs beyond the solar orbit. In this model we see no clear gaps in the arms. However, near the bar corotation radius $`R_{\mathrm{cr}}^{\mathrm{bar}}3.4\mathrm{kpc}`$, which coincides with the inner Lindblad resonance of the spiral pattern, the arms weaken, and the transition region appears to be more complicated than in model 40.

The comparison with the Galactic TVC is shown in Figure 14. All three models reproduce the observed data quite well, but the slope in the region $`10^{}|l|50^{}`$ is fit significantly better by the multiple pattern speed models than by the single pattern speed model. For $`|l|\mathrm{}>40^{}50^{}`$ the model TVCs start to fall below the observed TVC; as we have discussed in §3, for the assumed $`V_0=220\mathrm{km}\mathrm{s}^1`$ the dark halo starts to contribute to the observed velocities there, according to our NIR-based models.

Figure 15 shows ($`l,v`$) diagrams for the three models. To facilitate their interpretation, we have overplotted a number of observed features, and have highlighted particles belonging to specific arm features both in this figure and in the density plots of Fig. 13, labelled by ”A” and ”B”. Feature “A” corresponds to the 3kpc-arm in the observed CO ($`l,v`$) diagram (Fig. 4), feature “B” to its symmetric counter-arm in the models (see Fig. 13; in model 20 the neighbouring arm “A’ ” is included).

Particularly important is the existence of certain regions in the ($`l,v`$) diagram, e.g., next to the 3kpc arm (feature ”A”), where hardly any gas is found in the multiple pattern speed models. Such voids, designated by “V” in the plots, arise because the gas is aligned morphologically and kinematically with spiral arms there, which appear as well-defined ridges with adjacent voids in the ($`l,v`$) diagram. From Fig. 4 we see that similar voids are also visible in the observed 12CO ($`l,v`$) diagram, specifically next to features “A” and “B” in the bar corotation region. This suggests strongly that in the MW the gaseous arms go through bar corotation. On the contrary, in the single pattern speed model 60, the arms dissolve in the bar corotation region, the 3kpc-arm is thus imcomplete, and the gas is spread out approximately evenly over the corresponding parts of the ($`l,v`$) diagram. Models with a separate second pattern speed for the spiral arms in the MW are therefore preferred over single pattern speed models. Models with a growing bar amplitude also support spiral arms in the corotation region (Thielheim & Wolff ); however, it is likely that when self-gravity is included and the amplitude becomes non-linear, the growing spiral pattern will again develop an independent pattern speed.

Overall, the models in Fig. 13 provide a good match to the observed CO ($`l,v`$) diagram. The main features that are reasonably well represented are: the arm tangent at $`l=30^{}`$; the observed $`l=25^{}`$ tangent, although at smaller $`l20^{}`$ in models 20 and 40; the morphology of this arm; the location of the main spiral arm leading to the $`l=50^{}`$ tangent; the morphology of the ridges and voids at $`10^{}>l>25^{}`$ (apart from model 60); the arm morphology around $`(l,v)=(20^{},60\mathrm{km}\mathrm{s}^1)`$. The main weaknesses of the models, if we disregard the poorly resolved bulge region, are: the arm in the models which, returning from the $`l=20^{}`$ tangent, crosses the region $`5^{}\mathrm{}<l\mathrm{}<10^{}`$ at $`v>50\mathrm{k}\mathrm{m}/\mathrm{s}`$ and does not have a counterpart in the data, except perhaps if shifted to lower velocities; the missing tangent at $`l=30^{}`$ in models 60 and 20; the displacement of the 3kpc-arm towards lower velocities.

In the model 20, the spiral arm corresponding to the Centaurus tangent at $`l51^{}`$ is strongest when one of the spiral-driven arms coincides with the bar driven-arm (e.g., frames 3-4 in Fig. 12). Shortly before and after this evolutionary time this arm looks fragmented in the model, and its tangent point appears split in longitude. Indeed, there are observational indications that the Centaurus tangent is split into two parts at $`l50^{}`$ and $`55^{}(58^{})`$ (cf. Table 1 of Paper I).

The displacement of the 3kpc-arm in the model and the apparent absence of its counter-arm in the data might indicate a non-point-symmetric mass distribution in this region. This possibility is discussed further in Section 5.4. Also, in models where the spiral pattern speed differs from the bar pattern speed, we expect also the mass distribution in the transition region to be generally time-dependent. Thus in this region the mass distribution based on the NIR data can represent only one snapshot in time. In addition, there is not much information in the NIR data on the spiral arm heads in this crucial transition region, to constrain the luminosity model of Paper II. Hence the potential in this region is likely to be at best approximately correct.

In summary, the multiple pattern speed models reproduce the observed features in the bar corotation region better than the single pattern speed model, in particular the regions void of gas in the ($`l,v`$) diagramnear the 3 kpc arm. Comparing the lv-plots for models 20 and 40, we have found a slight preference for model 20 from the positions of the spiral arms; however, the differences between these two models are too small to determine the spiral arm pattern speed reliably. We take the model with $`\mathrm{\Omega }_{sp}20\mathrm{G}\mathrm{y}\mathrm{r}^1`$ as our standard model, and use this pattern speed in the following.

We now proceed to determine the best bar pattern speed $`\mathrm{\Omega }_p`$. To this end we consider two models bar50 and bar70 which have $`\mathrm{\Omega }_p=50\mathrm{G}\mathrm{y}\mathrm{r}^1`$ and $`70\mathrm{G}\mathrm{y}\mathrm{r}^1`$, respectively. For both models, we set $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$; since the differences between models 20 and 40 are small, the precise choice of $`\mathrm{\Omega }_{sp}`$ should not matter. Figure 16 shows ($`l,v`$) diagrams for the respective best evolutionary times of models bar50 and bar70. For both models the fit of the spiral ridges in the inner disk region of the ($`l,v`$) diagram is significantly worse than for the models with $`\mathrm{\Omega }_p60\mathrm{G}\mathrm{y}\mathrm{r}^1`$. In model bar50 the positions of the spiral arm ridges and tangents are worse (e.g., one of the $`l>0`$ tangents is missing), and the $`l<0`$ TVC and the 3kpc-arm are particularly badly represented. In model bar70 the spiral arm ridge with tangent position at $`l25^{}`$ is at significantly too small longitudes, and the same is true for the Centaurus tangent near $`l51^{}`$, which corresponds to one of the bar-driven outer arms (see discussion above). It is noteworthy that the position and shape of the Centaurus tangent are sensitive to $`\mathrm{\Omega }_p`$ despite a galactocentric radius of this tangent point of more than $`6\mathrm{kpc}`$. Note also that a pattern speed of $`70\mathrm{G}\mathrm{y}\mathrm{r}^1`$ would put corotation inside the end of the NIR bar (see §2 and Paper II).

We conclude that the best value for the bar pattern speed in the Milky Way is $`\mathrm{\Omega }_p=(60\pm 5)\mathrm{Gyr}^1`$. For the preferred scaling of the employed NIR bar model this corresponds to bar corotation at $`3.4\pm 0.3\mathrm{kpc}`$, equal to the length of the bar within the uncertainties. I.e., the Milky Way bar is a fast bar.

## 5 Bar orientation and influence of stellar spiral arms

### 5.1 Bar angle

In Paper II, luminosity models were generated from the COBE/DIRBE L-band data for bar angles $`\phi _{\mathrm{bar}}=10^{},15^{},20^{},25^{},30^{}`$ and $`44^{}`$ in the same way as for the standard $`\phi _{\mathrm{bar}}=20^{}`$ model. From considering the photometric residuals and the line-of-sight distributions of clump giant stars in the bulge, a preferred range for the bar angle $`15^{}\phi _{\mathrm{bar}}30^{}`$ was found. Here we obtain independent constraints on $`\phi _{\mathrm{bar}}`$ from corresponding gas flow models.

For each value of the bar angle, the luminosity model was converted to a mass model assuming constant M/L. As described above, the non-axisymmetric part of the potential was split into bar and spiral arm components, and gas models were computed with pattern speeds $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$ and $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$, respectively. In Figures 17 and 18 we show the TVCs and ($`l,v`$) diagrams of models incl10, incl15, incl25 and incl30 with bar angles $`\phi _{\mathrm{bar}}=10^{},15^{},25^{},30^{}`$ at their respective best evolutionary times, as judged by comparing their ($`l,v`$) diagrams with Fig. 4.

Model incl25 reproduces the terminal velocity observations with similar quality as the standard gas model 20, but model incl30 is clearly inferior – the terminal velocities are too large at positive longitudes and too small in modulus at negative longitudes. Obviously, this cannot be corrected for by a change of the velocity scaling factor $`\xi `$. The TVCs of models incl10 and incl15 are also inferior fits to the observed terminal velocities, albeit not as bad as model incl30.

The ($`l,v`$) diagrams of these models depend strongly on the assumed bar angle. Models incl10, incl15 and incl30 do not reproduce the observed spiral arm ridges. Models incl10 and incl15 also do not have a well-define 3 kpc arm at positive longitudes, and model incl30 shows very little arm structure at all in the region around bar corotation. Model incl25 reproduces the observations for a large range of longitudes similarly well as the standard model, but the positions of the spiral arm ridges at $`30^{}\mathrm{}<l\mathrm{}<10^{}`$ are not matched well and the non-circular velocities in the 3 kpc arm are smaller than in the standard $`\phi _{\mathrm{bar}}=20^{}`$ model. From the ($`l,v`$) diagrams model 20 is best, but we consider model 25 as still satisfactory.

One may ask, how much of the difference between these models is due to the different shapes of their underlying luminosity distributions and gravitational potentials, and how much of it is due simply to the different viewing geometries with respect to the gas flow? To answer this question, we have constructed the ($`l,v`$) diagram of the standard gas model 20 seen from a viewing angle of $`\phi _{\mathrm{bar}}=30^{}`$. The morphology of the arms in this ($`l,v`$) diagram differs little from that in the original $`\phi _{\mathrm{bar}}=20^{}`$ ($`l,v`$) diagram because of the tightly wound pattern. Thus the differences between the ($`l,v`$) diagrams in Fig. 18 must be caused mainly by genuine changes in the gas flow, originating from the different mass distributions corresponding to the COBE data for different $`\phi _{\mathrm{bar}}`$.

### 5.2 Spiral arm models

The MW very probably has four spiral arms in the gas distribution (e.g., Englmaier & Gerhard ). It is not clear, however, whether these also correspond to four stellar spiral arms, because some of the tangent points are not clearly seen in the near-IR light (see the discussion in the Introduction and in Drimmel & Spergel ). Here we investigate the ($`l,v`$) diagrams that result when different stellar spiral arm patterns drive the gas flow, and find evidence for a four-armed spiral pattern also in the distribution of luminous mass. Specifically, we have studied three models:

A two-armed model with the same spiral arm pitch angle $`13.8^{}`$ as in the standard (four-armed) mass model, but where the Sag-Car arm and its counter-arm are removed from the model. Such a model can be justified by the fact that the Sag-Car arm is hardly visible in the NIR.

Another two-armed model, but with approximately half the pitch angle of the standard four-armed model. This model reproduces approximately the same tangent point positions on the sky as the standard model.

Similar to the standard mass model, but the Sag-Car arm and its counter-arm are given only $`40\%`$ of the peak density amplitude of the other two arms. This is based on the result of Drimmel & Spergel , who found that they had to reduce the amplitude of the Sag-Car arm in their best-fit four-armed model for the COBE/DIRBE J and K-band NIR maps.

Luminous mass models with these patterns were derived using the algorithm described in Paper II, all for bar angle $`\phi _{\mathrm{bar}}=20^{}`$, by incorporating the respective spiral arm model in both the parametric initial model and in a penalty term for the non-parametric deprojection. As in the standard mass model the non-parametric algorithm changes the spiral arms somewhat, but the overall pattern stays intact. Because the NIR data constrain only the arm tangent points (Paper II), two-armed and four-armed models which reproduce the tangent point data fit the photometry with similar quality.

The TVCs of gas flow models computed in the corresponding gravitational potentials are shown in Figure 19. The overall slope with $`l`$ of all these model TVCs is similar, because this is dominated by the monopole term in the mass distribution. However, all three models open2, 2spi and mix do not fit the wavy structure of the TVC data as well as the standard model 20 (Fig. 6), with model mix the best among the three. For comparison we also include in Fig. 19 the TVC of a model noarms, whose gas flow was determined in a potential that includes only the perturbation from the bar, but not that from the spiral arms. This model has only two outer spiral arms in the gas distribution.

($`l,v`$) diagrams of these models are shown in Figure 20. Models open2 and 2spi compare poorly to the observed CO ($`l,v`$) diagram: the up-turning arm at $`l25^{}`$ is missing, and the spiral structure is generally wrong for $`15^{}\mathrm{}<l\mathrm{}<35^{}`$. In both these models the fit to the data is bad because the gas distribution of the model is only two-armed. We have checked that this is not a result of choosing the wrong bar and spiral arm pattern speeds, by computing gas models in the two-armed potential of model open2 for the additional combinations of $`(\mathrm{\Omega }_p=50\mathrm{G}\mathrm{y}\mathrm{r}^1,\mathrm{\Omega }_{\mathrm{sp}}=20\mathrm{G}\mathrm{y}\mathrm{r}^1)`$ and $`\mathrm{\Omega }_p=\mathrm{\Omega }_{sp}=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$ (single pattern speed). In both cases, the global, two-armed morphology of the gas flow is the same as in model open2. Model mix is the best of the three models in this section. However, the envelope of its lv-plot shows stronger bumps than the standard model, in particular near $`30^{}\mathrm{}15^{}`$. The gravitional potential of model mix is quite similar to our standard potential, and therefore the similarity of the gas flows is expected. From these tests we conclude that a four-armed spiral arm potential is preferred (see also Englmaier & Gerhard  and Fux ).

### 5.3 Stronger spiral arms

The standard gas model fits well most spiral arms in the observations. However, the famous 3kpc-arm, a prominent feature which extends from $`(l10^{},v=0\mathrm{k}\mathrm{m}/\mathrm{s})`$, through $`(l=0^{},v50\mathrm{k}\mathrm{m}/\mathrm{s})`$, to $`(l22^{},v120\mathrm{k}\mathrm{m}/\mathrm{s})`$ is still displaced towards lower non-circular velocities. Because the underlying mass model already contains spiral arms, assigning gravitating mass to the gas particles does not lead to an improvement (as a corresponding model confirmed).

However, compared to near-IR-observations of other spiral galaxies (cf. Rix & Zaritsky ), our standard mass model has rather narrow and weak spiral arms. This might be simply because only the positions of the spiral arm tangents are constrained by the COBE/DIRBE L-band data, while the spiral arm heads near the bar are not well-constrained. It is therefore possible that our standard mass model underestimates the strength of the Galactic spiral arms. We have therefore investigated a model strongarms in which we enlarged non-axisymmetric forces in the disk, by multiplying with $`1.5`$ the $`m2`$ multipoles of the spiral arm component (i.e., the component corresponding to the density outside of $`r_{\mathrm{cut}}=3.5\mathrm{kpc}`$) in the standard model, and computed the gas flow in this modified potential.

We show in Figure 19 the TVC and in Figure 21 the ($`l,v`$) diagram of this model strongarms. The fit to the TVC observations is not as good as for the standard model 20, but in the ($`l,v`$) diagram the spiral arm ridges are quite similar. Due to the stronger spiral arm gravity, the bumps in the terminal velocity are somewhat stronger, but still in the acceptable range. Interestingly, the 3kpc-arm now fits the observations nearly perfectly. This shows that the strength of the spiral arms is an important parameter for the observed kinematics of the 3 kpc arm, and the non-circular motions in this region of the Galactic ($`l,v`$) diagram cannot simply be used to determine the bar aspect angle (cf. Weiner & Sellwood 1999). It also suggests that our standard model can be improved when a better spiral arm model becomes available.

### 5.4 Asymmetric Models

There are indications in the HI and CO surveys that the inner Galaxy’s gas distribution deviates significantly from point-symmetry with respect to the centre. An example is the 3kpc-arm, which has no clear counter-arm in the observed ($`l,v`$) diagrams. That is not to say that there is no counter-arm; if asymmetric, its inner parts could for example appear at similar locations in Fig. 4 as the arm which reaches the TVC at $`l25^{}`$. All symmetric mass models, however, yield counter-arms with about the same absolute velocities as the 3kpc-arm. A way out of this dilemma has been shown by Fux . In his model, the 3kpc-arm and its counter-arm are significantly disturbed by strong non-axisymmetric modes in the Galactic centre, as well as in the outer disk. Such a mechanism may also help to explain that the peak in the terminal velocity at $`l=+2^{}`$ appears much higher than at $`l=2^{}`$ in the CO data. In addition, the x2-disk in the centre may be disturbed by this effect, and this might explain the uneven gas distribution seen in CS. There is no evidence for an asymmetric mass distribution in the NIR data, so in our mass models we can only introduce asymmetry in the density by hand and study the consequences.

To this end, we first created uneven $`m`$ modes in the disk by weakening one or two arms in the initial model, or by moving two spiral arms closer together. In such gas models we observed strong effects on the position of gas shocks and their relative strength in the outer disk. Some cases look similar to the result of Fux with an almost 3-armed outer disk structure. The inner arms, especially the 3kpc-arm, did not change much, however. Only in an extreme case have we been able to move the 3kpc-arm to higher velocities, while simultaneously the counter-arm was moved to lower velocities, but this model does not fit the observations well overall.

In a second class of asymmetric models, we let the stellar bar centre rotate on a circular orbit with radius $`R_{\mathrm{bar}}`$ and pattern speed $`\mathrm{\Omega }_c=60\mathrm{G}\mathrm{y}\mathrm{r}^1=\mathrm{\Omega }_\mathrm{p}`$, i.e., the centre rotates backwards with respect to the bar and with the same pattern speed as the bar. This introduces also a third parameter, the phase $`\alpha _{\mathrm{bar}}`$ of the bar centre rotation at model age $`0.00\mathrm{Gyr}`$. This approach was motivated by $`N`$-body simulations (Debattista, work in progress). In simulation tumblingbar ($`\mathrm{\Omega }_c=60\mathrm{G}\mathrm{y}\mathrm{r}^1,\mathrm{R}_{\mathrm{bar}}=800\mathrm{pc},\alpha _{\mathrm{bar}}80^{}`$), we see the 3kpc-arm and its counter-arm move in the right direction in the model’s ($`l,v`$) diagram (Figure 22) and the $`3\mathrm{kpc}`$-arm fits the observations well. However, we have not found a snapshot of a model at which the 3kpc-arm, its counter-arm and the overall spiral pattern all fit the data well.

Obviously, the available freedom in introducing deviations from point-symmetry is very large. We have only tested a few attractive possibilities, and these models show that asymmetries in the Galactic mass distribution may be important and could be at the root of some of the remaining problems in our standard model.

## 6 Summary and Conclusions

We have used new gas flow models to investigate the dynamics of the Milky Way (MW) Galaxy from observed ($`l,v`$) diagrams. Steady-state gas flows in rotating, point-symmetric gravitational potentials for the Galactic bar and disk were determined with SPH simulations. The potentials were derived from non-parametric estimates of the spatial near-infrared luminosity density (Bissantz & Gerhard ), based on the de-reddened COBE/DIRBE L-band map of Spergel et al. , but also incorporating clump giant star count data from Stanek et al. , . The luminosity models contain a spiral arm model for the disk, and were converted to mass models assuming constant mass-to-light ratio in the inner MW.

Our best gas flow model gives a very good fit to the Galactic terminal velocity curve for $`|l|>15^{}`$, and to the spiral arm ridges in the observed CO ($`l,v`$) diagram. This has enabled us to investigate a number of dynamically important parameters such as the bar and spiral arm pattern speeds, the multiplicity of the spiral structure in the potential, and the bar angle. The main results from this study are as follows.

1) In gas flow models with separate pattern speeds $`\mathrm{\Omega }_p`$ for the bulge/bar and $`\mathrm{\Omega }_{sp}`$ for the spiral pattern, the spiral arms go through the bar corotation region. Thus ($`l,v`$) diagrams of such models show well-defined spiral arm shocks (ridges) through corotation, next to areas which appear to be nearly void of gas. By contrast, in single pattern speed models the spiral arms dissolve in the bar corotation region, so that the gas fills this region of the ($`l,v`$) diagram approximately evenly, and no voids exist.

2) Similar voids are visible in the observed 12CO ($`l,v`$) diagram. From a comparison with model ($`l,v`$) diagrams with different $`\mathrm{\Omega }_{sp}`$ but similar $`\mathrm{\Omega }_p`$ we find evidence for separate pattern speeds in the Milky Way. The existence of self-consistent models with separate bar and spiral arm pattern speeds was demonstrated by Rautiainen & Salo  in a study of two-dimensional N-body simulations, some of which included a massless, dissipative gas component. Models with a growing bar amplitude also support spiral arms in the corotation region (Thielheim & Wolff ); however, it is likely that when self-gravity is included and the spiral arm amplitude becomes non-linear, the growing spiral pattern will again develop an independent pattern speed.

3) From a series of models the preferred range for the bar pattern speed in the MW is $`\mathrm{\Omega }_p=60\pm 5\mathrm{G}\mathrm{y}\mathrm{r}^1`$, corresponding to corotation at $`3.4\pm 0.3\mathrm{kpc}`$. This agrees well with previous pattern speed determinations by Englmaier & Gerhard , Dehnen , and Debattista et al. . The bar pattern speed is well constrained because it influences not only the inner spiral structure, but also the position of two outer spiral arms in the lv-plot. Models with $`\mathrm{\Omega }_p=50\mathrm{G}\mathrm{y}\mathrm{r}^1`$ and $`\mathrm{\Omega }_p=70\mathrm{G}\mathrm{y}\mathrm{r}^1`$ are inferior.

The spiral arm pattern speed is less well constrained. Our preferred value is $`\mathrm{\Omega }_{sp}20\mathrm{G}\mathrm{y}\mathrm{r}^1`$, but models with larger $`\mathrm{\Omega }_{sp}<\mathrm{\Omega }_p`$ give only marginally inferior fits to the observed ($`l,v`$) diagram.

4) Gas flows in models which include massive spiral arms clearly fit the observed 12CO ($`l,v`$) plot better than if the potential does not include spiral structure. Furthermore, comparing models with two and four arms in the gravitational potential, we found that only models with four massive arms reproduce the Galactic ($`l,v`$) diagram, while gas flows in two-armed potentials do not resemble the spiral arm pattern of the Milky Way.

In Galactic models with four-armed potentials and separate spiral arm pattern speed, the gas flow has two pairs of inner arms which rotate with the bar (lateral, and corresponding to the 3 kpc arm), and four outer spiral arms which exhibit a complicated, time-dependent back–and–forth oscillation in the bar frame. The outer and inner spiral structures are connected by a time-dependent transition region around bar corotation.

5) From a further series of gas models computed for different bar angles, using separately determined luminosity models and gravitational potentials as in Bissantz & Gerhard , we found a range of acceptable bar angles $`20^{}\mathrm{}<\phi _{\mathrm{bar}}\mathrm{}<25^{}`$. The models for $`\phi _{\mathrm{bar}}=15^{}`$ and $`30^{}`$ are clearly inferior, which is mainly due to differences in the inferred gravitational potential.

The model with ($`\phi _{\mathrm{bar}}=20^{}`$, $`\mathrm{\Omega }_p=60\mathrm{G}\mathrm{y}\mathrm{r}^1`$, $`\mathrm{\Omega }_{sp}=20\mathrm{G}\mathrm{y}\mathrm{r}^1`$) gives an excellent fit to the Galactic terminal velocity curve for $`10^{}|l|50^{}`$, and to the gaps and spiral arm ridges in the observed CO ($`l,v`$) diagram. There are still discrepancies in the bar corotation region where the potential has uncertainties: the 3 kpc arm has too low non-circular velocity, and its counterarm is missing in the data. In the bulge region, closed orbits reproduce the TVC well, while the gas model has lower velocities. This may be a resolution problem in the SPH model, but could in part also be due to uncertainties in the potential which influence the orbit shapes.

6) The 3 kpc arm non-circular velocities can be reproduced by a model in which we artifically increased the $`m2`$-multipoles of the spiral potential component by $`1.5`$ while keeping all other dynamical parameters fixed. This is well within the uncertainties. Guided by a number of asymmetries in the observed Milky Way gas distribution, we also investigated potentials which are no longer point-symmetric. Some of these models improved the fit to the 3kpc-arm and its counter-arm. Although we did not find a model which at the same time reproduces the entire ($`l,v`$) diagram as well as our standard model, these models suggest that such asymmetries may be important for better understanding the gas flow in the inner Milky Way.

## Acknowledgments

This work was supported by grant 20-64856.01 of the Swiss National Science Foundation.


format: text

text:
# Negative refraction and left-handed behavior in two-dimensional photonic crystals

## Abstract

We systematically examine the conditions of obtaining left-handed (LH) behavior in photonic crystals. Detailed studies of the phase and group velocities as well as the phase $`n_p`$ and group $`n_g`$ refractive index are given. The existence of negative refraction does not guarantee the existence of negative index of refraction and so LH behavior. A wedge type of experiment is suggested that can unambiguously distiguinsh between cases of negative refraction that occur when left-handed behavior is present, from cases that show negative refraction without LH behavior.
PACS numbers: 78.20.Ci, 41.20.Jb,42.25.-p, 42.30.-d

Recently, there have been many studies about materials that have a negative index of refraction $`n`$ smithpbg ; smithprl ; shelbyapl ; shelbysc ; mehmet . These materials, theoretically discussed by Veselagovesalago , have simultaneous negative permittivity $`ϵ`$ and permeability $`\mu `$. Pendry pendry ; pendry00  has studied structures that they can have both $`ϵ`$ and $`\mu `$ negative and suggested that they can be used to fabricate a perfect lens. In these materials, $`\stackrel{}{k}`$, $`\stackrel{}{E}`$, $`\stackrel{}{H}`$ form a left-handed (LH) set of vectors (i.e., $`\stackrel{}{S}\stackrel{}{k}<`$ 0, where $`\stackrel{}{S}`$ is the Poynting vector) and, therefore, are called left-handed materials (LHM). Futhermore Markoš and Soukoulis markos  used the transfer matrix technique to calculate the transmission and reflection properties of left-handed materials. These scattering data were used smithprb  to determine the effective permitivity and permeability of the LHM and the index of refraction was found smithprb  to be unambiguously negative in the frequency region where both of the retrieved values of $`ϵ`$ and $`\mu `$ are negative. All experiments smithpbg ; smithprl ; shelbyapl ; shelbysc ; mehmet  that showed left-handed behavior were performed in the microwave or millimeter regime.

In addition Notomi Notomi  studied light propagation in strongly modulated two dimensional (2D) photonic crystals (PCs). Such a photonic crystal behaves as a material having an effective refractive index $`n_{\mathrm{eff}}`$ controllable by the band structure. In these PC structures the permittivity is periodically modulated in space and is positive. The permeability $`\mu `$ is equal to one. Negative $`n_{\mathrm{eff}}`$ for a frequency range was found. The existence of negative $`n`$ was demonstrated Notomi  by a finite difference time domain (FDTD) simulation. Negative refraction on the interface of a 3D PC structure has been experimentally observed by Kosaka et al. kosaka . A negative refractive index associated to the negative refraction was reported. Also in Ref. 12, large beam steering was observed that was called the “superprism phenomenon”. Similar unusual light propagation was observed Russel  in 1D and 2D diffraction gratings. In addition, in Luo et. al. it was shown that the phenomenon of negative refraction in the PC can be utilized to make a superlensLuo . Finally a theoretical work Gralak  predicted a negative refraction index in PCs.

Both in the PC and LHM literature there is a lot of confusion about what are the correct definitions of the phase and group refractive index and what is their relation to negative refraction. It is important to examine if left-handed behavior can be observed in photonic crystals at optical frequencies. In addition, it is instructive to see how the latter relates with the sign of the phase and group refractive index for the PC system.

In this paper we will attempt to clarify the differences between phase refractive index $`n_p`$ and group refractive index $`n_g`$ in PC structures. We show that left-handed behavior in photonic crystals is possible and a comparison of this LH behavior to regular LH materials is given. We will also show that a frequency range exists for which negative refraction can occur while it does not correspond to regular LH behavior. A wedge type of experiment is suggested that can unveil the sign of $`\stackrel{}{S}\stackrel{}{k}`$note  (i.e. the “rightness”) in the PC and so the correct sign for the “effective” index (phase index $`n_p`$ as we will see in the following). This way we can unambiguously distinguish between the case of “negative refraction” that corresponds to LH behavior in the PC structure ($`n_p`$ and $`\stackrel{}{S}\stackrel{}{k}`$note  are negative), from cases that show “negative refraction” without LH behavior ($`n_p`$ and $`\stackrel{}{S}\stackrel{}{k}`$ are positive). FDTDtavlove  simulations are presented which support our claims that left-handed behavior can be seen in the the PC only in the cases where the predicted value for $`n_p`$ (for the infinite system) is negative. In all FDTD simulations we used perfectly matched layer (PML) Berenger  boundary conditions.

Our results are general but for completeness and for comparison with the results of Notomi Notomi  we studied a 2D hexagonal lattice of dielectric rods with dielectric constant 12.96 (e.g., GaAs or Si at 1.55 microns) for the H (TE) polarization. The radius of the dielectric rods is $`r=0.35a`$, where $`a`$ is the lattice constant. In Fig. 1a the band structure of the system is shown. In the figure as well as in the following the quantity $`\stackrel{~}{f}`$ is the frequency normalized as $`\omega a/2\pi c`$, where $`c`$ is the velocity of light. To study the system the surfaces (contours in the 2D system) of equal frequency in $`\stackrel{}{k}`$-space,– called equifrequency surfaces (EFS)–, are needed. The EFS basically consists of the allowed propagation modes within the crystal for a certain frequency. For both the band structure and the EFS plots, the plane wave expansion method chan ; villen  with 1003 plane waves was used. The EFS can have different shapes for various frequencies. We notice that in some cases their shape can be almost circular (i.e., $`\omega Ak`$, where $`A`$ a constant). Basically, these cases are for frequencies close to the band edges and we will refer to their dispersion relation $`\omega (k)`$ as normal dispersion from now on.

Before we proceed in presenting the FDTD simulation results it is useful to discuss the appropriate expressions for the phase ($`\stackrel{}{v}_p`$) and group ($`\stackrel{}{v}_g`$) velocities and their associated phase ($`n_p`$) and group ($`n_g`$) refractive indices for the infinite photonic crystal system. The analysis that follows is not used in the interpretation of the FDTD results, that we discuss later in this paper. However it does serve as a theoretical prediction to be compared with the FDTD data. For any general case $`\stackrel{}{v}_p=\frac{c}{|n_p|}\widehat{k}`$ , with $`\widehat{k}=\stackrel{}{k}/k`$ . Also $`|\stackrel{}{v}_g|`$=$`|_k\omega |`$=$`\frac{c}{|n_g|}`$. The value of $`|n_p|`$ for a certain angle of incidence $`\theta `$ will be $`c|\stackrel{}{k}_f(\theta )|/\omega `$ (the refracted wave-vector $`\stackrel{}{k}_f`$ can be determined from the EFS as it will be discussed later). The sign of $`n_p`$ is determined from the behavior of EFS. If the equifrequency contours move outwards with increasing frequency then $`\stackrel{}{v}_g\stackrel{}{k}>`$0; if they move inwards $`\stackrel{}{v}_g\stackrel{}{k}<`$0. It can be proven analytically that for the infinite PC system the group velocity coincides with the energy velocity fotein ; sakoda . Therefore the sign of $`\stackrel{}{v}_g\stackrel{}{k}`$ is equivalent to the sign of $`\stackrel{}{S}\stackrel{}{k}`$note  – where $`\stackrel{}{k}`$ is in the $`1^{st}`$ BZ–. So, to be in accordance with the LH literaturevesalago  the sign of $`n_p`$ will be the sign of $`\stackrel{}{v}_g\stackrel{}{k}`$. In the special case of normal dispersion, $`n_p`$ is independent of the angle of incidence . For the cases with dispersion close to normal $`n_p`$ vs. $`\stackrel{~}{f}`$ for up to the $`6^{th}`$ band has been calculated and the results can be seen in Fig. 1b. We observe in Fig. 1b that $`|n_p|`$ can have values less than one (i.e., superluminal phase velocity). However $`|n_g|`$ has to be always greater than one, because the group velocity is the energy velocity in the PC, which is less than $`c`$. We have calculated the group velocity with the $`\stackrel{}{k}\stackrel{}{p}`$ perturbation method Busch  and checked that indeed $`|n_g|>`$1 . For the cases with normal dispersion $`\stackrel{}{v}_g`$=$`\frac{c}{n_g}\widehat{k}`$ with $`n_g`$=$`\omega d|n_p|`$/$`d\omega +|n_p|`$ proof . From this relation and from the graph of $`n_p`$ versus frequency (see Fig. 1b), it can be checked that indeed, $`n_g>`$0 when $`n_p>`$0 and $`n_g<`$0 when $`n_p<`$0. This means that for the normal dispersion cases negative refraction occurs when the phase index is negative. We note this finding is different of that stated in Ref. 26. For the anisotropic dispersion cases $`\stackrel{}{v}_g`$ and $`\stackrel{}{k}`$ are not along the same direction. In these cases we choose the sign of $`n_g`$ to manifests the sign of refraction at the interface and therefore $`n_g`$ will have the sign of $`(\stackrel{}{v}_g\stackrel{}{k}_i)_x`$ (x being the direction along the interface and $`\stackrel{}{k}_i`$ the incident wavevector). We would like to stress that $`n_g`$ cannot be used in a Snell-like formula to obtain the signal’s propagation direction in the crystal. However, $`n_p`$ can be used in Snell’s law to determine the signal’s propagation angle in the normal dispersion case, under certain conditions fotein .

There are two characteristic cases for which an EM wave refracts in the “wrong way” (negatively) when it hits the PC interface. One can predict the direction of the refracted signal from the equifrequency contours in the air and PC media. This methology has been used widely in 1D diffraction gratingsconfined . The EFS for two frequencies, –$`\stackrel{~}{f}=0.58`$ and $`\stackrel{~}{f}=0.535`$ – that are representative for these two cases, are seen in Fig. 2a and Fig. 2b. The $`\stackrel{}{k}_{//}`$ conservation condition given by the dashed vertical line in the figures determines the allowed refracted wavevector $`\stackrel{}{k}_f`$. In both cases there are two choices for $`\stackrel{}{k}_f`$ (to point either towards A or towards B). For every choosen $`\stackrel{}{k}_f`$, the corrresponding $`\stackrel{}{v}_g`$ will be perpendicular to the EFS at that point and point towards increasing values of frequency. The correct choice for $`\stackrel{}{k}_f`$ will be the one that gives a $`\stackrel{}{v}_g`$ that points away from the source. After checking the EFS for a region around $`\stackrel{~}{f}=0.58`$ and $`\stackrel{~}{f}=0.535`$ we determine that the frequency increases inwards for the case of Fig. 2a (i.e. $`\stackrel{}{v}_g\stackrel{}{k}<0`$ and $`n_p<`$0 ) while it increases outwards for the case of Fig. 2b (i.e. $`\stackrel{}{v}_g\stackrel{}{k}>0`$ and $`n_p>0`$). Taking this into consideration the correct choice for $`\stackrel{}{k}_f`$ is to point towards A for the case of Fig. 2a and towards B for the case of Fig. 2b. For an infinite PC the group and the energy velocity are equal, so the group velocity vector (orange vector) in the figures represents the direction of propagation for the EM signal in the PC. We have confirmed the theoretical prediction of negative refraction with a FDTD simulation for both cases. In other words, we have determined both theoreticaly and numerically that the group index $`n_g`$ for both the cases is negative.ngnote . However this type of FDTD simulation cannot decide what is the “rightness” of the system.

Therefore we consider an incoming EM wave incident normally on a wedged structure, as was the case of the UCSD experiment shelbysc . We should mention that couplings to higher order Bragg waves sakoda  when the wave hits the wedged interface are unavoidable in some cases. However, the sign of the angle of the $`0^{th}`$ order outgoing Bragg wave coincides with the sign of $`\stackrel{}{S}\stackrel{}{k}`$ inside the wedge, and determines the “rightness”vesalago  (positive for the right and negative for the left handed case) and so the sign of $`n_p`$ of the PC.

In particular, we show the results for this type of simulation for $`\stackrel{~}{f}=0.58`$ (case of Fig. 2a) and $`\stackrel{~}{f}=0.535`$ (case of Fig. 2b) in Figs 3b and 4b respectively. For each of these cases the symmetry directions of the first and the wedged interface are choosen appropriately fotein . Naturally, the $`k_{//}`$ wavevector will be conserved across both interfaces. In Fig. 3b the final outgoing beam is in the negative hemisphere (negative in respect to what is expected for a regular homogeneous material with positive index). This means that the perpendicular to the first interface component of the wavevector has reversed sign when crossing the air PC interface. Since $`\stackrel{}{S}`$ is always pointing away from the source $`\stackrel{}{S}\stackrel{}{k}`$ will be negative in the PC for the case of Fig. 3b. We note that the second outgoing beam in Fig. 3b is due to consequent multi-reflection on the upper corner (one can determine that by looking at the field inside the wedge). To extract the sign of $`\stackrel{}{S}\stackrel{}{k}`$ in Fig. 4b, because of the presence of the $`2^{nd}`$ beam (higher order Bragg wave proof2 ) we need the additional knowledge of the theoretically expected magnitude of the angle for the $`0^{th}`$ order Bragg wave. We use the EFS contours for the infinite system to determine the magnitude of the refracted wavevector and from that the magnitude of the $`0^{th}`$ order Bragg wave. (We obtained a value of $`57^0`$). Although the actual magnitude of the refracted wavevector may differ from the predicted one, no dramatic difference is expected. We then observe in the FDTD simulation of Fig. 4b that the beam with deflection angle of magnitude close to $`57^0`$ lies in the positive hemisphere. So, in the case of Fig. 4b $`\stackrel{}{S}\stackrel{}{k}>0`$.

The ray schematics in Figs 3a and 4a – that accompany the simulations of Figs. 3b and Fig. 4b– show the theoretically predicted refracted and final outgoing beams by using the properties of the infinite system. The beams are determined in a similar manner as in the ray schematics of Fig. 2. In the following we compare the theoretical expected angles (magnitude and sign) with the ones measured in the FDTD simulation. The expected outgoing angle for the case of Fig. 3b is $`36^0`$, while the FDTD one is $`30^0`$. Likewise for the case of Fig. 4 the formula in Ref. 28 predicts an outgoing wave at $`57^0`$ and a $`1^{st}`$ order Bragg wave at $`14^0`$. These values are close to those of $`60^0`$ and $`10^0`$ measured in the FDTD simulation.

By looking on other bands we can find cases with $`n_p>0`$ and $`n_g>0`$, as well as $`n_p<0`$ and $`n_g>0`$ , so all combinations of signs for the group and phase index are possible. By performing FDTD simulations (of both types (slab and wedge)) we have shown that in all cases, that the “rightness” of the PC – i.e. the sign of $`\stackrel{}{S}\stackrel{}{k}`$ note –, coincides with the the theoretical prediction for the sign of $`n_p`$ , –as defined earlier in the paper for the infinite PC–, and is irrelevant to the sign of $`n_g`$. For the case of Fig. 2a, the PC is left-handed and the perpendicular component of the wavevector reverses sign in the air PC interface. In addition it has close to normal dispersion and a negative effective index. So one could say that the negative refraction in this case is similar to the one occurring in a regular LH material that has both $`ϵ`$ and $`\mu `$ negative.smithpbg ; smithprl ; shelbyapl ; shelbysc ; mehmet ; vesalago ; pendry ; pendry00 ; markos ; smithprb

In conclusion, we have systematically examined the conditions of obtaining left-handed behavior in photonic crystals. We have demonstrated that the existence of negative refraction is neither a prerequisite nor guarantees a negative effective index and thus LH behavior. Although we have focused on the presentation of 2D photonic crystals, the results are general enough and are expected to be correct for 3D PCs. Special care must be taken for the couplings of the TE modes and TM modes in the 3D systems. We have suggested experiments to check our theoretical predictions both in the microwave regime as well as in the optical regime. The present work will help in the understanding and design of the optical devices.

###### Acknowledgements.

We thank K. Busch for useful discussions regarding the group velocity calculation with the $`\stackrel{}{k}\stackrel{}{p}`$ perturbation method. We would also like to thank Mario Agio for useful discussions. Ames Laboratory is operated by the U.S Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was partially supported by DARPA.


text:
# Double-periodic blue variables in the Magellanic Clouds

## 1 Introduction

Over the past years, the microlensing projects (OGLE, MACHO, EROS) have monitored millions of stars in the Magellanic Clouds and Galactic bulge for variability. The resulting huge photometric databases are very well suited not only for microlensing studies but also for many other issues of modern astrophysics, including the distance scale, variable stars, star clusters etc. In particular, the OGLE-II project (Udalski, Kubiak and Szymanski, 1997), has provided accurate BVI measurements for about 6.5 million stars from the central parts of the Magellanic Clouds (Udalski et al. 1998, 2000). Based on this same material, a unique catalog containing about 68.000 variable stars has just been released (Żebruń et al. 2001).

One example of the scientific information that it is possible to extract from such a catalogue is shown by the work of Mennickent et al. (2002, hereafter M02), who report the discovery of Be star candidates in the Small Magellanic Cloud showing unexpected photometric variations. Basically, these authors found four types of variability in targets within the luminosity-colour box of typical Be stars: Type-1 stars showing outbursts, Type-2 stars showing sudden magnitude jumps, Type-3 stars showing periodic variations and Type-4 stars showing random variability. According to these authors, possible causes for Type-1 and Type-2 stars include blue pre-main sequence stars surrounded by thermally unstable accretion discs and white dwarfs accreting from the Be star envelope in a Be+WD binary. Type-4 stars could be classical Be stars. On the other hand, Type-3 stars could not be linked to the Be star phenomenon at all, mainly due to their rather red colours and strictly periodic behaviour. The above study is complemented by the work by Keller et al. (2002) on blue variable stars from the MACHO database in the LMC. These authors report basically the same photometric behaviour in a sample of LMC blue variables. These authors also find emission lines in the spectra of many of these variables.

Further inspection of the misterious Type-3 stars found in the SMC and the LMC (Mennickent et al. in preparation) revealed 30 objects showing double-periodic photometric variations: a long-period sinusoidal variation and a short-period modulation. Some of these short-term light curves are typical of eclipsing Algol variables. The study of this new double-period photometric variability is the subject of this paper.

## 2 Light curve analysis and the long-term oscillation

Table 1 lists the Type-3 stars in the SMC and LMC where two photometric periods were found while Fig. 1 shows some examples of light curves. MACHO light curves were also inspected to confirm the long-term variations and search for colour variability. MACHO identifications are also provided along with the OGLE-II names in Table 1. In general, the long-term periodicity was evident from the first inspection of the light curve. The period $`P_2`$ in these cases was calculated using algorithms based on the $`F`$ (variance ratio) statistics and fourier transform giving consistent results. The results shown in Table 1 reveal that the long-term variations have periods between 140$``$960 d. The period error was calculated as the half width at half maximum of the periodogram peak. Trial and error tests showed that these are upper limits for the error and the accuracy of $`P_2`$ in Table 1 is likely better than shown by a factor 2 at least. After finding $`P_2`$, we applied a non-linear least squares fit to the data of the type:

$`I=I_0A_Isin(2\pi (HJDHJD_0)/P_2)`$ (1)

where $`I_0,A_I`$ and $`HJD_0`$ were parameters to be determined while the period $`P_2`$ remained fixed. The half-amplitudes $`A_I`$ found were between 0.03 and 0.26 mag. The residuals to the best fit were examined with the same period searching routines employed for the analysis of the long-term oscillation. The surprising result was that in all cases listed in Table 1 we find a second period $`P_1`$, in the range of 2$``$16 days, which is listed along with their error in Table 1. For non-eclipsing stars this error was found in the same way that for $`P_2`$, but for 5 eclipsing stars the accuracy was usually much better and the error was estimated interactively.

We also studied the MACHO data for all the variables in our sample, constructing the $`b`$, $`r`$ and $`br`$ light curves, and searched for periodicities. We confirmed the periods found on the basis of the OGLE-II database, and found that for all cases the long-term oscillation is redder at maximum (Fig. 2). In general, the colour variations are of very low amplitude (Table 1).

## 3 On the short-term light curves

Once the long-period sinusoidal modulation is subtracted, the phase curve of the residuals (Fig. 1) turned out to be of three types: eclipsing, double-wave and single-wave. The light curve of the eclipsing stars are typical of those found in Algol type variables. On the other hand, the fact that the stars showing double-wave or single-wave light curve share similar long-term variability with the eclipsing stars and exhibit similar correlation between the short-term and long-term period (see next section) suggests that we are in presence of an homogenous group of short-period variables, and that these non-eclipsing stars could also be Algol-like stars but seen under higher orbital inclinations. The period range along with the fact that emission lines have been observed in Type-3 stars (M02 and Keller et al. 2002) are consistent with this interpretation. In this view, double-wave light curves arises from the changing aspects of a non-spherical Roche-lobe filling secondary star, as occurs in the suspected Algol-type binary V 1080 (Šimon et al. 2000). It is worth to mention that long-term periodic oscillations have been repported in a few Algol systems but with very low amplitude (e.g. Walter 1981). The fact that Algols contain B-A type primaries and cooler secondaries should explain the distribution of Type-3 stars in their colour-colour diagram (M02). In 5 cases we found single-wave light curves. It is possible that the true period is twice the reported period in these cases. This is confirmed by the similar power shown by $`P_1`$ and 2$`\times P_1`$ in the periodogram and also by the good fit of $`2\times P_1`$ by the $`P_1P_2`$ linear correlation (see next section).

## 4 The $`P_2P_1`$ relationship

In Fig. 3 we show what we consider the most surprising result of this paper. The $`P_2`$ and $`P_1`$ periods seems to be correlated. All five deviant points correspond to single-wave light curves, so it is possible that the true orbital period has been missed by a factor 2. The value $`2\times P_1`$ is also indicated in the figure. A linear least squares fit passing by zero gives:

$`P_2=35.17(75)P_1`$ (2)

with a correlation coefficient 0.97 and standard deviation of 28 days. The existence of such a relationship suggests that the phenomenon causing the long-term oscillation could be directly related to the binary nature of the system.

## 5 Discussion: on the nature of the long-term oscillation

The long-term well behaved sinusoidal modulations observed in our sample is not known to occur in other types of blue variable stars. Here we propose as a possible cause for the long-term oscillation the precession of an elliptical disc around a blue star in a semi-detached binary system. The origin of the disc ellipticity and precession could be the tidal interaction between the disc and the Roche-lobe filling secondary star, in a similar way that occurs in short-orbital period dwarf novae of the SU UMa class during superoutburst (e.g. Patterson 2001).

The 3:1 resonance between a disc particle orbiting the primary and the binary system occurs at $`R_{disc}0.46a`$, where $`a`$ is the binary separation. Elliptical orbits at this radius will experience a dynamical apsidal advance with a period given by:

$`P_p=\frac{\sqrt{1+q}}{0.37q}(R_{disc}/0.46a)^{2.3}P_o`$ (3)

(Murray 2000) where $`q`$ is the ratio between the mass of the secondary star and the primary star and $`P_o`$ the orbital period. If we interpret $`P_1`$ as the orbital period and $`P_2`$ as the precession period, then the above equation imposes a strong correlation between both periodicities, as effectively seen in Fig. 3. The observed slope of $`P_2/P_1`$ $``$ 35 should imply a mass ratio $`q`$ 0.08 for $`R_{disc}=0.46a`$. There are few Algols with $`q<`$ 0.1 (e.g. Richards & Albright 1999), and they do not show the long-term oscillations described in this paper. However it is exciting that only in these low $`q`$ systems two important conditions are fulfilled: the primary radius is smaller than the stream circularization radius, being possible the formation of an accretion disc (e.g. Richards & Albright 1999) and the 3:1 resonance radius is below the tidal truncation radius, so the disc could in principle grown beyond the 3:1 resonance and may eventually start to precess (Fig. 4). Long-term spectroscopic observations and detailed modeling are required to confirm this scenario.

###### Acknowledgements.

REM acknowledges support by Grant Fondecyt 1000324 and DI UdeC 202.011.030-1.0. MPD thanks CNPq support under grant \# 301029. We thank the OGLE collaboration team for making their data public domain. This paper also utilizes public domain data originally obtained by the MACHO Project.

Figure Captions

Fig. 1: Eclipsing (upper graph) and double-wave (below graph) light curves. The three pannels in every graph show, from bottom to top, the long-term light variability, the residuals from a sinusoidal fitting and the residuals folded using the short-term period.

Fig. 2: MACHO $`b`$ and $`br`$ light curves folded with the long-term period.

Fig. 3: The long-term period $`P_2`$ versus the short-term period $`P_1`$. The errors in $`P_1`$ are smaller than the used symbols. Dashed horizontal lines connect two possible solutions for stars with single-wave short-term light curve. The linear fit given by Eq.(2) is also shown.

Fig. 4: The primary radius in units of the the binary separation for the Algols shown by Richards & Albright 1999 (filled circles). The circularization radius, the tidal truncation radius and the 3:1 resonance radius, all of them in units of the binary separation, are shown as a function of the mass ratio. Theoretical expressions for these quantities have been taken from Warner (1995). It is possible that only Algols with very low mass ratios could maintain (precessing) accretion discs beyond the 3:1 resonance.