Dynamics of thick, open spirals in Perlas potentials

TL;DR

This study explores the three-dimensional stellar dynamics supporting open spiral arms in galaxies using the PERLAS potential, identifying orbit families that sustain spiral structures and their stability across different mass ratios.

Contribution

It extends previous two-dimensional studies by analyzing 3D orbit families and their stability in spiral galaxy models with varying spiral-to-disk mass ratios.

Findings

Optimal support for 25° pitch angle spirals occurs at M_s/M_d between 0.03 and 0.07.

Strong spiral patterns are supported between the 2:1 and 4:1 resonances.

A bar-like structure appears inside the 2:1 resonance in all models.

Abstract

The PERLAS potential has been successfully used in many studies related with the dynamics of the spiral arms \textit{on} the equatorial plane of normal (non-barred) spiral galaxies. In the present work we extend these studies by investigating the three-dimensional dynamics of the spiral arms in the same type of potential. We consider a typical open, logarithmic, spiral pattern of pitch angle 25$^{\circ}$ and we examine the stellar orbits that can support it as the ratio of the masses of the spiral over the disk component ($M_{s}/M_{d}$) varies. We indicate the families of `three-dimensional' periodic orbits that act as the backbone of the spiral structure and we discuss their stability in the models we present. We study further the quasi-periodic and non-periodic orbits in general that follow spiral-supporting orbits as the $M_{s}/M_{d}$ ratio increases. We find that a bisymmetric spiral with 25$^{\circ}$ pitch angle is better supported by orbits in models with $0.03\lessapprox M_{s}/M_{d} \lessapprox 0.07$. In these cases a strong spiral pattern is supported between the radial 2:1 and 4:1 resonances, while local enhancements of the imposed spirals are encountered in some models between 4:1 and corotation. A characteristic bar-like structure is observed in all models at radii smaller than the radius of the 2:1 resonance.

Tables

Table 1. Table 1: Parameters of the galactic models.

Galaxy type	Spiral Locus	Arms Number	Pitch Angle $i^o$	$\mu =M_{sp}/M_D$	Scale length (kpc)	${\Omega }_p$ ($\text{km}/\text{s}/\text{kpc}$)
Sc	logarithmic	$2$	${25}^o$	$0.01-0.10$	$3.7$	$-20$
Axisymmetric Components
$M_D/M_H$	$M_c/M_D$	Rotation Velocity222Maxima of the rotation velocity. (${\text{km}\text{s}}^{-1}$)	$M_D\left({10}^{10}{M}_{\odot }\right)$	$M_c\left({10}^{10}{M}_{\odot }\right)$	$M_H\left({10}^{11}{M}_{\odot }\right)$	Disk Scale-length (kpc)
$0.1$	$0.2$	$170$	$5.10$	$1.02$	$4.85$	$3.7$
	$b_c$ (kpc)	$A$ (kpc)	$B$ (kpc)	$a_H$ (kpc)
	1	5.32	0.25	12

Table 2. Table 2: Forcing at the main resonances for models M1, M4, M7 and M10.

model	$M_s/M_d$	ILR	4:1	corotation
		(2.03 kpc)	(5.23 kpc)	(8.63 kpc)
M1	0.01	0.024	0.023	0.016
M4	0.04	0.097	0.094	0.064
M7	0.07	0.170	0.165	0.112
M10	0.10	0.243	0.236	0.175

Equations

$${\Phi }_{cm}(x,y,z)=-\frac{G{M}_c}{\sqrt{x^2+y^2+z^2+b_c^2}},$$

$${\Phi }_D(x,y,z)=-\frac{G{M}_D}{\sqrt{x^2+y^2+(A+\sqrt{z^2+B^2}{)}^2}},$$

$${\Phi }_H(r)=-\left(\frac{M(r)}{r}\right)-\left(\frac{M_H}{1.02a_H}\right)\times$$

$${\left[-\frac{1.02}{1+(r/a_H{)}^{1.02}}+\ln (1+(r/a_H{)}^{1.02})\right]}_r^{100},$$

$$M(r)=\frac{M_H(r/a_H{)}^{2.02}}{1+(r/a_H{)}^{1.02}}.$$

$${\nu }^2(R)={\left(\frac{{\partial }^2{\Phi }_0}{\partial z^2}\right)}_{z=0},$$

$$\frac{\nu }{\Omega -{\Omega }_p}=\frac{n}{m}.$$

$$H=\frac{1}{2}\left(p_x^2+p_y^2+p_z^2\right)+\Phi (x,y,z)-{\Omega }_p(x{p}_y-y{p}_x)=E_j$$

$$\dot{x}=p_x+{\Omega }_{b}y$$

$$\dot{y}=p_y-{\Omega }_{b}x$$

$$\dot{z}=p_z$$

$${\dot{p}}_x=-\frac{\partial \Phi }{\partial x}+{\Omega }_b{p}_y$$

$${\dot{p}}_y=-\frac{\partial \Phi }{\partial y}-{\Omega }_b{p}_x$$

$${\dot{p}}_z=-\frac{\partial \Phi }{\partial z}$$

$$\mathbf{\xi }=M{\mathbf{\xi }}_0,$$

$${\lambda }^4+\alpha {\lambda }^3+\beta {\lambda }^2+\alpha \lambda +1=0$$

$${\lambda }_i,1/{\lambda }_i=\frac{1}{2}[-b_i\pm (b_i^2-4{)}^{1/2}]$$

Full text

chaos—galaxy: kinematics and dynamics —galaxy: structure

††journal: ApJ

1 Introduction

The PERLAS potential has been extensively used in many works that study the dynamical properties of the spiral arms on the equatorial plane of disk galaxies (Pichardo et al., 2003; Martos et al., 2004a, b; Allen et al., 2008; Bellini et al., 2010; Antoja et al., 2011; Pichardo et al., 2012; Pérez-Villegas et al., 2012, 2013; Moreno et al., 2014; Martinez-Medina et al., 2015; Moreno et al., 2015; Pérez-Villegas et al., 2015; Martinez-Medina et al., 2016a, b). All these works shed light to the basic orbital dynamics in normal (non-barred) spiral galaxy models. Despite the fact that the regions of galactic disks where the spirals exist are thin, rendering the two-dimensional (2D) modeling as reliable, galaxies are intrinsically three-dimensional (3D) objects. Thus, additional dynamics introduced by orbital instabilities due to vertical perturbations, the role of the vertical resonances present in the disk, as well as the regular or chaotic character of the orbits that can be used for building 3D spiral arms, should be investigated taking into account the existence of the third dimension.

The PERLAS potential is suitable for such a study. It is 3D by construction and its density remains positive everywhere in the parameter space (Pichardo et al., 2003). This is an important advantage of the present model, compared with other 3D potentials used in the past for the same purpose (Patsis & Grosbol, 1996), since our goal is to trace the differences in the dynamics of models when one or more parameters vary. In PERLAS the spiral arms are introduced in a different way than in the models typically employed in literature. Namely, in this case, the spiral arms potential is not an ad hoc perturbation represented by a simple function (or the addition of several of them). Instead, PERLAS is based on a mass density distribution that shows the intricacies and complications of a more realistic large scale spiral mass distribution.

There are many issues that have to be addressed for constructing a 3D spiral pattern. The basic property of the model is to harbour spiral supporting orbits with the appropriate morphology. Such orbits, when projected on the equatorial plane, should precess with respect to each other as their energy varies, in such a way that their apocentra remain close to the imposed potential minima. This way the density will be enhanced along the spiral arms. Typical figures, frequently used to describe this configuration, can be found in Kalnajs (1973, his figure 3) and in Contopoulos & Grosbol (1986, their figure 9). In this consideration of the spiral structure, the orbits are in agreement with the basic principle of the classic density wave theory (Lin & Shu, 1964). However, when the amplitude of the perturbation increases, non-linear phenomena become important and orbital instabilities, as well as considerable morphological deviations of the periodic orbits from ellipses may appear. According to Contopoulos & Grosbol (1988) and Patsis et al. (1991), large amplitudes are inevitable for modeling strong, open, spiral patterns. In the present paper we want to investigate the limitations that are imposed by the strength of the perturbation for building 3D spiral arms. The goal of our study is to explore the orbital dynamics of normal, open, thick spirals, that can be considered of Sc morphological type. Another constrain for building 3D spiral arms is their thickness. Since we do not see morphological features attributed to the spirals exceeding the disks of spiral galaxies observed edge-on, one has to assume a maximum height when considering orbits. The orbits should not exceed the thickness of the disk at any distance from the center.

Investigating the 3D spiral structure in a case of an open, bisymmetric spiral using PERLAS potentials is ideal for a detailed study of the effects introduced in the system when the strength of the spiral varies. In our work we study orbits that support a spiral pattern starting at the 2:1 resonance, inside corotation. Thus, in the present paper we do not investigate in detail the role of chaotic orbits that start at the unstable Lagrangian points region and extend beyond corotation (see Patsis, 2006; Voglis et al., 2006; Romero-Gómez et al., 2006, and all relevant papers on the subject by these groups that followed). However, a possible collaboration of the two mechanisms should be further investigated, as in some models the presence and synergy of both of them can reproduce morphologies encountered in some grand-design spirals (Patsis & Tsigaridi, 2017).

We briefly present the potential components in Section 2 and the algorithms we use to calculate the periodic orbits and their stability in Section 3. Then, in Section 4 we study successively for models with $M_{s}/M_{d}=0.01,0.04$ and 0.07, the 2D and 3D periodic orbits that are the backbones of each model and we discuss their origin and their stability. The orbital content of the models based on quasi-periodic and non-periodic orbits, as well as response models describing their overall morphology are presented in Section 5. In Section 6 we present PERLAS response models that summarize and verify the orbital analysis. Finally we enumerate our conclusions in Section 7.

2 Potential and parameters used

The potential model includes an axisymmetric component that has three parts. For the central mass distribution (bulge) we adopt the spherically symmetric version of the Miyamoto & Nagai (1975) potential, namely, in usual Cartesian coordinates,

[TABLE]

where $M_{c}$ is the mass of the central bulge, $b_{c}$ is a scale length parameter, and $G$ is the gravitational constant.

For the 3-dimensional disk we use again the Miyamoto & Nagai (1975) model, this time in its general form, i.e.

[TABLE]

where $M_{D}$ is the mass of the disk and $A$, $B$ are scale lengths.

For the massive halo we used the potential of Allen & Santillan (1991), which at radius $r$ is given by

[TABLE]

where

[TABLE]

Here $M(r)$ has mass units, $M_{H}$ is the mass of the halo, and $a_{H}$ is a scale length.

Superposed to the axisymmetric components, a non-axisymmetric component is included that represents the spiral arms. For this component, we used the PERLAS potential (Pichardo et al., 2003). It is a bisymmetric, three-dimensional potential and is shaped by a density distribution formed by individual, inhomogeneous, oblate Schmidt spheroids (Schmidt, 1956), superposed in this study along a logarithmic spiral locus of constant pitch angle $i$. The spirals are considered to unwind clockwise.

The spheroids used to model the spirals have constant semi-axes ratio and their density falls linearly outwards, starting from their centers on the spiral locus. Their total width and height are 2 and 1 kpc, respectively; the separation among the centres of the spheroids is 0.5 kpc. The spiral arms formed this way begin at the ILR and end at a distance $1.5$ times the corotation radius. Their density along the locus falls exponentially, as the one of the disk does. We consistently add the spiral arms mass subtracting it from the disk to maintain the model invariable in mass. For details on the spiral arms model PERLAS, see for instance Pichardo et al. (2003); Antoja et al. (2009); Pérez-Villegas et al. (2012, 2013).

Having in mind to model the open spirals of an Sc type galaxy, we adopted the parameters used for that purpose by Pérez-Villegas et al. (2015). We briefly summarize them in Table 1. We give the parameters for the axisymmetric components, as well as those for the spiral perturbation.

The rotation curve and its decomposition is given in bottom panel of figure 1 in Pérez-Villegas et al. (2015). The two most characteristic parameters of our model are the pattern speed $\Omega_{p}=-20$ km s*-1* kpc*-1* with which our system rotates and the pitch angle of the logarithmic spirals, for which we have adopted the value $i=25^{\circ}$, typical for an Sc galaxy. The amplitude of the spiral perturbation is determined by the ratio $M_{s}/M_{d}$. This is the parameter we have varied in our study in the range $0.01\lesssim M_{s}/M_{d}\lesssim 0.1$.

An informative diagram about the properties of the potential is Fig. 1. In this figure we give the effective potential in the equatorial plane $(x,y)$ in the cases with $M_{s}/M_{d}=0.04$ (a) and $M_{s}/M_{d}=0.1$ (b). The first case corresponds to one of the standard models of this study, which will be discussed in detail in the following sections. The second is the model with the most massive spirals we studied, so the properties of the imposed potential are displayed pronounced.

Estimating the basic radial resonances of the system from the axisymmetric part of the model, corotation is at $8.63$ kpc, the 2:1 resonance at $2.03$ kpc and the 4:1 resonance at $5.35$ kpc. (Thus, in our models the spiral arms end at 12.9 kpc, which is 1.5 times the corotation radius). However, small displacements are expected to occur when the perturbation is introduced, especially when the spiral is strong. This will be discussed for each model separately. In addition, in our model we have also vertical frequencies, in the $z$ direction, defined by

[TABLE]

where $\Phi_{0}$ is the axisymmetric potential, with the corresponding vertical resonances determined by

[TABLE]

In practice the location of the radial and vertical resonances should be always taken from the variation of the stability curves of the spiral-supporting families of periodic orbits, since these curves have specific properties at the resonances (see Sections  3 and 4 below).

3 Periodic orbits and stability indices

In order to study the orbital behavior of our galactic models, we use the Hamiltonian formalism. As in the previous, 2D, PERLAS studies, the galaxy is modelled as an autonomous Hamiltonian system, rotating with pattern speed $\Omega_{p}$, in which case the Hamiltonian can be written as

[TABLE]

where $x,y,z$ are Cartesian coordinates in the rotating frame, $p_{x},p_{y},p_{z}$ are the canonically conjugate momenta, $\Phi(x,y,z)$ is the total potential (including axisymmetric and non-axisymmetric terms) and $E_{j}$ is the Jacobi’s constant (in the text we may also refer to it as the “energy”).

The equations of motion are

[TABLE]

In order to calculate orbits we start in the plane $y=0$ with $p_{y}>0$. Due to the existence of the $E_{j}$ integral the needed initial conditions can be reduced in four $(x_{0},p_{x_{0}},z_{0},p_{z_{0}}$). Our system rotates clockwise, thus the value of the initial $x_{0}$ coordinate will be negative. For finding periodic orbits we use an iterative Newton method with an accuracy of at least $10^{-10}$. The equations of motion for this purpose are integrated using a Runge-Kutta integrator of order four. For the integration of chaotic orbits and in the calculations of Poincaré surfaces of section we have used also a Runge-Kutta Fehlberg 7-8th order, scheme (Fehlberg, 1968).

The computation of the stability indices of the periodic orbits is based on the theory of variational equations. We consider small deviations from its initial conditions, and then the orbit is integrated. Therefore the initial and final deviation vectors are related as:

[TABLE]

where $M$ is the monodromy matrix. The characteristic equation of this matrix has the form

[TABLE]

Its solutions obey the relations $\lambda_{1}\lambda_{2}=1$, $\lambda_{3}\lambda_{4}=1$ and for each pair we can write

[TABLE]

where $b_{i}=\frac{1}{2}(\alpha\pm\Delta^{1/2})$ and $\Delta=\alpha^{2}-4(\beta-2)$. $b_{1}$ and $b_{2}$ are the stability indices. If $\Delta>0$, $|b_{1}|<2$ and $|b_{2}|<2$, the four eigenvalues are on the unit circle and the orbit is stable. If $\Delta>0$, $|b_{1}|>2$ and $|b_{2}|<2$ or $|b_{1}|<2$ and $|b_{2}|>2$, two eigenvalues are on the real axis and two on the unit circle, and the orbit is simple unstable. If $\Delta>0$, $|b_{1}|>2$ and $|b_{2}|>2$, all four eigenvalues are on the real axis and the orbit is called double unstable. Finally if $\Delta<0$, all four eigenvalues are complex numbers but off the unit circle and then the orbit is characterized as complex unstable. For a detailed description of the method the reader is referred to Broucke (1969) and Hadjidemetriou (1975).

4 2D and 3D periodic orbits as building blocks

Sketches like the one by Kalnajs (1973, his figure 3) illustrate the assumed stellar (and to a large degree also gaseous) flows in a barred and in a spiral case. In 2D models the drawn ellipses correspond to elliptical periodic orbits (hereafter p.o.) that belong to the x1 family (Contopoulos & Grosbol, 1989). In cases in which the elliptical p.o. have their apocentra along an axis, they support a bar, while in cases they have their apocentra close to spiral loci, they support a spiral pattern. Essentially, in both cases, it is the same family of p.o. found in different potentials. A bisymmetric spiral pattern can be seamlessly reinforced in such a configuration between the radial 2:1 and 4:1 resonances. Beyond that resonance the ellipses obtain a rectangular character that leads either to boxy bars (Contopoulos, 1980; Patsis et al., 1997a), or to bifurcations of the arms and/or an overall deterrence of the evolution of open, bisymmetric stellar spirals in normal spiral potentials (Contopoulos & Grosbol, 1986; Patsis et al., 1991).

In addition, in a 3D case, the x1 family is substituted by a tree of families of p.o., called by Skokos et al. (2002) the “x1-tree”. The x1-tree includes, besides the planar x1 family on the equatorial plane, also the vertical bifurcations of x1, which are introduced in the system in pairs at the vertical resonances, starting from the vertical 2:1 resonance. In Skokos et al. (2002) they have been named x1v1, x1v2 (the two families that are associated with the vertical 2:1 resonance), x1v3, x1v4 (the two families that are associated with the vertical 3:1 resonance), etc. This is the usual succession of the 3D families of p.o. in rotating galactic models. Exceptions have been found in particular cases (Patsis & Harsoula, 2018), but usually the families introduced as stable in the system are those that have as last digit in their name in the Skokos et al. nomenclature an odd number (x1v1, x1v3, etc.). Looking for 3D stable families to support the spiral arms, these are the obvious first candidates to be examined as possible building blocks of a thick spiral. These families, being bifurcations of the planar x1 family at a transition of stability from stable to simple unstable at a vertical resonance, have by definition exactly at the bifurcation point a morphology identical to x1. However, as soon as we depart from the bifurcating point towards larger energies (Jacobi constants) they develop a typical morphology for each family. As in the case of barred models (Skokos et al., 2002) the p.o. of the x1-tree are expected to have at a certain energy range elliptical projections on the equatorial plane, similar to the x1 ellipses. Then, gaining in height, they will have projections that will not be able to be combined with the x1 p.o. and support the same spirals. We have to underline though that the typical edge-on morphologies associated with each family in barred potentials is because of their orientation with respect to the major and minor axes of the bar, which coincide with the corresponding axes of the elliptical orbits. In a spiral potential however, since the axes of the 3D periodic orbits of a family projected on the equatorial plane precess as their $E_{j}$ varies, their projections on the $(y,z)$ and $(x,z)$ planes are not expected to be in general the known, recognizable shapes of the corresponding side-on and end-on views of the barred cases.

We note that in a non-axisymmetric system the most reliable way of locating the resonances, is by means of its periodic orbits. At the radial and vertical resonances new families are introduced. Thus, one has to find the bifurcating points along the characteristic curve, or along the stability diagram, of the main family x1 and associate them with a $n:1$ resonance. The location of a resonance specified this way, is expected to be different from a direct estimation from the $\Omega\pm\lambda/n$ vs $r$ curve (where $\lambda$ is either the epicyclic or the vertical frequency of the radial and vertical resonances respectively), which is based on the axisymmetric part of the potential. In weak spiral models this difference can be practically ignored, but it has to be taken into account in strong spiral cases.

In Table 2 we give for the models we will refer to in the following sections the forcing, i.e. the maximum tangential over the axisymmetric force, at the resonances.

4.1 Model $M_{s}/M_{d}=0.01$

The less massive spiral we have examined has $M_{s}/M_{d}=0.01$. We will refer to it as Model “M1”. Fig. 2 gives the $(E_{j},x_{0})$

(for a definition see Contopoulos & Grosbol, 1989) for the x1-tree families of this model. It is composed by the individual projections of the characteristic curves, while in 3D systems the p.o. demand in general four initial conditions $(x_{0},z_{0},p_{x_{0}},p_{z_{0}})$ in order to be uniquely specified and be plotted against $E_{j}$. Even for the characteristic curve of the planar x1 family we need two initial conditions $(x_{0},p_{x_{0}})$, because in general the x1 ellipses in a spiral model have $p_{x_{0}}\neq 0$. Nevertheless the $(E_{j},x_{0})$ projection gives information about the Jacobi constant in which a family is introduced in the system and its extent, so we will use it, together with the stability diagrams, for the description of the models.

In Fig. 2 the black curve is the characteristic of the x1 family, the red curve that of the 3D x1v1 family, the green curve is the x1v3 characteristic, the blue one corresponds to the x1v5, while the curves depicted in light blue are 2D orbits close to corotation, bifurcated at radial $n:1$ resonances with $n\geq 4$. $E_{j}$ is given in units of $(10km/s)^{2}$. The same holds for all figures $E_{j}$ appears, throughout the paper.

We restrict ourselves to the study of the 3D families bifurcated from x1 as stable, since they are expected to attract around them regular orbits that will support the thick spiral arms. The way these families are introduced in the system is usually presented with a stability diagram, which gives the variation of the stability of the p.o. by means of the two stability indices, $b_{1}$ (for the radial perturbations) and $b_{2}$ (for the vertical perturbations). In Fig. 3 we show the stability diagram for model M1. We focus ourselves at the intersections of the x1 vertical stability index $b_{2}$ with the $-2$ axis. At these points new 3D families, of the same multiplicity with x1, are introduced in the system (Contopoulos & Magnenat, 1985). Usually the $b_{2}$ index dives below the b$=-2$ axis for a certain $\Delta E_{j}$ interval and then it returns back to values larger than $-2$, giving rise to successive $S\rightarrow U\rightarrow S$ transitions. The x1v1, x1v3, x1v5 families are bifurcated from x1 at the $S\rightarrow U$ stability transitions, as stable.

In order to have a better view of the variation of the stability indices, Fig. 3 is split in two parts. In Fig. 3a we give the stability curves for $-1600\lesssim E_{j}\lesssim-1186.45$ and in Fig. 3b for $-1186.45\lesssim E_{j}\lesssim-1022.71$. The indices of x1 are drawn with black, while its 3D bifurcations are drawn with the colours used for the same families in Fig. 2. They are also indicated with arrows. In the $E_{j}$ intervals in which a family becomes complex unstable we draw a straight line segment at $b_{1,2}=0$. The vertical red line in figures 2 and 3 denotes the location of the 4:1 resonance in model M1.

We used these p.o. in order to find a reasonable orbital backbone that could support our imposed 25*∘* bisymmetric spiral. However, as we can observe in Fig. 4, such a backbone is not provided by model M1. In Fig. 4 the p.o. belonging to the x1-tree (x1, x1v3, x1v5) are plotted in black. Firstly we observe that inside the radial 2:1, Inner Lindblad Resonance (ILR), denoted with a magenta circle, the members of the x1-tree p.o., belonging to x1 and x1v1 families, are elongated red ellipses. Despite the fact that the arms of the spiral potential end at ILR, at $r=2.03$ kpc, the particles inside this radius feel an m=2 component, due to the presence of the arms. These ellipses have their apocentra almost along an axis. Thus, they could support, a non-imposed, central bar-like component in the central parts of the model.

Beyond ILR, we find orbits with elliptical projections on the equatorial plane up to the radial 4:1 resonance. However, their apocentra are not aligned with the imposed spiral. They tend to form a tighter and weaker one. The 2D families beyond 4:1, depicted in light blue in Fig. 2, have polygonal shapes with apocentra not aligned with the imposed spirals. They are plotted again in light blue in Fig. 4. We find also 3D p.o, which have a polygonal shape in their projections on the equatorial plane. These projections appear at smaller radii than the 2D orbits for the same energy, since they have also a certain thickness. Such orbits, belonging to the x1v3 family are plotted with green colour in Fig. 4. It is evident that the overall response of the model does not provide a skeleton of stable p.o. with which we could proceed in building a spiral with regular orbits at any region and between any resonances of the model.

4.2 Model $M_{s}/M_{d}=0.04$

The next model we present has $M_{s}/M_{d}=0.04$ and we will refer to it as model M4. In Fig. 5 we give for M4 the characteristic curves of the families of the x1-tree, projected in the $(E_{j},x_{0})$ plane. Colors are as in Fig. 2 and arrows point to each family. For this model we give also the loop of the x2, x3 characteristic, while one more 3D family, x1v7, is included, drawn in magenta color. In Fig. 6 we present the stability diagram for the x1-tree families of this model. We keep the same colors for the families as in Fig. 5. In Fig. 6b, the $E_{j}$ for which the x1 orbits start having a rectangular character is indicated with a dashed vertical line as “4:1”. We indicate this $E_{j}$ also as “4:1”, since it is the shape of the orbits that acts as obstacle in the continuation of the support of a spiral pattern along the loci of the spiral potential by the periodic orbits of the system. This way we have an “effective” 4:1 resonance at a smaller energy than the one predicted for the 4:1 location from the axisymmetric part of the potential.

In Fig. 7 we isolated x1-tree p.o. that offer a backbone for supporting the $25^{\circ}$ spiral of our model. It includes x1, as well as projections on the equatorial plane of p.o. belonging to the families x1v3, x1v5 and x1v7. In the case of x1, the last spiral supporting orbit was at energy $E_{j}\approx-1100$. Beyond that energy the morphology of the x1 orbits was rectangular-like. This is the $E_{j}$ value indicated with an arrow in Fig. 6. The circles drawn in Fig. 7 show the locations of the ILR (innermost), the location of the 4:1 region as defined by the appearance of rectangular-like x1 p.o., the 4:1 radius as expected from the axisymmetric part of the potential and corotation (outermost). We observe that the included p.o. of the x1-tree offer a backbone that tends to support the spiral (plotted with red solid line) up to the 4:1 resonance. The x1 p.o. close to 4:1, but at a slightly higher $E_{j}$ do not help the spiral extending further out (Fig. 8). Contrarily, they tend to support density maxima off the spiral (red curve) that can be described as bifurcations of spiral arms with a pitch angle smaller than $25^{\circ}$, as well as a boxy feature in the middle of the way to corotation. All these are conspicuous in Fig. 8 that shows the orbital backbone offered by x1 p.o. for $E_{j}$’s larger than $E_{j}\approx-1100$.

The projections of the 3D families in Fig. 7 are also subject to a morphological evolution towards rhomboidal shapes as their $E_{j}$’s increase. However, simultaneously, they gain in height, causing their rhomboidal projections to appear at smaller radii in general, than the x1 p.o. Taking into account the morphological evolution of all orbits involved in the building of an orbital skeleton of p.o. that supports a thick spiral, we find that their thickness at the 4:1 resonance region ($r\approx$ 5 kpc) is about 0.3 kpc. The corresponding isodensity contour of the axisymmetric part of the model will be given in subsequent figures with orbital edge-on profiles. This is in agreement with a reasonable geometry for an Sc galaxy. The last 3D spiral enhancing p.o. are found at $E_{j}\approx-1179.523$ for the x1v3 family, at $E_{j}\approx-1143.313$ for x1v5, and at $E_{j}\approx-1100$ for x1v7 (close to the last x1 spiral-enhancing p.o.). We will return to this point in Section 6.

If we consider all p.o. discussed so far and plot them in a single figure we end up with Fig. 9. The p.o. that provide a skeleton for supporting the spiral pattern are plotted in black, while the non-supporting in light blue. It is clear that on the equatorial plane of the model there are p.o. to support the spiral pattern up to the 4:1 resonance region. It is not obvious at all that the light-blue periodic-orbits could provide a backbone that would help the pattern continue towards corotation.

The inner limit of the spiral structure is by construction at the ILR. Appart from the x1 orbits, in all previous figures the 3D families we discussed are the x1v3, x1v5 and x1v7, i.e. stable orbits bifurcated from x1 at the vertical 3:1, 4:1 and 5:1 respectively. The family x1v1, which seems to be very important for the edge-on profiles of barred galaxies (Patsis et al., 2002) has elongated elliptical projections, similar to the x1 p.o. inside the ILR. In the PERLAS models its orbits do not contribute to the reinforcement of the spiral. They support a bar-like structure in the central parts despite the fact that PERLAS does not include an explicit central bar component. This bar-like response of our non-barred, normal spiral model was found in all models we examined. Thus the PERLAS spiral keeps its 25*∘* pitch angle only beyond the inner 2:1 resonance. Crossing this resonance towards the center of the galaxy the pitch angle opens and the response has been always bar-like. Interesting is also to observe the $(y,z)$ edge-on view of a set of x1v1 orbits. We realize that the shape of the orbits supports a peanut-shaped profile (Fig. 10). The y-axis is not aligned with the major axis of the bar. The $\infty$-shape figure observed in Fig. 10 appears in this particular projection.

Models with $0.03\lessapprox M_{s}/M_{d}\lessapprox 0.05$ give similar results with those discribed for model M4. We only observe that the x1 p.o. become rectangular or rhomboidal-like at an earlier energy, $E_{j}$, which deviates more from the location of the axisymmetric 4:1 resonance.

4.3 Model $M_{s}/M_{d}=0.07$

The orbital behavior of our models started changing when we reach the ratio $M_{s}/M_{d}=0.07$, which is our model M7. As we can observe in Fig. 11 that shows the evolution of the stability indices, the vertical resonances appear now shifted to smaller energies than in the previous models. Families x1v1, x1v3 and x1v5, introduced as stable in the vertical 2:1, 3:1 and 4:1 resonances respectively, have larger complex unstable parts, compared to models with smaller $M_{s}/M_{d}$ ratios.

The backbone of p.o. for model M7 is given in Fig. 12. We plot again 2D orbits of the planar families and representatives of the stable 3D families. Apart from the fact that the x1 orbits become rhomboidal at a shorter distance from the center than before, we observe that the candidate periodic orbits to support the spirals have their apocentra ahead, in the sense of rotation, of the potential minima (indicated with a continuous red line). Nevertheless, they are not far away from them.

5 Contribution of quasi- and non-periodic orbits

The comparison of the backbones of p.o. we calculated for models with $0.01\leq M_{s}/M_{d}\leq 0.1$ led to the conclusion that the best choice for building a $25^{\circ}$ spiral is in the range $0.03\leq M_{s}/M_{d}\leq 0.07$. As a typical case we have chosen model M4 to present the contribution of quasi- and non-periodic orbits to the modeled spiral pattern.

5.1 Orbits in model M4

For this purpose we chose a typical $E_{j}$ that harbors a stable x1 p.o. from the backbone of the M4 model spiral. The $(x,p_{x})$ Poincaré surface of section at $E_{j}=-1286.277$ for planar orbits in the case of model M4 is presented in Fig. 13. The location of the x1 p.o. on the surface of section is indicated with the arrow labeled with “a” and its morphology is given in Fig. 14a. The location of the rest of the orbits depicted in Fig. 14 is indicated with arrows in the surface of section in Fig. 13. The letter that characterizes an orbit on the surface of section corresponds to the label of each panel in Fig. 14. Their initial conditions deviate from those of x1 by $\Delta x$ = 0.14, 0.4, 0.5, 0.7, 0.84 and 0.98 kpc for orbits in (b), (c), (d), (e), (f) and (g) respectively. All orbits have been integrated for five x1 orbital periods. As we can observe in Fig. 13, orbits (b), (c) and (d) belong to invariant curves around x1, orbit (e) to an orbit on the chain of islands that follows, while orbits (f) and (g) are chaotic orbits. Orbit (f) is chosen to be at the edge of the stability island and orbit (g) well inside the chaotic sea. By comparing the size, the morphology and the orientation of these orbits with respect to the loci of the 25*∘* spiral, we realize that those in Fig. 14b to e, support during their integration the enhancement of the density of the model at the spiral arms region. Orbit (f) supports the spiral only partly, while orbit (g) does not.

As a next test for estimating the contribution of the 2D quasi- and non-periodic orbits to the enhancement of the spiral structure we applied the same $\Delta x$ perturbations successively to the set of x1 p.o. used in Fig. 7. These p.o, each at a different $E_{j}$, are typical orbital building blocks of the spirals. They have the right orientation and shape, so that they enhance the spirals, due to the relative position of their apocentra with respect to the potential minima of the model. The result is presented in Fig. 15. The x1 p.o. in Fig. 7 have been perturbed by $\Delta x=$ 0.14, 0.4, 0.5, 0.7, 0.84 and 0.98 kpc and three of the resulting sets of non-periodic orbits are given in Fig. 15a to 15c. Namely, we plot the orbits perturbed by $\Delta x=$ 0.14, 0.5 and 0.98 kpc respectively. We realize that the orbits contribute less as we deviate from the initial conditions of the p.o. In all orbits of the specific sample we use to describe this property, as well as in many other cases of orbits we have examined, the quasi-periodic orbits around x1 had the largest contribution to the imposed spiral. The enhancement of the spirals from orbits at the borders of the stability islands on the surfaces of section, or from chaotic orbits, like the one in Fig. 15c, was less evident.

Additional perturbations have been applied to radially perturbed p.o. like those discussed in the previous paragraph, this time in the vertical direction. Typical cases are presented in Fig. 16, which are the orbits depicted in Fig. 15 perturbed also in the $z$ direction. Initially the radially perturbed p.o. orbits have $z_{0}=0$. We started increasing the $z$ coordinate of the orbit and we have checked their projection on the equatorial plane. We stopped increasing $z_{0}$ as long as the projection did not support the spiral any more. The edge-on profiles of these sets of spiral-supporting orbits were similar to those we obtained when we have considered just the corresponding periodic orbits. Namely, the thickness of the profiles did not exceed $z=$0.3 kpc above or below the equatorial plane at $r\approx$5 kpc, while we have $z\approx 0.54$ kpc at $r=0$. The orbits in panels Fig. 16a and Fig. 16b support the spiral. However, all perturbations in the $z$ direction we tried for the orbit in Fig. 15c gave orbits that do not reinforce the spiral pattern. We note that by perturbing a chaotic orbit, this does not necessarily imply that the perturbed orbit will be chaotic as well, since by the displacement of the initial conditions we may reach regions of the phase space occupied by tori.

Along the largest parts of the spiral of the models, besides the stable x1 p.o., exist also stable 3D families of periodic orbits. These orbits will trap again around them regular orbits and it is expected that they will also participate in the enhancement of the density of the spirals. For this purpose we apply perturbations to all spiral supporting p.o., 2D and 3D, depicted in Fig. 7. By doing so, we can find a thick spiral pattern supported by 3D regular orbits up to the radial 4:1 resonance, at the region where the x1 p.o. become rhomboidal. The face-on view of this spiral pattern is presented in Fig. 17a and its edge-on view in Fig. 17b. In Fig. 17b we plot also an isodensity contour of the axisymmetric part of the potential that does not exceed the 0.3 kpc height above or below the equatorial plane at the 4:1 resonance distance. We observe that the total edge-on profile of the set of orbits is included within this contour.

By perturbing the p.o. of the backbones of the models we investigated we find similar results for the cases $0.03<M_{s}/M_{d}<0.07$. The 3D spirals are reinforced by orbits similar to those we described in the previous paragraphs for model M4. In the cases with $M_{s}/M_{d}<0.03$ the ellipticity of the x1 p.o. and that of the projections of its vertical bifurcations on the equatorial plane is small. Considering quasi-periodic orbits around them does not improve the relation between the orbits of the model and the imposed spiral. For $M_{s}/M_{d}>0.07$, apart from the large complex unstable parts of the 3D families the size of the stability islands around x1 decreases, which makes the support of the spiral structure problematic.

6 Response models

Having studied the basic orbital dynamics, we investigate further the overall dynamics of PERLAS models by integrating sets of $4.6\times 10^{4}$ initial conditions. We start by distributing randomly test particles in a disk of 13 kpc radius in order to obtain an initially homogeneously populated disk. The velocities of these particles correspond to those of circular motion in the axisymmetric part of the potential ($M_{s}=0$). Then we increase linearly $M_{s}/M_{d}$ within two system rotations reaching at the end of this period the $M_{s}/M_{d}$ ratio that characterizes each model. Following this procedure, the particles obtain the initial conditions with which they will be integrated in the final full potential. We follow the response of the model to the imposed potential for about 30 system rotations. Besides the 2D responses, we calculated 3D response models as well, by placing randomly the particles in cylindrical configurations, which had a 13 kpc radius and a height of 0.3 kpc above and below the equatorial plane. We have chosen the radius to be the same with the 2D disk, while for the height we have taken into account the total thickness indicated by the orbital models, namely 0.6 kpc. This is a height within which the regular orbits can easily support 3D spirals. The number density of the particles in the initial cylinder is homogeneous. The initial velocities have been chosen to be those for circular motion on the equatorial plane at the same radius. This rather arbitrary choice is used just as a starting point, which leads to the initial conditions in the full potential after the transient two pattern rotations period.

The snapshots of the response models have been converted to images by means of the ESO-MIDAS package, taking into account the local number density of the particles. The results of this exercise pointed again to best matching between the imposed potential and the response model for the cases with $0.03\leq M_{s}/M_{d}\leq 0.07$. If we weight the images taking into consideration that the disk of the imposed potential is exponential with a scale length 3.7 kpc and that the central density of the Schmidt spheroids falls also exponentially with radius with the same scale length as the disk, we end up to figures like Fig. 18, which is our “density map” for model M4. Such images describe the morphology and the relative importance of the morphological features supported by the orbits of the model.

Fig. 18, describes the 3D response of model M4, projected on the equatorial plane. The three drawn circles indicate the ILR (innermost), the 4:1 resonance (middle) and corotation (outermost). The loci of the imposed spirals are drawn with heavy black dots. The overall response is in agreement with the results of past studies on the subject based on orbital theory (Contopoulos & Grosbol, 1986, 1988; Patsis et al., 1991; Patsis & Grosbol, 1996). Namely, we have a strong, bisymmetric spiral pattern that starts existing at the ILR region and it is almost in-phase with respect to the imposed spiral. It broadens and bifurcates just before the radial 4:1 resonance at the points indicated with red arrows. The response morphology of this model is established already after about 10 pattern rotations. The pattern depicted in the model, with the two, main, symmetric spiral arms that are split at larger distances, is frequently encountered in the morphology of non- or weakly-barred spiral galaxies (Grosbol et al., 2004).

Faint extensions of the spirals are observed in the region between the 4:1 resonance and corotation. They are faint, since they are located in low density regions of the disk, due to its exponential profile. In general, extensions of the spiral arms beyond 4:1 are found in the responses of gaseous models. Typical examples are given in Patsis et al. (1997b, see e.g. figures 3 and 4). The difference between those extensions and what we find here, is that now we have mainly fragments of stellar spiral arms in- or nearly in-phase with the imposed spiral. To better describe their properties we present in Fig. 19 the images of the 2D and 3D responses of model M4 without weighting them with the exponential decrease of the disk and spiral densities. Thus, morphological features appear equally important, regardless of their distance from the center of the system. In Fig. 19a it is depicted the response of the 2D, while in Fig. 19b that of the 3D response model (the projection on the equatorial plane). In both panels we observe a gap at the 4:1 resonance region. However, in both cases we find again local density enhancements between 4:1 and corotation, close to the imposed spiral loci of the PERLAS potential. In the 3D case (panel b) the gap is larger (indicated with two arrows along each arm) but the segments that follow are longer than the chunks of local density enhancements we find in the 2D model (panel a). Nevertheless, in neither of the cases we find segments of spiral arms crossing corotation.

The orbits of the particles populating the ring between the resonances 4:1 and corotation have large $E_{j}$’s. Most of them are in the range $-1056<E_{j}<-1036$. A typical $(x,p_{x})$ surface of section for $E_{j}=-1046$ in model M4 is given in Fig. 20a. We observe that the phase space is dominated by a chaotic sea. The stable periodic orbits at this $E_{j}$ do not provide the necessary backbone for supporting the local density enhancements by regular orbits. Thus the segments of spirals observed in the models in this region can be only features reinforced by sticky chaotic orbits. Indeed, the stability islands in the region are tiny. In Fig. 20a we give in enlargement the area of the surface of section that includes the stability island of a “4:1” family with rhomboidal orbits (lower right island). As an example, we give in Fig. 20b quasiperiodic orbits of this family and sticky chaotic orbits from the region around it, that show how the segments could be supported. The detailed mechanism, as well as model dependencies, will be investigated in a subsequent paper. Nevertheless, we note that in this case we have a similar situation like in the barred-spiral model for NGC 4314, studied by Patsis et al. (1997a), where the sticky chaotic orbits around tiny islands affect the morphology of the model at the radial 4:1 resonance region.

At the level of orbital analysis of the present study, it is difficult to name a specific $M_{s}/M_{d}$ ratio that gives the best matching between imposed potential and response model. All 2D and 3D PERLAS response models with $0.03\leq M_{s}/M_{d}\leq 0.07$ have similar overall morphologies. Outside this range the models show conspicuous deviations from this morphology. For example the 3D response of model M1 ($M_{s}/M_{d}=0.01$) does not show a well organized spiral structure. This can be seen in Fig. 21 in which we avoid weighting the image taking into account the exponential profile of the disk. This way we show that there is an intrinsic discrepancy between the orbits of the model and the imposed spiral. In the response we do not find well organized spiral arms. Only parts of them can found along the imposed spirals. Scattered arcs of enhanced local densities, in general not along the loci of the imposed PERLAS spirals, are observed in other regions. This is expected from the calculated backbone of periodic orbits provided by the model (Fig. 4). Nevertheless, the characteristic quadruple gap at the 4:1 resonance is discernible in Fig. 21 along the 4:1 resonance circle.

When the mass ratio increases, the alignment of the response density maxima with the imposed spirals improves and we have the level of agreement we described for model M4. However, for $M_{s}/M_{d}>0.07$, the sizes of the stability islands around the x1 p.o. in the $(x,p_{x})$ surfaces of sections are considerably reduced, even at small energies. We find though sticky zones, which provide to the system sticky-chaotic orbits that behave like regular for several rotational periods and support a spiral pattern with $25^{\circ}$ pitch angle in the response models. However, as time increases, the spiral arms become broader and less well defined. A typical example is given in Fig. 22, which is the 3D response model for $M_{s}/M_{d}=0.1$ (hereafter called model M10). In (a) we present all particles projected on the equatorial plane at time corresponding to about 25 pattern rotations. We have used the GALI2 stability index (Skokos et al., 2007, 2008) to characterize the chaoticity of the orbits and based on this we plot with black the particles on regular and with red the particles on chaotic orbits. We use the GALI2 value of the orbit at the time of the snapshot to characterize it and we plot with black the particles on orbits with $\log_{10}$(GALI2)$>-8$. However, at later times the index of an orbit may well fall to smaller values, in which case the particle will be plotted red. This is reflected in Fig. 22a, where we observe that we have more red points closer to the center of the system, because in this region the system is dynamically more evolved. Gradually the red region expands towards larger radii. The fact that the region between the 4:1 resonance and corotation is not red after 25 patter rotations, indicates that many of the chaotic orbits are sticky, mimicking a regular behavior and keeping their GALI2 index in the range with $\log_{10}$(GALI2)$>-8$. Building a spiral pattern by means of such orbits has as a result mainly to obtain broad spiral arms with fuzzy boundaries in the responses. This can be observed in Fig. 22b, where we give the image of model M10, non-weighted for its exponential character. The response spirals are still formed close to the loci of the imposed spiral, with a tendency to have their main part ahead of them, but their boundaries have an irregular character. Also the characteristic bifurcations of the arms at the 4:1 resonance are much less conspicuous.

Although models with $0.03\lessapprox M_{s}/M_{d}\lessapprox 0.07$ ratios match better the imposed spirals than models with $M_{s}/M_{d}>0.07$, the “strong” models of our sample show that the role of sticky-chaotic orbits can be in some cases important for the dynamics of spiral galaxies.

In Fig. 22 we observe two more features that we would like to discuss. With a spiral forcing $F_{\theta}/F_{r}\approx 0.175$ at the corotation region, M10 is the strongest model we have studied. Nevertheless, this forcing is much smaller than what we encounter in the corotation region of fast, strong bars (Block et al., 2001). The banana-like periodic orbits around the stable Lagrangian points L4 and L5 (Contopoulos & Grosbol, 1989) are stable in large fractions of their characteristic in M10. Thus, the regions around L4, L5 are populated and we observe the formation of banana-like features. These regions do not reach the neighborhood of L1 and L2, because the p.o. of the banana-like family at larger $E_{j}$ are unstable. As a result the regions close to the unstable Lagrangian points are rather empty (unlike with what we observe in Fig. 21 for model M1 and in Fig. 19 for model M4). Despite the fact that banana-like orbits have been associated with the reinforcement of spiral arms in some models (see e.g. the model of Lépine et al, 2017, for the Local Arm), such morphological features are hardly observed in real galaxies, at least as strong features. Two reasons that could lead to the absence of these structures in real galaxies are (i) large corotation distances from the center or (ii) strong perturbations. As we have noted, Fig. 22b is not weighted taking into account the exponential character of the imposed disk and spirals. However, in images like Fig. 18 the corotation region is a low density region and the banana-like features are not discernible. The second option is that the families around L4 and L5 are unstable. This can happen in barred-spiral models, where the forcing at the corotation region is strong. In the latter case, particles located at the L4, L5 regions, are not trapped by regular orbits. Their motion will be determined firstly by the presence of the stable branch of the unstable manifold of the unstable p.o. around L1, L2. This will lead them to the neighborhood of either L1 or L2 and then, the unstable branch of the manifold will lead them beyond corotation, reinforcing “chaotic” spirals. This mechanism for populating or depleting the areas around L4 and L5 in a response model by changing the forcing has been presented in Patsis & Tsigaridi (2017). In Fig. 22b we can observe weak “chaotic” spirals emerging from the L1 and L2 regions. The dynamical mechanisms that act in this case are similar to those in the models in Patsis & Tsigaridi (2017), where we have two sets of spirals supported by regular (inside corotation) and chaotic orbits (beyond corotation) coexisting (see e.g. their figure 2).

7 Conclusions

We have investigated the orbits that can support an open, thick (3D) spiral structure in the PERLAS potential. We have taken advantage of the PERLAS property to have everywhere positive density as the parameters of the spiral potential part vary and we investigated models in the range $0.01\leq M_{s}/M_{d}\leq 0.1$. Our main conclusions are the following:

A thick spiral pattern in the PERLAS potential is supported by orbits associated with the stable vertical bifurcations of x1 at the vertical $n:1$ resonances with $n\geqq 3$. The thickness of the spirals supported by such orbits in our models is about 0.6 kpc (i.e. 0.3 kpc above and below the equatorial plane). Fig. 23 illustrates how the periodic orbits in our model shape an appropriate backbone for reinforcing a thick spiral. Around these p.o. can be build a 3D spiral.

By varying the mass ratio $M_{s}/M_{d}$ in our models we concluded that the best relation between the imposed spiral and the orbital content of the models is for $0.03\lessapprox M_{s}/M_{d}\lessapprox 0.07$. 3. 3.

We find a thick, strong, bi-symmetric spiral pattern extending between the ILR and the radial 4:1 resonance of the model in agreement with previous studies by Contopoulos & Grosbol (1986, 1988), Patsis et al. (1991), Patsis et al. (1994), Patsis & Grosbol (1996) and Patsis et al. (1997b). The characteristic morphological feature associated with the 4:1 resonance in all models, was the bifurcation of the arms. 4. 4.

Inside the ILR, although we do not have an explicit bar component in the potential, the particles feel an m=2 term due to the existence of the massive spirals starting, by construction, at the ILR. This leads to a bar-like response. We find x1v1-like orbits supporting a 3D peanut structure in the central region of the models. The backbone of the p.o. inside the ILR, projected on the equatorial plane, gives always a continuation of the spirals in the central part with an abrupt increase of their pitch angle. 5. 5.

Just after the radial 4:1 resonance we observe in general a gap at the spiral region in the response models. However, between 4:1 and corotation we find in the PERLAS potential local enhancements of the spiral arms, i.e. segments of arms, close to the imposed spiral loci. Similar local enhancements are not found beyond corotation. If we take into account the exponential profile of the PERLAS model these segments of arms are found in sparsely populated regions of the galaxies. Thus, it is unlikely to be part of the strong symmetric arms observed in grand-design, normal spirals. Until now, extensions of the arms beyond the 4:1 resonance were found in gaseous response models (see e.g. figure 4 in Patsis et al., 1997b), while weak extensions have been identified in $N$-body stellar models (see e.g. figure 10 in Patsis & Kaufmann, 1999). 6. 6.

In our response models, especially in those with large $M_{s}/M_{d}$ ratios, we observe “chaotic” spirals emerging from the L1, L2 regions. In such a case in the models act the dynamical mechanisms described in Patsis & Tsigaridi (2017) that support in the same model an inner and an outer spiral structure. These outer spirals are in regions of low density in the models.

L.C.V. and I.P. thank the Mexican Foundation CONACYT for grants that supported this research. L.C.V. thanks the Research Center for Astronomy (RCAAM) of the Academy of Athens for its hospitality during his visit there, when part of this work has been completed. LCV thanks IA-UNAM for their hospitality during his visit there. This work is part of the project “Study of stellar orbits and gravitational potentials in galaxies, with numerical and observational methods” in which researchers from RCAAM and INAOE participate.