Novel results from an algebraic approach to molecular bending dynamics

Novel results from an algebraic approach to molecular bending dynamics Francisco P´erez-Bernal1, Osiris ´ Alvarez-Bajo1, Jos´e M Arias2, Miguel Carvajal1, Jos´e E Garc´ıa-Ramos1, Danielle Larese3, and Pedro P´erez-Fern´andez2 1Departamento de F´ısica Aplicada, Fac. de CC. Experimentales, Universidad de Huelva, Huelva 21071 (SPAIN) 2Departamento de F´ısica At´omica, Molecular y Nuclear, Facultad de F´ısica, Universidad de Sevilla, Apdo. 1065, Sevilla E-41080 (SPAIN) 3Department of Chemistry, Yale University, P.O. Box 208107, New Haven, CT 06520-8107 (USA) E-mail: [email protected] Abstract. We present a brief review of research topics of current interest that depend on an algebraic approach to molecular bending dynamics. This approach is based on a u(3) spectrum generating algebra. In particular, we briefly present results on three topics: the calculation of finite-size analytical corrections to mean field results, the application of the model to the largeamplitude vibrational bending mode of the NCNCS molecule, and the analysis of the influence of quadratic Casimir operators on excited state quantum phase transitions. 1. Introduction The modeling of n-dimensional systems using a u(n+ 1) Lie algebra as the system’s spectrum generating algebra (SGA) has proved successful in an ample set of physical situations [1]. This approach has caused a major impact in nuclear structure studies, where the interacting boson model (IBM), based on a u(6) SGA and formulated by Iachello and Arima in the seventies, has become a standard theoretical tool [2]. Following the trail of the IBM, the Vibron Model was born in the eighties, treating molecules as bosonic systems where bosons represent vibrational and rotational excitations (vibrons) [3]. The SGA in this case is a u(4) Lie algebra, due to the three-dimensional nature of the problem. The Vibron Model, though useful, proved quite cumbersome when applied to molecular species composed by more than four atoms. In the cases that do not require a simultaneous treatment of molecular rotations and vibrations, a possible alternative approach is the onedimensional (1D) limit of the vibron model, based on a u(2) SGA [4]. In this limit rotational degrees of freedom are neglected and each local vibrational degree of freedom has an associated u(2) Lie algebra. Therefore this model considers anharmonicity from the outset, due to an isomorphism between the u(2) Lie algebra and the 1D Morse potential [5]. In addition, this approach naturally embeds discrete symmetry considerations [6]. The bending vibrational dynamics of linear and quasi-linear molecules imply the interplay of vibrational and rotational degrees of freedom. In this case an appropriate model is the twodimensional (2D) limit of the vibron model, with a u(3) Lie algebra as the system’s SGA [7]. GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 Published under licence by IOP Publishing Ltd 1 The 2D limit of the vibron model has successfully modeled spectral term energies and intensities of various molecular species, from rigidly-linear to rigidly-bent molecules, including quasi-linear, and non-rigid intermediate situations [8, 9, 10]. A thorough description of both quantal and classical aspects for this model can be found in Ref. [11]. In recent years the study and characterization of quantum phase transitions (QPTs) in algebraic models has attracted much attention, a line of research that can be traced back to the seminal work of Gilmore [12]. These transitions, also called shape phase transitions and ground state transitions, are zero-temperature phase transitions between different geometric limits of the system’s ground state upon variation of one or several control parameters in the Hamiltonian. This research line has been recently bolstered up with the discovery of precursors of QPTs in nuclear systems. Reference [13] is a excellent and recent review on this topic. More recently, the study of excited state quantum phase transitions (ESQPTs), transitions that take place upon the variation of the system’s excitation energy instead of the Hamiltonian parameters, has been proposed [14]. In this case the transition happens for some particular excited states of the system. This is especially relevant in the field of molecular spectroscopy, where novel spectroscopic techniques allow the access to highly-excited states. In fact, quantum monodromy effects, whose appearance in several molecular systems has been experimentally confirmed [15, 16], can be regarded as an ESQPT from the algebraic model perspective [11]. We finish this section with a brief outline of the present work. Section 2 is devoted to a short presentation of the bosonic algebraic approach to 2D systems and its classical limit. Analytical corrections to the classical limit of the model, going beyond its mean field limit are presented in section 3. An application to a realistic case can be found in section 4, where we fit vibrational bending energies of a large amplitude bending mode of NCNCS. Section 5 contains some interesting results concerning the influence of quadratic Casimir operators on the model ESQPT. In section 6 we finish by giving some concluding remarks. 2. Algebraic approach to 2D systems In this section we present, in an abridged form, the 2D limit of the vibron model and its classical limit. 2.1. Basic results The building bricks for the u(3) Lie algebra used for modeling 2D systems are three bosons of two different types. One of them is a scalar boson, σ†, and the other two are Cartesian bosons, {τ† x, τ† y}. The operators follow the usual bosonic commutation relations: [σ, σ†] = 1 and hτi, τ† ji=δi,j with i, j =x, y. All other commutators are zero. For the sake of convenience we introduce circular bosons [11] τ† ±=∓τ† x±iτ† y √2, τ±=∓τx∓iτy √2.(1) The nine u(3) generators are built with the possible bilinear products of σ†and τ† x,τ† ycreation and annihilation operators. They are conventionally expressed as [7]: ˆn=τ† +τ++τ† −τ−,ˆns=σ†σ ˆ l=τ† +τ+−τ† −τ− ˆ D+=√2(τ† +σ−σ†τ−),ˆ D−=√2(−τ† −σ+σ†τ+) ˆ Q+=√2τ† +τ−,ˆ Q−=√2τ† −τ+ ˆ R+=√2(τ† +σ+σ†τ−),ˆ R−=√2(τ† −σ+σ†τ+). (2) GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 2 The operator ˆ lis the 2D angular momentum, as one can easily see once it is expressed in terms of Cartesian boson operators. There are two possible dynamical symmetries that conserve 2D angular momentum, starting in the u(3) SGA and ending in the so(2) symmetry algebra. u(3) ⊃u(2) ⊃so(2) Chain I ,(3) u(3) ⊃so(3) ⊃so(2) Chain II .(4) Numerical calculations in the present article have been carried out using the basis associated with dynamical symmetry I. A detailed discussion of both dynamical symmetries can be found in [11]. Another ingredient of the algebraic approach is the Casimir (or invariant) operators associated to each subalgebra chain [5]. The first and second order Casimir operators of the subalgebras in chains (3) and (4) are ˆ C1[u(2)] = ˆn , ˆ C2[u(2)] = ˆn(ˆn+ 1) , ˆ C2[so(3)] = ˆ W2= ( ˆ D+ˆ D−+ˆ D−ˆ D+)/2 + ˆ l2, ˆ C1[so(2)] = ˆ l , ˆ C2[so(2)] = ˆ l2. (5) The so(3) Casimir operator ˆ W2can be replaced by the pairing operator ˆ P=N(N+1)−ˆ W2, where the total number operator, ˆ N= ˆns+ ˆn, has been replaced by its value, N, due to the fact that we consider systems with a fixed number of bosons. The most general rotation and parity-invariant, one and two-body Hamiltonian can be written as a linear combination of the four possible first and second order Casimir operators (5) ˆ H=E0+ǫˆ C1[u(2)] + αˆ C2[u(2)] + βˆ C2[so(2)] + Aˆ C2[so(3)] .(6) Matrix elements in the two possible bases for the Casimir operators (5) can be found in Ref. [11]. 2.2. Mean field limit of the model The mean field limit of the model, also called the thermodynamical limit, the large Nlimit, or the classical limit of the model, can be found following an algorithm established by Gilmore [17] known as the coherent (or intrinsic) state approach. In making use of the coherent state approach, we use projective coherent states that define an intrinsic ground state |[N]; r, θi=1 √N!b† cN|0i,(7) where rand θare polar coordinates associated to Cartesian coordinates xand y, and b† cis the boson condensate b† c=1 √1 + r2σ†+xτ† x+yτ† y.(8) The coherent state (7) is the number projected generalized coherent state of u(3) [2]. The coherent state approach provides an energy functional that is an approximation for the ground state energy of a given Hamiltonian. This approximation becomes exact for N→ ∞. The existence of a second order QPT between the two possible dynamical symmetries can be easily proven using the aforementioned coherent state approach [11]. In order to do so, it is convenient to define a simplified model Hamiltonian, with one operator for each chain and properly normalized ˆ H=ε(1 −ξ)ˆn+ξ N−1ˆ P.(9) GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 3 The parameter εis an energy scale and ξis a control parameter that takes values between zero (u(2) dynamical symmetry) and one (so(3) dynamical symmetry). The energy per boson for the essential Hamiltonian (9) is [9, 11] Eξ(r)≡1 Nh[N]; r, θ|ˆ H|[N]; r, θi h[N]; r, θ|[N]; r, θi=ǫ (1 −ξ)r2 1 + r2+ξ 1−r2 1 + r2!2 .(10) The extremization of the energy functional (10) reveals that there are two different geometric limits: symmetric (or linear) and deformed (or bent). In the symmetric case the energy minimum is at re= 0, while in the second case re=p(5ξ−1)/(3ξ+ 1) 6= 0. The symmetric phase takes place for control parameter values ξ≤ξc= 0.2, and the deformed phase for ξ > ξc= 0.2. When evaluated at r=re, the energy functional is Eξ(re) = (ξ0≤ξ≤ξc −9ξ2+10ξ−1 16ξξc< ξ ≤1.(11) By evaluating the derivatives of Eξ(re) with respect to ξ, one finds that, according to Ehrenfest’s classification scheme, the u(2) −so(3) phase transition is of second order [11]. Another observable of interest is the expectation value of the u(2) number operator in the system ground state h[N]; r, θ|ˆn|[N]; r, θi=Nr2 e 1 + r2 e =(0 0 ≤ξ≤ξc 5ξ−1 8ξξc< ξ ≤1.(12) This operator can act as a classical order parameter for the transition between the symmetric and deformed phases [11]. 3. Beyond mean field analytical corrections We evaluate analytical corrections to the mean field limit of the algebraic approach to go beyond the results presented in the previous section. The mean field or classical limit is only exact in the large Nlimit, and the calculated corrections allow us to calculate the N0corrections to the ground state energy (11) and the number of τbosons (12). In order to compute these corrections we perform a Holstein-Primakoff expansion and a shift transformation, followed by a Bogoliubov transformation, following Ref. [18]. This reference is a general work, aimed at bosonic two-level Hamiltonians with SGA u(2L+ 2) with L= 1,2,.... In the symmetric region Dusuel and collaborators go a step further than us, making use of the continuous unitary transformations approach (CUTS) [19]. The Holstein-Primakoff expansion implies the definition of a new pair of Cartesian bosons, b† i(bi), i=x, y, with the usual bosonic commutation relations, [bi, b† j] = δij. The inclusion of this boson removes the dependence on the scalar σ†boson of the Hamiltonian τ† iτj=b† ibj, τ† iσ=√N b† iq1−ˆnb/N =σ†τi†,(13) ˆnσ=σ†σ=N−ˆnb, where i, j =x, y, and ˆnb=b† xbx+b† yby. A third set of bosons, c† i(ci), with i=x, y is defined via a shift transformation b† i=√Nλδix +c† i;i=x, y . (14) GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 4 The parameter λ, known as the shift parameter, is zero in the spherical phase and nonzero in the deformed phase. Once the model Hamiltonian (9) is transformed, we obtain an expansion in powers of N, with terms ˆ H=ˆ H1+ˆ H1/2+ˆ H0+O(1/√N),(15) where ˆ Hiincorporates the terms with a Nidependence. The first term is ˆ H1=Nhξ+ (1 −5ξ)λ2−4ξλ4i.(16) Setting λ=r/√1 + r2this term is nothing more than the mean field result (10). The next order in the expansion of the model Hamiltonian, ˆ H1/2, is zero once the equilibrium values of r(or λ) are substituted. This is due to the relation ˆ H1/2=2 √N dˆ H1 dλ . Therefore, the first finite-size correction to the mean field limit is the ˆ H0term, that is quadratic in the cboson operators H0=λ2(4λ2−1)ξ+ (1 −3ξ+ 14λ2ξ)c† xcx+ (1 −3ξ+ 4λ2ξ)c† ycy + (5λ2−1)ξc† xc† x+cxcx+ (2λ2−1)ξc† yc† y+cycy.(17) The Hamiltonian (17) can be diagonalized with a Bogoliubov transformation c† i=uia† i+viai, ci=uiai+via† i,(18) that should be independently performed in the symmetric and deformed phases [18, 20]. The final result provides the correction to the ground state energy per particle in the symmetric and deformed regions Esym 0=ξ+N−1h3ξ−1 + Ξsym(ξ)1/2i,(19) Edef 0=−9ξ2+ 10ξ−1 16ξ+N−1"1−6ξ−27ξ2+ 8ξΞdef (ξ)1/2 16ξ#,(20) where Ξsym = 5(ξc−ξ)(1 −ξ) and Ξdef (ξ) = 5(ξ−ξc)(1 + 3ξ). In the deformed phase the y coordinate contribution is a spurious Goldstone boson, associated with a ground state rotation. Compared to the mean field limit, the deduced finite-size corrections greatly improve the agreement with numerical results. This can be clearly seen in Fig. 1 where the mean field and beyond mean field (BMF) corrections for the ground state energy are compared to a numerical calculation for N= 20. The finite-size correction can be calculated for other observables. We include results for the expectation value in the ground state of ˆn(5). In this case we can make use of the Hellman- Feynman theorem, defining a new control parameter xsuch that ξ= 1/(1 + x), we obtain for the symmetric and deformed phases hˆni N=d dx [(1 + x)Esym 0(x)] = 1 N 1−3ξ−Ξsym(ξ)1/2 Ξsym(ξ)1/2+O(N−2),(21) hˆni N=d dx h(1 + x)Edef 0(x)i=5ξ−1 8ξ+1 N 4ξ(ξ−1) + (1 −3ξ)Ξdef (ξ)1/2 8ξΞdef (ξ)1/2+O(N−2).(22) We compare this result, the mean field limit and the numerical calculation for N= 20 in Fig. 2. As in the previous case, the obtained correction notably improves the mean field result, in particular in the vicinity of the critical control parameter. GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 5 0.0 0.2 0.4 0.6 0.8 1.0 Control Parameter (ξ) 0.00 0.05 0.10 0.15 0.20 0.25 Ground State Energy per Particle (ε0=E0/N) Mean Field BMF Symmetric N = 20 BMF Deformed N = 20 Numerical, N = 20 Figure 1. (Color online) Ground state energy per particle (E0) in arbitrary units as a function of the control parameter ξ. 0.0 0.2 0.4 0.6 0.8 1.0 Control Parameter (ξ) 0.0 0.1 0.2 0.3 0.4 0.5 Normalized Order Parameter (〈n(ξ)〉/N) Mean Field BMF Symmetric N = 20 BMF Deformed N = 20 Numerical N = 20 Figure 2. (Color online) Expectation value of the number operator ˆnin the ground state of Hamiltonian (9) as a function of the control parameter ξ. The results presented are encompassed in an open line of research [21]. We are currently working on computing the BMF correction for other observables, on analytically deriving finitesize scaling exponents for the u(2) −so(3) second order phase transition, and on establishing a connection between Ref. [18] results and ours. 4. Application to the ν7bending mode of cyanogen isothiocyanate The bending dynamics of several molecular species [7, 9, 10, 22, 23, 24, 25] have been modeled using the two-dimensional limit of the vibron model. We are currently interested in the analysis of the spectra of non-rigid molecules, as these probably provide examples of ESQPT. As the excitation energy increases, the non-rigid molecule’s excited states have to overcome a potential hump in the origin, and change from a bent-like to a linear-like character. This fact was known long ago [26], but nowadays such systems have focused an increasing degree of attention due to quantum monodromy effects and its implications [16]. We have performed a fitting to the bending energy spectrum of cyanogen isothiocyanate (NCNCS), a non-rigid molecule whose large amplitude bending spectrum displays monodromy effects and is experimentally accessible [27]. We fit Hamiltonian (6) to the NCNCS term values obtained with the General Semirigid Bender Hamiltonian model in Ref. [27]. The final rms is 2.2 cm−1, in a fit to 70 term energies. The optimized parameters are given in Table 1. The calculations were performed using the FORTRAN code triat u3 [28]. Table 1. Optimized parameters of the one- and two-body Hamiltonian (6) and their associated uncertainties obtained in the fit to ν7term values of NCNCS [27]. All parameters, excepting N, are in units of cm−1. The fit root mean square deviation is rms = 2.2 cm−1. N ǫ α β A 70 203.6(19) -2.58(3) 1.496(10) 0.813(6) The obtained results are promising, especially if the simplicity of the model is considered. GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 6 The results in Ref. [27] and our fit are compared in a quantum monodromy plot in Fig. 3 where the good agreement between both approaches is remarkable. We reproduce the change in slope in the origin that characterizes the critical monodromy point. -8 -6 -4 -2 0 2 4 6 8 Vibrational Angular Momentum (K) 0 200 400 600 800 1000 Energy (cm-1) This work PRL 95 243002 Figure 3. (Color online) Quantum monodromy plot of calculated ν7mode bending vibrational levels of NCNCS. The vibrational angular momentum lis labelled Kfollowing the usual spectroscopic notation. We are currently working in the fit of the vibrational bending spectrum of water, where quantum monodromy effects have been experimentally reported [15]. 5. Influence of the quadratic u(2) Casimir operator in ESQPT As shown in Table 1, the second order Casimir operator of the u(2) subalgebra plays an important role when fitting experimental data. This also happens in other cases [9, 10] and has suggested us to explore the role of this anharmonic term in the model, in particular its effects on ground state QPT and ESQPT [29]. This can be accomplished with a new model Hamiltonian, conveniently scaled, that incorporates a second control parameter ˆ H=ε(1 −ξ)ˆn+α N−1ˆn(ˆn+ 1) + ξ N−1ˆ P.(23) The coherent state approach, when applied to this case, gives as a result that the ground state transition is mostly unperturbed by the addition of the new term. The energy functional E0(ξ, α) in the mean field limit becomes E0(ξ, α;r) = εξ+ (1 −3ξ)r2+ (1 + α)r4 (1 + r2)2.(24) The minimization of the energy functional (24) provides the equilibrium values re= 0 ,q5ξ−1 3ξ+2α+1 that implies a lower bound α > −(1+3ξ)/2 if ξ > 0.2. The critical value of the control parameter is, once more, ξc= 0.2. For ξ≤ξcthe energy functional minimum lies at the origin, while for values of ξ > ξca new minimum appears and the minimum at the origin becomes a maximum. When evaluated for r=reand ε= 1, the energy functional is E0(ξ, α;re) = (ξ0≤ξ≤ξc −9ξ2+10ξ+4αξ−1 16ξ+4αξc< ξ ≤1,(25) which has a discontinuous second order derivative in ξ=ξc. Although the second order ground state transition is basically unaffected by the new control parameter, there are conspicuous effects in the system’s ESQPT for negative αvalues. This GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 7 is specially noticeable once the excitation energy diagram for the Hamiltonian (23) is plotted for different αvalues. Instead of one separatrix between the configuration with linear and bent character, like in the α= 0 case [11], there are two separatrices. In Figs. 4 and 5 we plot the excitation energy diagram of zero angular momentum (l= 0) states for and α= 0 and −0.4, respectively. The number of bosons in both cases is N= 100 and the energy is plotted as a function of the control parameter ξ. The color code in the curves is assigned in accordance to the magnitude of the maximum squared-component of the eigenvector in the two bases of choice. If the eigenvector is close to dynamical symmetry (3), the maximal component will be close to one in the u(2) basis and the point will be plotted in red color. If the maximal component is in the so(3) basis the mapped color is blue. This way of representing the data clearly marks the states that correspond to the separatrix in the monodromy plot, in the high level density region that separates the linear- and bent-like states. The appearance of two separatrices in Fig. 5 is markedly clear. Figure 4. (Color online) Excitation energy diagram in arbitrary units (ε= 1) for zero angular momentum eigenstates of Hamiltonian (23) as a function of the control parameter ξ for α= 0. The number of bosons is N= 100 and energies are normalized by N. The colorcode indicates the state proximity to the u(2) or so(3) dynamical symmetries (see text). Figure 5. (Color online) Excitation energy diagram in arbitrary units (ε= 1) for zero angular momentum eigenstates of Hamiltonian (23) as a function of the control parameter ξfor α=−0.4. The number of bosons is N= 100 and energies are normalized by N. The colorcode indicates the state proximity to the u(2) or so(3) dynamical symmetries (see text). The position of the two separatrices in Fig. 5 depends on the values of the energy functional (25) maximum at the origin and its asymptotic value [29]. The value at the origin is only a function of ξ, while the asymptotic value depends only on α. In particular, there is a critical value of α,αc=ξ−1, where the two separatrices cross and it is associated with equal values of the energy functional at zero and for large rvalues. In the case of ESQPT there is not a clear definition of an order parameter [14]. The expectation value of the number operator ˆnin the different eigenstates provides a possible approximation. This observable changes abruptly in the critical excitation energy but it does not go to zero in any of the phases [14]. We plot in Fig. 6 the expected value of ˆnfor the eigenstates of Hamiltonian (9) as a function of the normalized state excitation energy. We have GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 8 fixed the number of bosons (N= 800) and the control parameter ξ(ξ= 0.6). For α= 0 (red curve in Fig. 6) we recover the results published in [14], where the critical energy is marked by the warped minimum in the plotted function. For 0 > α > αc, in addition to the displacement to larger energies of the initial minimum, a maximum appears for larger excitation energies (orange curve in Fig. 6). As αtends to αcboth extremes approach. We found the surprising result that for α=αcthe expectation value of ˆnbecomes constant for the full energy range except a sudden variation around the critical energy that corresponds to the crossing point of the two separatrices (green curve in Fig. 6). 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Excitation Energy (E/N) 0.2 0.4 0.6 0.8 1.0 Order Param. [〈n(E)〉/Ν] N = 800, ξ = 0.6 α = 0.0 α = -0.25 α = -0.4 Figure 6. (Color online) Expectation value of the number operator ˆn in the eigenstates of model Hamiltonian (23) for N= 800 and ξ= 0.6 as a function of the excitation energy of the state normalized by N for three different values of α. 6. Concluding remarks From our point of view, the two dimensional limit of the vibron model is a simple algebraic model that still encompasses enough complexity to be applied to interesting real physical systems. The calculation of finite-size corrections to the mean field limit greatly deepens the knowledge of the model. These corrections are a fundamental help to characterize the QPTs between the different geometric limits of a model and to analytically extract the scaling exponents for the model, as it will be shown for the model under study in a forthcoming publication [21]. We have applied the model to the bending dynamics of a non-rigid molecule. Despite the simplicity of the approach, it is capable of coping with complicated situations, as it has been shown in the case of NCNCS, where quantum monodromy effects appear. We are especially interested in the description with this model of water’s vibrational bending levels. Being the simplest bosonic two-level model with a non-trivial angular momentum, this limit of the vibron model is a very convenient aid in the study of ESQPTs and their implications. The experimental access to bending excitation energies above the monodromy critical energy makes this feature especially attractive. The obtainment of experimental data for such states is not possible in other fields, as nuclear spectroscopy, which have been the traditional test ground for ground state transitions. Acknowledgments This work was supported in part by the Spanish Junta de Andaluc´ıa under projects P07- FQM-02962, P07-FQM-03014, and P07-FQM-02894 and by the Spanish MICINN and the European regional development fund (FEDER) under projects FIS2008-04189 and CPAN- Ingenio (CSD2007-00042). PPF acknowledge the Spanish MEC for a FPU grant. The authors thank Jorge Dukelsky, Francesco Iachello, Pieter van Isacker, and Amiram Leviatan for valuable comments and suggestions. GROUP 28: Physical and Mathematical Aspects of Symmetry IOP Publishing Journal of Physics: Conference Series 284 (2011) 012049 doi:10.1088/1742-6596/284/1/012049 9