Lauric Thiault
email: lauric.thiault@gmail.com

Stefan Gelcich

Nadine Marshall

Paul Marshall

Frédérique Chlous

Joachim Claudet

Operationalizing vulnerability for social-ecological integration in conservation and natural resource management

Keywords: Co-benefits, human-nature interactions, interventions portfolio, social-ecological fit, trade-offs, resilience

come

Introduction

Achieving sustainability on our overexploited planet is one of the grand challenges of our time [START_REF] Rockström | A safe operating space for humanity[END_REF]. This global challenge has local expressions that are both social and ecological in form, because people and nature are linked and interdependent (Fischer et al. 2015). Such strong social-ecological relationships are especially apparent in resourcedependent settings such as forestry communities or coastal fisheries, where unsustainable use of natural resources can lead to serious and tangible impacts on both ecosystems and the people that depend on them [START_REF] Ostrom | A general framework for analyzing sustainability of social-ecological systems[END_REF]IPBES 2019).

Many governmental agencies and NGOs are beginning to embrace a more nuanced view of sustainability that sits at the nexus between social and ecological perspectives [START_REF] Bakker | The changing landscape of conservation science funding in the United States[END_REF][START_REF] Díaz | The IPBES Conceptual Framework -connecting nature and people[END_REF]. As a result, strategies aiming to improve conservation and social outcomes increasingly incorporate both elements in design and implementation [START_REF] Mace | Whose conservation? Science[END_REF]. Indeed, through initiatives such as multiple-use protected areas and ecosystem-based management, social considerations are now embedded in the design of many "eco-centric" measures, hence broadening a predominantly ecological view of conservation and natural resource management [START_REF] Ban | A social-ecological approach to conservation planning: embedding social considerations[END_REF][START_REF] Kittinger | A practical approach for putting people in ecosystem-based ocean planning[END_REF]. Correspondingly, the sustainable livelihood approach illustrates how the human development community, whose "socio-centric" entry-point has been predominantly centered around reducing poverty or fostering development opportunities, increasingly recognizes good environmental status as part of the conditions affecting the success of interventions [START_REF] Krantz | The Sustainable Livelihood Approach to Poverty Reduction[END_REF][START_REF] Roe | Are alternative livelihood projects effective at reducing local threats to specified elements of biodiversity and / or improving or maintaining the conservation status of those elements ?[END_REF][START_REF] Wicander | Can the Provision of Alternative Livelihoods Reduce the Impact of Wild Meat Hunting in West and Central Africa ?[END_REF].

Integration of a social-ecological science perspective into human development, conservation and natural resource management has enhanced the long-term equitability and effectiveness of the initiatives of each. Moreover, decades of applications of eco-and sociocentric strategies in various settings have offered important insights and experience that provide valuable foundations upon which more integrated, cross-disciplinary approaches can be built. Although still imperfect, we now have a better understanding of what can work and what cannot, in what contexts, why, and how to avoid potentially undesirable outcomes [START_REF] Barrett | Institutional arrangements for rural poverty reduction and resource conservation[END_REF][START_REF] Cox | A Review of Design Principles for Communitybased Natural Resource Management[END_REF][START_REF] Wright | Reframing the concept of alternative livelihoods[END_REF][START_REF] Barnes | Understanding local-scale drivers of biodiversity outcomes in terrestrial protected areas[END_REF][START_REF] Wicander | Can the Provision of Alternative Livelihoods Reduce the Impact of Wild Meat Hunting in West and Central Africa ?[END_REF]. Despite these positive developments, responses to sustainability problems continue to be dominated by strategies focusing mostly on either the human or environmental elements of the social-ecological systems.

Successfully dealing with conservation and sustainability requires a diverse portfolio of interventions. Therefore, the challenge now is to stop striving for eco-centric or socio-centric strategies, and instead seek synergies of the two. Indeed, and although they may diverge in many ways, eco-centric and socio-centric approaches are often complementary: when well designed, eco-centric interventions can enhance elements of human well-being, and sociocentric interventions can improve ecological condition [START_REF] Mcclanahan | Conservation action in a changing climate[END_REF][START_REF] Roe | Are alternative livelihood projects effective at reducing local threats to specified elements of biodiversity and / or improving or maintaining the conservation status of those elements ?[END_REF][START_REF] Ban | Well-being outcomes of marine protected areas[END_REF][START_REF] Naidoo | Evaluating the impacts of protected areas on human well-being across the developing world[END_REF]). Yet, neither intervention is likely to provide a "silver bullet" [START_REF] Ostrom | Going beyond panaceas[END_REF]. Instead, we should be looking for a "silver buckshot", where several tools in the box are used.

Insights offered by social-ecological systems thinking, and the extensive and mature knowledge supporting human development, natural resource management, and conservation together provide momentum for developing and institutionalizing a new generation of management practices that positions the links between people and nature at its core. Here we aim to address the narrower, but still difficult challenge of improving integration across independent but complementary sustainability-seeking strategies while ensuring relevance to decision-makers and practitioners. To do so, we have developed an approach based on "vulnerability profiles", which represent the system's social and ecological elements that are favoring or undermining sustainability, thus revealing the internal features that can most effectively be targeted by sustainability interventions. This approach ultimately makes apparent a portfolio of interventions that can help realize cobenefits across goals relating to conservation, food security and human well-being. We illustrate our approach using the case of a small-scale coral reef fishery in French Polynesia, where fishing activity represents both an invaluable source of benefits for local communities and an important pressure on the ecosystem.

A vulnerability-based approach for integrated management of social-ecological systems

The approach we present here draws on recent developments in vulnerability and socialecological system thinking [START_REF] Cinner | Evaluating social and ecological vulnerability of coral reef Fisheries to Climate Change[END_REF]Thiault et al. 2018a; Appendix A). It involves a four-step procedure that leads to the identification of practical interventions that most appropriately echoes the needs and opportunities of a particular social-ecological system (Fig. 1). It is intended to serve as an operational guide for the place-based management of resource and resource user interactions, where ecological vulnerability refers to the vulnerability of the resource (e.g. water, wild food, landscape) to use by the resource users (e.g. farmers, fishers) and social vulnerability refers to the vulnerability of the resource users to use-induced resource degradation. Therefore, it does not necessarily aim to address all drivers of change in the social-ecological system of interest. It assumes that the system's boundaries have been identified and that analysts aspire to achieve social and ecological outcomes.

Step 1: Assessing resource and resource user interactions through the lens of vulnerability

The approach first guides analysts to independently assess each key dimension of socialecological vulnerability [START_REF] Cinner | Evaluating social and ecological vulnerability of coral reef Fisheries to Climate Change[END_REF]Thiault et al. 2018a), namely resource exposure, sensitivity and adaptive capacity to exploitation (ecological vulnerability), and users' exposure, sensitivity and adaptive capacity to resource decline (social vulnerability). Social exposure is determined by ecological vulnerability (Appendix A; Fig. S1) and thus does not need to be assessed explicitly. Ecological sensitivity and adaptive capacity can be difficult to untangle because they are determined by similar processes. Here, we refer to their combination as "intrinsic resilience" but acknowledge that resilience entails far more complex processes that are not captured by this model. The four remaining dimensions can then be combined to allocate ecological and social components to one of four quadrants, hereafter referred to as vulnerability profiles (Fig. 2). Profiles are labelled as "lower concern", "potential adapter", "high latent risk" and "greater concern" and characterize the main elements that best determine social and ecological vulnerabilities, highlighting what needs to be targeted to reduce vulnerability. Analysists can characterize vulnerability profiles in the way that it best suited to social-ecological context (Appendix B).

Step 2: Selecting interventions that can reduce source(s) of vulnerability

Step 2 involves identifying relevant interventions that target the features identified in the 1). Analysts might look to this template as a starting point, adapting and rearranging as necessary.

Step 3: Prioritizing interventions that can advance social-ecological co-benefits

Step 3 entails being critical of negative collateral impacts that some interventions might have, as well as employing those initiatives that benefit each system [START_REF] Sayer | Ten principles for a landscape approach to reconciling agriculture, conservation, and other competing land uses[END_REF][START_REF] Howe | Creating win-wins from trade-offs? Ecosystem services for human well-being: A meta-analysis of ecosystem service tradeoffs and synergies in the real world[END_REF]. Indeed, in order to be successful and balanced, management interventions identified in step 2 must be appropriately positioned in the social and ecological context in such a way that they do not further undermine any component of the system. Instead, they should be employed to reduce negative impacts and/or induce positive change. In Table 1, we summarize how various types of commonly used interventions implemented on one component may have indirect effects on others, and how this can be interpreted using the ecological and social vulnerability profiles from step 1. Like in the previous step, this template can be further tweaked to local contexts.

Step 4: Developing an interventions portfolio that suits the broader social-ecological landscape The last step captures the wider context in which the local interactions between resource and resource users are embedded to ensures the feasibility and viability of previously identified interventions . This step includes documenting the social norms, values, cultural practices, aspirations, place attachment and historical and environmental features that can facilitate or hinder specific interventions [START_REF] Ostrom | A general framework for analyzing sustainability of social-ecological systems[END_REF][START_REF] Armitage | Environmental governance and its implications for conservation practice[END_REF]. To ensure the sustainability and legitimacy of selected interventions, traits of the governance structure need to be accounted for (Leenhardt et al 2015). Mixed methods approaches and triangulation of qualitative and quantitative data from various sources (Appendix B; [START_REF] Game | Cross-discipline evidence principles for sustainability policy[END_REF]) can create a cohesive picture that will help analysts assess whether each candidate intervention is appropriate, equitable and legitimate [START_REF] Kittinger | A practical approach for putting people in ecosystem-based ocean planning[END_REF]).

Our vulnerability-based approach thus consists of four steps eventually leading to the selection of one or more interventions that are important and actionable to reduce social and ecological vulnerabilities (Fig. 3). By effectively considering the linkages between key social and ecological components, it enables to identify management strategies that are more likely to meet both ecological and social goals than if one criterion only was considered. As such, it offers practical insights that can inform integrated management strategies and planning in a broad range of contexts.

Illustrating the approach: a coral reef fishery case study

We use the coral reefs and the associated small-scale fishery of Moorea, French Polynesia, to illustrate the application of the approach described above. Overall, the Moorea fishery is highly challenging to manage due to inextricable yet diffuse links between people and the reef [START_REF] Leenhardt | Complexities and uncertainties in transitioning small-scale coral reef fisheries[END_REF]. The marine spatial plan in which fisheries management is embedded was under revision at the time this study was conducted [START_REF] Hunter | Maneuvering towards adaptive co-management in a coral reef fishery[END_REF]), and our pilot assessment was conducted in parallel of the revision process.

In order to consider linked social-ecological vulnerabilities in the specific context of fish (the resource) and fishing households (the resource users), we compiled data on marine resource dependency (i.e. social sensitivity) and adaptive capacity from 6,698 households, and combined it with reef-wide models of target fish assemblages, characterized by their intrinsic resilience and exposure to fishing. The combination of each dimension of social and ecological vulnerabilities was represented spatially to visualize the vulnerability profiles (step 1; Fig. 4). We then applied the general typology of eco-and socio-centric management interventions (Table 1) into the context of small-scale fisheries (steps 2-3; Table 2). Finally, we used a combination of archival research, semi-structured interviews from key informants, and participant observations to gain insights into the broader context, thus allowing us to better capture elements that could facilitate or hinder each potential intervention (step 4; Table 2). See Appendix B for a full description of the methods.

The current management approaches implemented in Moorea to manage local fisheries are not aligned with the approaches suggested by our approach. For example, the fore reef generally shows high intrinsic resilience and relatively low exposure to fishing (Fig. 1; Appendix C, Fig. S2). Our results suggest that such configurations may support the development of fully protected areas because these ecologically efficient but socially restrictive measures are easier to implement and represent lower opportunity costs for local households. Yet, despite the large permanent fisheries closure system (20% of the total reef area), the fore reef only represents 7.7% of the total area protected (Appendix C, Fig. S3). In contrast, lagoon areas closed to fishing are in some cases located in front of poorly adaptive, and sometimes highly sensitive households (Fig. 3; Appendix C, Fig. S3), creating a policy setting that could exacerbate social vulnerability and certainly lead to challenges for compliance. Given the criticisms against the current network of fully protected areas and their lack of ecological effectiveness (likely due to, in part, to poaching; [START_REF] Thiault | Ecological evaluation of a marine protected area network: a progressive-change BACIPS approach[END_REF], such conservation measures should be prioritized on the fore reef or in lagoon areas where associated people are weakly sensitive can adapt to the loss of fishing grounds (Fig. 4).

Where households are most vulnerable (e.g. Fig. 4b; Appendix C, Fig. S2), less restrictive interventions such as size and species regulations, or temporal closures, could be used to reduce fishing effort (i.e. ecological exposure) at a lower opportunity cost for users.

Although these types of interventions can be more difficult to enforce, and the perceptions on which are the best modalities can differ among stakeholders, they are generally supported by users and can be underpinned by preexisting legislation (Table 2). In parallel to addressing ecological exposure, ecological intrinsic resilience needs to be enhanced, particularly within the lagoon (Fig. 4). Although managers may for instance replicate previous stock enhancement interventions of targeted herbivores [START_REF] Taiarui | Is stock enhancement the best option to manage fisheries? A case study from Taiarapu (French Polynesia)[END_REF], improving the management of land-based activities are likely to have the greatest positive impact on ecological intrinsic resilience [START_REF] Leenhardt | Exploring social-ecological dynamics of a coral reef resource system using participatory modeling and empirical data[END_REF]. This type of approach is in line with principles from traditional 'ridge-to-reef' management, but its implementation would require greater collaboration among relevant agencies (Table 2).

In various locations around Moorea (e.g. Fig 4b-c), it is particularly relevant to couple the above eco-centric interventions with socio-centric ones focusing on the root cause of social vulnerability. This implies moving beyond stakeholder consultation processes to also investing in strategies that directly tackle social adaptive capacity and/or sensitivity. This may entail livelihood-focused measures such as incentives to diversify occupations (e.g. agriculture, tourism, aquaculture) and catch, although challenges regarding socio-cultural barriers need to be anticipated to avoid discrepancies between expectations and actual outcomes (Table 2). Community buy-in may for instance be leveraged via churches and other stakeholder groups, while land tenure issues can for example be overcome through enabling local community members to lease land cheaply for agricultural purposes. If well designed, and if new livelihoods are effectively created as alternatives rather than supplementary sources of outcome [START_REF] Wright | Reframing the concept of alternative livelihoods[END_REF], such interventions have the potential to reduce both social (reduced dependency, enhanced flexibility) and ecological vulnerably (released pressure on the resource). Enhancing adaptive capacity through social capital building, and encouraging learning and cooperation may, in Moorea, build on established stakeholder groups like cultural associations, while the recently created decentralized management committees provide an obvious forum for discussion on reef-related issues and solutions within the community (Table 2). Investments in market-based interventions and insurance schemes do not seem applicable for Moorea due to the absence of markets and the difficulty of identifying individual fishers (Table 2). Island-wide, and in particular in high socially sensitive areas (Fig. 4b-c), it is essential to develop strategies that do not make local communities more dependent on reef-based resources that are already highly vulnerable. This is why island-scale incentives for motorized boats or new fishing gear should be avoided.

Instead of constraining decision-makers to a single strategy-focused approach defined a priori (i.e. eco-or socio-centric), our results compel decision-makers to consider multiple entry points. Although many challenges remain to ensure the success of Moorea's management [START_REF] Hunter | Maneuvering towards adaptive co-management in a coral reef fishery[END_REF], our results suggest that the current strategy could be upgraded by shifting from a focus exclusively on the resource to account more specifically for social-ecological linkages in each location, and embracing a broader range of management options that include eco-and socio-centric interventions.

Reflections on the vulnerability-based approach

Our four-step process represents a significant departure from more mainstream approaches to vulnerability conceptualization and practice. First, the framing is new. At its core, it builds on, and brings together insights from social-ecological science and vulnerability, moving the latter from its original perspective (natural hazards and climate) towards a sustainability one that includes people, both as a factor affecting environmental outcomes, and as a recipient of environmental benefits that require human well-being to be improved (Thiault et al. 2018a). Second, it builds on previous applications and uses of the vulnerability construct, providing guidance not only for prioritization, but also for real, pragmatic, and balanced interventions. Third, our framework effectively bridges independent, yet complementary approaches such as ecosystem-based management [START_REF] Levin | Integrated ecosystem assessments: Developing the scientific basis for ecosystem-based management of the ocean[END_REF], resilience-based management [START_REF] Mcleod | The future of resilience-based management in coral reef ecosystems[END_REF], and the sustainable livelihood approach [START_REF] Krantz | The Sustainable Livelihood Approach to Poverty Reduction[END_REF].

Based on our experience we believe that this approach is likely to infuse a more comprehensive vision into conservation and natural resource management [START_REF] Guerrero | Achieving the promise of integration in social-ecological research: a review and prospectus[END_REF], and empower practitioners to develop more integrated management strategies.

The spatial representation of the vulnerability profiles revealed potential interventions best suited to each location around the island, thus allowing local managers to examine previously unexplored, yet locally-relevant management possibilities. The approach leading to the selection of the interventions portfolio is transparent and can be replicated through time, providing a structure for implementing an adaptive management process that supports responsive strategies [START_REF] Kaplan-Hallam | Adaptive social impact management for conservation and environmental management[END_REF].

The use of vulnerability in a resource management context is relatively recent (REF), and its use as practical tool is still controversial (Appendix A, Table S1). Most commonly, critiques relate to the potential 'top-down' process of vulnerability assessments, where local communities' input into the process can be left aside [START_REF] Cameron | Securing indigenous politics: A critique of the vulnerability and adaptation approach to the human dimensions of climate change in the canadian arctic[END_REF]. In our approach, each step can rely on a community-based, participatory process, for instance by engaging stakeholders into the design and collection of indicators (step 1), the identification of candidate interventions (steps 2-3), or analysis of the overall context (step 4). This would not only be critical for improving the quality of the assessment, but may also promote opportunities to reflecting a richer knowledge that aligns with local people's perspectives and insights [START_REF] Reed | Stakeholder participation for environmental management: A literature review[END_REF][START_REF] Dacks | Developing biocultural indicators for resource management[END_REF]. In addition, this process -like any vulnerability assessment-underplays other elements such as design considerations and complex interactions between and among social and ecological components that also shape the success of a particular intervention [START_REF] Agrawal | Poverty, development, and biodiversity conservation: shooting in the dark? WCS Work[END_REF][START_REF] Wright | Reframing the concept of alternative livelihoods[END_REF][START_REF] Selig | A typology of fisheries management tools: using experience to catalyse greater success[END_REF]). Our approach is not meant to be prescriptive and should rather be the support of discussions within a planning process.

Conclusions

Achieving social and ecological objectives are intimately linked goals. Therefore, they should be integrated within the same framework. Our proposed approach illustrates that there is much scope for improved integration of data, ideas and management practices across various fields and disciplines. Although this will not solve all the challenges facing conservation and natural resource management, it broadens the range of options, offering more transparency in decision-support tools. Step 3).

  previous step. They could include interventions focusing on the resource ("eco-centric" interventions such as ecological engineering, permanent closures, or output controls), on resource users ("socio-centric" such as livelihood-focused interventions, market-based approaches, or assets enhancement), or a combination of those depending on the elements that need addressing. Analysts may be interested in implementing participatory mechanisms to develop this initial list of interventions. All options should be explored carefully for holistic management. To help in the screening process, we propose a typology of interventions commonly used by development, natural resource management and conservation communities, and describe their expected impacts on ecological and social vulnerability profiles (Table

  social and ecological vulnerability profiles. Symbols indicate the effect of interventions (• positive; ⊗ negative; ◯ no effect) on each vulnerability profile (Step 1: greater concern: purple; potential adapter: yellow; high latent risk: blue; lower concern: grey; see Fig. 2). Intervention types a-e: "ecocentric" interventions. Intervention types 1-8: "socio-centric" interventions. Clear boxes indicate direct effects on the component (e.g. effect of eco-centric interventions on the resource; Step 2) and shaded boxes indicate indirect effects (e.g. effect of eco-centric interventions on the resource users;

  ) to a small-scale coral reef fishery. Examples of interventions are presented, together with how they would be filtered in Moorea according the island's broader context (Step 4), which may facilitate or prevent successful implementation of particular interventions.

Figure legends Figure 1 |

 legends1 Figure legends

Figure 2 |

 2 Figure 2 | Social and ecological vulnerability profiles and associated management targets (Step 1). Each profile is identified through the combinations of exposure and intrinsic resilience gradients (ecological vulnerability), or sensitivity and adaptive capacity gradients (social vulnerability), and thus reveals the internal elements that can most effectively be targeted by sustainability interventions. See Appendix A for full description of vulnerability profiles.

Figure 3 |

 3 Figure 3 | Flowchart illustrating the key steps of the approach proposed. Step 1: identify the key vulnerability driver(s) to address through social and ecological vulnerability profiles. Step 2: for each component, determine a set of potential interventions to reduce each component's driver(s) of vulnerability. Step 3: consider the vulnerability profile of the associated component and determine a portfolio of potential interventions that minimizes trade-offs and promotes co-benefits. Step 4: ensure the viability of the interventions portfolio by reviewing identified interventions in the local context (e.g. institutional, management capacity, past experience, community aspirations). This generic framework can be adapted to each context by identifying specific interventions falling into

Figure 4 |

 4 Figure 4 | Assessment of the coral reef fishery of Moorea, French Polynesia, using spatially-explicit profiles of social and ecological vulnerability (Step 1). Since households mostly depended on resource located on adjacent reefs for provision and cultural services associated with fishing, combinations of social and ecological vulnerability profiles were spatially linked. Insets highlight different

Table legends Table 1 |

 legends1 Typology of interventions to manage resource-user dependencies, and implications for

Table 2 |

 2 Application of the generic typology of eco-and socio-centric interventions (Table1

a fishery-specific application of the general approach).

Acknowledgements.

We thank CRIOBE Service d'Observation CORAIL, and ISPF for providing the ecological and social monitoring data, respectively. LT PhD grant was funded by Pierre and Marie Curie University (PDI-MSC grant). This project received financial support from Agence Nationale de la Recherche (ANR-14-CE03-0001-01).

Supporting Information.

Appendix A. Additional conceptual information on vulnerability.

Appendix B. Case study methods.

Appendix C. Additional results.


Dominika Viglaska
email: dominika.viglaska@univ-reims.fr

Michael Rey
email: michael.rey@univ-reims.fr

Andrei V Nikitin

Vladimir G Tyuterev

Sorensen Curvilinear

Szalay Normal

Wang & Carrington Curvilinear

Makarewicz Curvilinear

Al Pearson

Normal

Elvibrot-Tnum

Derivation of ρ-dependent coordinate transformations for nonrigid molecules in the Hougen-Bunker-Johns formalism

Keywords: Hougen-Bunker-Johns Hamiltonian, nonrigid molecules, large amplitude motions, ρ, ρ, normal coordinates, coordinate transformations, internal axis system

In this paper, we report a series of transformations for the construction of a Hamiltonian model for nonrigid polyatomic molecules in the framework of the Hougen-Bunker-Johns formalism (HBJ). This model is expressed in normal mode coordinates for small vibrations and in a specific coordinate ρ to describe the large amplitude motion. For the first time, a general procedure linking the "true" curvilinear coordinates to ρ is proposed, allowing to express the potential energy part in the same coordinate representation as the kinetic energy operator, whatever the number of atoms. A Lie group-based method is also proposed for the derivation of the reference configuration in the internal axis system. This work opens new perspectives for future high-resolution spectroscopic studies of nonrigid, medium-sized molecules using HBJ-type Hamiltonians. Illustrative examples and computation of vibrational energy levels on semirigid and nonrigid molecules are given to validate this method.

I. INTRODUCTION

The astronomical observations based on new generation of space telescopes, satellites and exploratory missions provide exciting information about molecular universe 1 . In this context, the analyses of astronomical spectra are vital for the determination of chemical abundances and temperature conditions. There is an increasing interest in the investigation of molecular systems detected in both cold (Titan 2 , interstellar media 3 ) and hot atmospheres (exoplanets, brown dwarfs, and cool stars [4][5][6][7] ), requiring the development of accurate theoretical models for the modelling, analysis and interpretation of spectra. One of the challenging problems in computational molecular spectroscopy 3 is the study of nonrigid "floppy" molecules characterized by one or several large amplitude motions (LAMs). The use of adapted vibrational coordinates is thus required to describe LAMs such as internal rotation (e.g. C 2 H 6 , HNO 3 , HONO) [8][9][10][11] , floppy bending (e.g. C 3 , CH 2 ) [12][13][14] or inversion (e.g. NH 3 ) 15,16 . Because the study of high-lying LAM vibrational states lies in the frontier between spectroscopy and chemical dynamics, floppy molecules have attracted much attention from both theoretical and experimental sides. As stated by Bačić and Light in a review paper 17 , "studying floppy molecules is challenging because it forces us to reexamine, modify, and often abandon many of the basic concepts, formalisms, and computational methods applicable to small amplitude vibrations". These past few decades, a large number of theoretical studies 10,[18][19][20][21][22][23][24][25][26][27][28] has been devoted to the investigations of LAMs using different kinds of vibrational coordinates implemented in various computer codes (see Section III.).

The present work essentially focuses on the approach proposed in 1970 by Hougen, Bunker and Johns (HBJ) 29 . The HBJ model can be considered as the nonrigid counterpart of the semirigid Watson formalism 30 which over the years became the standard toolbox in highresolution molecular spectroscopy [START_REF] Papousek | Molecular Vibrational-rotational Spectra: Theory and Applications of High Resolution Infrared, Microwave, and Raman Spectroscopy of Polyatomic Molecules[END_REF][START_REF] Aliev | Higher-order effects in the vibration-rotation spectra of semirigid molecules[END_REF] . The conceptual difference with the commonly-used curvilinear nonrigid models is this similarity with the Watson formalism, augmented by an efficient combination of rectilinear, normal-mode coordinates to describe the small vibrations with curvilinear ρ ("real") and ρ ("effective") LAM coordinates and including their couplings with the rotational motion.

Here, the geometrically-defined ρ coordinate is used to build the potential energy function while ρ is involved in the kinetic energy operator (KEO). Although the link between ρ and ρ must be clearly established in order to express both the KEO and PES in the same set of coordinates, the combination of Q with (ρ, ρ) is threefold. (i) It turns out quite natural because each all the vibration characterized by Q will be small along the reference configuration defined by ρ. (i) It allows fast analytical computation of the matrix elements in variational calculations. (iii) It could further serve to derive "polyad" spectroscopic effective Hamiltonians for high-resolution analyses of nonrigid, medium-sized molecules. So far, explicit relations between ρ and ρ were derived only for triatomic molecules by Jensen [START_REF] Jensen | [END_REF] .

A generalization to an arbitrary number of atoms through a series of transformations is the purpose of the present paper that will be structured as follows. The strategy of the proposed approach as well as its main potential advantages are summarized in Sec. II. An overview of relevant theoretical studies devoted to nonrigid molecules using or not the HBJ philosophy is given in Sec. III. In the next section, we will recall general definitions about LAM coordinates and propose a method for the construction of the ρ-dependent reference configuration in the internal axis system. The procedure for all coordinate transformations as well as the construction of the HBJ Hamiltonian will be presented in Sections V. and VI. As numerical validations, the computation of vibrational energy levels for some 3-and 4-atomic semirigid and nonrigid molecules will be given in Section VII. This is a necessary step before considering larger nonrigid systems like CH 3 OH, C 2 H 6 , etc. detected in many astrophysical environments (e.g. giant planets, Titan 2,34 ) and for which the computation of highly-excited ro-vibrational states is still challenging.

II. KEYS OF THE PROPOSED APPROACH

Based on a series of coordinate transformations, our strategy consists in converting the ab initio curvilinear potential V naturally expressed in the geometrically-defined ρ and internal (or symmetry) coordinates to a function of 3N -7 normal coordinates Q and of ρ. The construction scheme of HBJ-type Hamiltonians for the nuclear motions in molecules can be outlined as

H HBJ = T HBJ (ρ, Q; -id/dρ, -id/dQ; J x , J y , J z ) + [V (ρ, R) → V (ρ, Q)],
where J α are the components of the total angular momentum in the molecular frame and T HBJ is the original kinetic energy operator derived by Hougen et al. 29 Let us summarize the main potential advantages of an extension of this approach for polyatomic molecules:

• Contrary to most of the KEO expressions in curvilinear coordinates, the form and complexity of T HBJ remain the same whatever the number of atoms. Similarly to the Watson KEO for semirigid molecules, a unique and simple implementation of T HBJ in a computer code would allow one to tackle the problem in the general case of nonrigid polyatomic molecules.

• The use of the normal mode Q and effective ρ coordinates in T HBJ allows an optimal separation between the small and large amplitude nuclear motions. This is achieved by the Sayvetz condition 35 .

• A polynomial expansion in Q and ρ in appropriate elementary functions for both KEO and PES may lead to an efficient algebraic solution of Schrödinger equation with the HBJ Hamiltonian, making thus possible the use of our reduction-compression tools [36][37][38] initially designed for variational nuclear motion calculations of semirigid molecules. A full account of symmetry via the use of irreducible tensor operators 36,39,40 can be also implemented.

• The algebraic formulation of the HBJ model could allow us to derive an ab initio effective model by applying a series of contact transformations 41,42 . In this context, an empirical optimisation of only a small set of Hamiltonian parameters will be sufficient to match high-resolution observed spectra 43 .

We expect that a systematic development and computer implementation of such an approach could pave the way for the computation of highly-excited states involved in the construction of complete lists of rovibrational transitions in high-temperature conditions. This is in spirit of the ongoing work on the construction of theoretical spectroscopic databases [44][45][START_REF] Tennyson | Proc. Int. Astron. Union[END_REF] for new molecular species relevant for astrophysical applications.

Methods based on the HBJ approach

The study of nonrigid systems is a long-standing topic which started in 1939 with the work of Sayvetz 35 . He was the first to introduce the notion of "anomalous" molecules with one or several vibrational LAM requiring a specific treatment that goes beyond conventional "semirigid" approach. However, no advanced calculations were carried out until the pioneering work by Hougen, Bunker and Johns 29 in 1970 who derived a prototype Hamiltonian and applied it mostly to triatomic molecules. In their approach, the LAM is treated as part of the rotational problem using a specific coordinate ρ. A so-called reference configuration playing the role that has the equilibrium configuration in a standard scheme is introduced.

It is chosen to minimize the coupling between LAMs and the small vibrations, achieved by the Sayvetz condition 35 . These past five decades, the HBJ model inspired a substantial number of studies aimed at the derivation of effective Hamiltonians using perturbation theory. A part of them can be divided into 3 classes, namely: the rigid bender (1970-1972), the semirigid bender (1977-1983) and the nonrigid bender (1982 and later) models. The first two ones are simple models based on some crude approximations while the development of the more elaborate nonrigid bender Hamiltonian was motivated by the needs of interpretation of experimental CH 2 spectra in its ground electronic state. The rigid bender Hamiltonian was used in the original HBJ paper 29 for the calculation of HCN and DCN vibrational levels and later by Bunker and Stone [START_REF] Bunker | [END_REF]48 for studying H 2 O, HCN and their deuterated isotopologues. An extension to the four-atomic ammonia molecule was reported out by Špirko and co-workers [49][50][51][52][53][54][55] in nonrigid invertor model. See table I for more detailed bibliographic references. An extension of the Contact Transformation method for nonrigid water and ammonia-type molecules has been considered in Refs. [56][57][START_REF] Starikov | Intramolecular ro-vibrational interactions and theoretical methods in the spectroscopy of non-rigid molecules[END_REF] . More recently, the BELGI computer code has been developed by Kleiner et al. [59][60][61][62][63][64] to study inversion or internal rotation of methyl groups and to compute both line positions and line intensities of nonrigid molecules.

Variational calculations, taking advantage of increasing computational resources, were performed for 3-atomic molecules [65][66][67][68][69][70][71] using the MORBID computer program developed by Jensen 21 . More recently, Yurchenko et al. have developed the TROVE code 20 based on the Sorensen's approach 72 and designed for general polyatomic molecules. The use of TROVE led to the construction of extensive ammonia line lists 73,74 from accurate ab initio potential energy surfaces 15,20,[75][76][77] .

Contrary to the computer codes like MORBID or TROVE employing curvilinear internal coordinates, even for the description of the small vibrations, we adopt here the original HBJ formulation based on the use of both rectilinear and curvilinear coordinates. This makes possible an algebraic formulation of the HBJ Hamiltonian for each point of the numerical grid in ρ.

Other types of approach

There exist other approaches to compute energy levels of nonrigid molecules without using the HBJ "philosophy". As stated above, the HBJ approach allows a maximum separation between the different types of molecular motions but in turn requires to truncate power series expansion in Q. This can be avoided by using exact KEO built in terms of curvilinear coordinates that are a natural choice for the description of floppy systems, combined with variational calculations (see for example Carter, Handy 78,79 , Schwenke, Huang et al. 18,[80][81][82][83] , Csaszar et al. 19,22,23,84 , Lauvergnat et al. 9,24,85 , studies in the framework of MCTDH 10,11,25,[86][87][88][89] , Carrington, Wang [90][91][92] , etc.). In such numerical calculations, no specific representation as e.g. Taylor expansion of the actual PES is required. A review of some earlier works can be found in Refs. [START_REF] Laane | Structures and Conformations of Non-Rigid Molecules[END_REF][START_REF] Smeyers | Structure and Dynamics of Non-Rigid Molecular Systems[END_REF][START_REF] Meinander | [END_REF][96] . Ab initio effective rotational and rovibrational Hamiltonians for nonrigid systems via curvilinear second order vibrational Møller-Plesset perturbation theory was considered in Ref. 97 .

To summarize, we give in Tab. I a non-exhaustive list of studies devoted to the treatment of nonrigid molecules, using or not the HBJ formalism.

IV. DEFINITIONS OF THE LAM COORDINATES-REFERENCE CONFIGURATION

As already mentioned, the HBJ approach inspired a series of works (see Tab. I). Their common element is the use of two classes of coordinates: large amplitude coordinates ρ i and normal mode coordinates Q k for the description of the small vibrations. The major difference a Extension of the TENSOR computer code initially designed for semirigid molecules 37,39 between the semirigid and HBJ approaches lies in the introduction of a so-called reference configuration a(ρ i ) following LAMs instead of the fixed, equilibrium geometry a e . If we consider a molecule with L large amplitude motions, the vibrational part will be described as follows:

Small amplitude motions

↔ (3N -6 -L) -like normal coordinates Q k Large amplitude motions ↔ ρ i , ρi with i = 1, • • • , L
In this paper, we will limit ourself to only one LAM (L = 1) associated with 2 specific coordinates, namely ρ 1 ≡ ρ and ρ1 ≡ ρ but the strategy would be similar for L > 1. As stated by Szalay 143 , ρ can be considered as a kind of "effective" mass-dependent coordinate for which all other coordinates are kept fixed. The HBJ kinetic energy operator depends explicitly on this variable while the potential function is built from the geometrically defined angle ρ (i.e. bond angles, torsional/dihedral angles, inversion angle) that will differ from ρ beyond the reference configuration (see Fig.

and we can easily notice that ρ and ρ coincide at the reference configuration.

A. Determination of the reference configuration in the Internal Axis System

By definition, the reference configuration depends explicitly on ρ and is defined with respect to a specific molecular axis system. It is common to use the Internal Axis System (IAS) for which the internal angular momentum due to ρ vanishes and where the zero-order cross terms

I (0)
αρ (α = x, y, z) of the 4×4 generalized inertia matrix I are zero. This condition of vanishing implies

I Rot-LAM = N i=1 m i U(ρ)a i × dU(ρ)a i dρ = 0, (2)
where N denotes the number of atoms. U(ρ) is a LAM-dependent rotation matrix that will transform the position vectors a i defined in an initial molecular axis system (MF 1 ) to those in the IAS for each value of ρ. In the standard HBJ derivation, the initial reference configuration was defined in the literature for 3 and 4-atomic nonrigid molecules as a function of an auxiliary angle (ρ) [START_REF] Jensen | [END_REF]116 that must satisfy Eq. ( 2). This angle is generally defined as the angle between the z-axis and one of the equilibrium bond length (see e.g. Fig. 1 of Ref. [START_REF] Jensen | [END_REF] or Fig. 5 of Ref. 144 ). Szalay and Ortigoso 145 proposed a way to derive the rotation matrix from the Floquet theory 146 . Recently, an alternative method to determine U(ρ) was considered 147 . It consists in integrating the kinematic equation

dU(ρ) dρ = U(ρ)Ω(ρ), (3)
using methods by Lie group and Lie algebra where Ω(ρ) is a skew-symmetric matrix involving the LAM dependent angular velocity vector components (ω x , ω y , ω z ) defined by ω(ρ) ≡

ω exact = -I Rot-Rot I Rot-LAM with I Rot-Rot (ρ)
dP dρ = 1 2 (p × ω) • S + ω • S 2 p cot p 2 - (p • ω)(p • S) 2 p 2 p cot p 2 -2 , ( 4
)
where p is the Euclidean norm. From the initial reference configuration of asymmetric ABC molecules (see Fig. 2 of Ref. 120 ), the condition (2) for HDO is not fulfilled. For NH 2 D (see Fig. 5 of Ref. 144 for the definition of MF 1 ), the condition (2) is not satisfied as well.

For example, we obtain the following rotation matrices (evaluated here at ρ = π/4) after solving (4)

U HDO (π/4) =      1 0 0 0 0.9607 0.2776 0 -0.2776 0.9607      , U NH 2 D (π/4) =      0.9792 0 0.2030 0 1 0 -0.2030 0 0.9792      . ( 5
)
which are nothing but rotations of angle (π/4) (see Eqs. (2.7) and (3.28) of Refs. [START_REF] Jensen | [END_REF]144 ) about the x and y axes, respectively. Note that no rotation is required for H 2 O and NH 3 when using the axis system of Ref. 120,144 . Another interesting example is nitric acid (HNO 3 ) possessing one internal rotation where ω(ρ) changes of direction (see Fig. 2)O(1) (see Fig. 3). As

Our strategy for the derivation of the ρ-dependent reference configuration in the IAS can be summarized as follows: (i) determination of the angular velocity vector ω(ρ) in the MF 1 , (ii) integration of Eq. ( 4) leading to the determination of the (p x , p y , p z ) components, (iii) calculation of the ρ-dependent rotation matrix U and (iv) computation of a

IAS i (ρ) = U(ρ)a MF 1 i (ρ).
Convergence of rovibrational energy levels computed in both IAS and MF 1 will be discussed in Section VII for the asymmetric HDS molecule.

V. CONSTRUCTION OF THE HBJ HAMILTONIAN: PHILOSOPHY OF THE PRESENT APPROACH

We describe in this part the strategy to build the HBJ Hamiltonian in a convenient form in order to transpose our tools of variational calculations initially designed for the treatment of rigid molecules [36][37][38] . The reference configuration a(ρ) defined in the IAS as well as the relation linking the coordinates ρ and ρ are key elements for the success of the proposed approach. A general procedure for the derivation of the relation ρ = f (ρ) will be proposed in the next section but for the time being it will be assumed known. Note that unless specified otherwise, all running indices i, j, k, rσr where Γ r is the irreducible representation for a mode r and σ r stands for a component to distinguish degenerate vibrations.

• • • will take values 1, • • • , 3N -7,
As a starting point, we assume that the potential function can be developed in power series of small-amplitude vibrational displacements

V = V 0 (ū(ρ)) + i Fi (ū(ρ))S i + 1 2 ij Fij (ū(ρ))S i S j + 1 6 ijk Fijk (ρ)S i S j S k + • • • , (7)
where Fi , Fij , Fijk are ρ-dependent linear, quadratic and cubic force constants. The first term in Eq. ( 7)

V 0 (ρ) ≡ V 0 (ū(ρ)) = r g r [ū(ρ)] r , (8)
is the 1D LAM potential defined at the reference configuration (V 0 (ρ) = V 0 (ρ) if all S i = 0). The notation ū stands for appropriate elementary functions of LAM ρ coordinate i.e.

ū(ρ) = ρ -ρe , sin ρe -sinρ, etc. which will be used for the modelling of the potential energy surface. In order to properly solve the GF problem, it is quite straightforward to show from Eqs. ( 7) and ( 8) that ū(ρ) takes necessarily the form

ū(ρ) = u(ρ) + i C i (u(ρ))S i + ij C ij (u(ρ))S i S j (9)
where the ρ-dependent C i (u(ρ)) and C ij (u(ρ)) coefficients are to be determined (see Section VI) and u is not necessarily the same function as ū. We can note that only terms of degree 2 in S in ( 9) are required to fully determine a set of normal coordinates. This is in contract with the MORBID approach where higher-order coefficients C ijk , C ijkl , etc., are needed 21 and have been determined analytically for triatomic molecules up to the fourth order (see Tab. 1 of Ref. [START_REF] Jensen | [END_REF] ). Assuming the C i and C ij coefficients known in this section, we can substitute Eq. ( 9) into ( 7) and solve the eigenequation

GF(ρ m )L(ρ m ) = λ(ρ m )L(ρ m ) ( 10
)
to finally obtain the eigenvalues λ i (ρ m ) (≡ ω i (ρ m )) as well as the elements L ij (ρ m ) of the corresponding eigenfunctions. This is done in the same fashion as in the Born-Oppenheimer approximation (separation of "slow" and "fast" motions) where the electronic equation is solved on a grid of nuclear geometries. The elements of the 3N ×3N -7 orthogonal transformation l(ρ m ) linking the mass-weighted Cartesian displacements to the normal coordinates are simultaneously computed using the same grid of points.

At this stage, it is shown that both the normal-mode KEO and the potential parts of the HBJ Hamiltonian can be completely determined from l(ρ m ). In order to further express the Hamiltonian as a power series in the (Q k , u(ρ)) coordinates, our strategy consists in fitting the 3N × 3N -7 elements l sα,k (ρ m ) to nth order polynomial functions as

l sα,k (ρ) = n r=0 lsα,k r [u(ρ)] r , (s = 1, • • • , N ). ( 11
)
The Coriolis parameters can be either computed directly from l(ρ) or empirically-adjusted from the same grid. We thus write

ζ Θ ij (ρ) = n r=0 ζΘ ij,r [u(ρ)] r , (Θ = x , y , z , ρ). (12)
Generally, n=8, 10 or 12 in Eqs. (11) and ( 12) is sufficient to obtain reasonably accurate rootmean-squares fit. This fitting procedure, though numerical, will make possible analytical computation of the matrix elements.

B. Transformation of the potential energy surface

In order to build an extended HBJ Hamiltonian, the initial PES (7) needs to be conveniently transformed to the normal mode and LAM coordinates. Unfortunately, the use of Eq. ( 9) is only valid for the solution of the GF problem within the quadratic approximation.

In this context, and by analogy with the Hoy, Mills and Strey formalism 151 , all quantities are expressed in Cartesian displacements making possible the use of the orthogonal transformation l(ρ) ≡ l(u(ρ)) defined in Eq. (11). We finally write

ū(ρ) = f (S, u(ρ)) Cartesian ------→ ū(ρ) = f (S(d), u(ρ)) d=M -1 2 l(u(ρ))Q --------→ ū(ρ) = g(Q, u(ρ)). ( 13
)
From this expression, we can easily convert the potential function as :

V (S, ū(ρ)) Cartesian, g(Q,u(ρ)) ----------→ V (S(d), g(u(ρ), Q)) d=M -1 2 l(u(ρ))Q --------→ V (Q, u(ρ)). (14)
The passage ( 14) from the ab initio potential V (S, ū(ρ)) to V (Q, u(ρ)) can be made from two different manners. The first way to proceed is to perform a Taylor-series expansion of all quantities around Q = u = 0. However, this procedure may suffer from a lack of convergence and thus fail to describe accurately the LAM if the form of the elementary function u is not appropriately chosen. Another route consists in performing this transformation numerically by generating a grid of points (Q h , ρ m ) from a(ρ m ) and l(ρ m ) and fitting a potential function in the full 3N -6-dimensional nuclear space. These two approaches will be discussed and compared in Section VI.

C. Formulation of the Hougen-Bunker-Johns Hamiltonian

Though only applied to tri-and tetra-atomic molecules 33,49,50 so far, a general and compact expression of the KEO has been proposed by Hougen, Bunker and Johns 29 in 1970. It is given by

T = 1 2 µ 1/4
α,β=x,y,z,ρ

(J α -π α )µ αβ µ -1/2 (J β -π β )µ 1/4 + 1 2 µ 1/4 3N -7 k=1 P k µ -1/2 P k µ 1/4 , (15)
where P k = -i ∂/∂Q k is the conjugate momentum of the normal coordinates Q k and J α (α = x, y, z) and J ρ = -i d/dρ are the components of the total angular momentum and LAM coordinate, respectively. All other terms are well-described elsewhere 29,[START_REF] Jensen | [END_REF]144 . It is important to note that according to Eqs. ( 11) and ( 12), all quantities involved in (15) will be ρ-dependent.

Finally, using Eqs. ( 14) and ( 15), the extended HBJ Hamiltonian can be written as

H(Q, P, ρ, J ρ , J α ) = T (Q, P, ρ, J ρ , J α ) + [V (S, ū(ρ)) → V (Q, u(ρ))]. (16)
Contrary to the previously published HBJ-based studies (see Tab. I) where a numerical treatment for the LAM coordinate was usually required (e.g. numerical integration using Numerov-Cooley), this Hamiltonian is written as a polynomial expansion in Q k and u(ρ) making thus possible the use of tensor operators for a full account of symmetry and computation of matrix elements from the Wigner-Eckart theorem. To this end, it is convenient to put the Hamiltonian ( 16) in a better suited normally-ordered form where each term is expressed as f (Q, ρ)g(P )J ρ w(J α ). Here f and g are not necessarily commuting functions and f (Q, ρ) does not commute with J ρ . In order to apply our reduction technique 39 for efficient variational calculations previously developed for semirigid molecules, creation and annihilation operators of vibrational modes could be also introduced.

VI. GENERAL PROCEDURE FOR THE DERIVATION OF ū(ρ) = f (S, u(ρ))

At this stage, the key relation (9) linking the effective and "real" LAM coordinates allowing the construction of the matrix F(ρ m ) on a grid, remains undetermined. We propose in this section a general method for computing the C i (ρ) and C ij (ρ) coefficients through a two-step procedure.

A. Transformations between quadratic forms in internal and Cartesian coordinates

The first step is the determination of a set of equations d = d(S, u(ρ)) relating the Cartesian displacements to the 3N -7 curvilinear and LAM coordinates. In other words, we have to invert the following set of nonlinear equations

S i = r B v i;rα (ρ)d rα + 1 2 rs B v i;rα,sβ (ρ)d rα d sβ + • • • (r, s = 1, • • • , N ), (17)
where the B v i;rα (ρ) ≡ B v i;rα (u(ρ)) and B v i;rα,sβ (ρ) ≡ B v i;rα,sβ (u(ρ)) coefficients depend on ρ via the reference configuration a(ρ). Bearing in mind that in the linear approximation, the B v (ρ)

T α = 1 √ M N i=1 m i d iα = 0, R α = 1 I ref αα (ρ) N i=1 m i αβγ a iβ (ρ)d iγ = 0, S = 1
I ref ρρ (ρ) N i=1 m i ∂a iα (ρ) ∂ρ d iα = 0. ( 18
)
Here, M is the total mass of system, corresponds to the Levi-Civita tensor and I ref αα , I ref ρρ are the diagonal elements of 4 × 4 inertia matrix 29 . Finally, the new matrix B lin (ρ m ) of dimension 3N ×3N defined such as (S, T α , R α , S) t = B lin (ρ m )d can be numerically computed on the grid of points and inverted. We thus write

d = B -1 lin (ρ m ) (S, T α = 0, R α = 0, S = 0) t . ( 19
)
However, quadratic terms S i S j are also involved in Eq. ( 9) making insufficient the use of Eq. ( 19). To properly compute C i (ρ) and C ij (ρ), we have to construct all the quadratic forms X.Y with X i = {S i , T β , R β , S} and Y j = {S j , T β , R β , S} from Eqs. ( 17) and ( 18), each of these forms being expanded at order 2 in Cartesian displacements to get a linear system of equations. This strategy closely follows the "linearisation" procedure proposed in Ref. 40 to invert the set of nonlinear equations S = S(Q). This leads to

(S, T α , R α , S, X.Y) t = B quad (ρ m ) (d, d.d) t , (20)
where the new matrix

B quad (ρ m ) is of dimension 9N (N +1) 2 × 9N (N +1) 2
. We have thus transformed the initial nonlinear problem ( 17) to a linear one. This transformation is schematically represented on Fig. 4 where B quad (ρ m ) is numerically evaluated for each value ρ m ∈

[ρ min , ρ max ]. For example, this matrix will be of dimension 54 × 54, 90 × 90 and 189 × 189 for 3-, 4-and 8-atomic molecules.

B. Computation of the linear Cartesian displacements and determination of

Eq. ( 9)

From Eq. ( 20), we are now able to express the linear and quadratic Cartesian displacements for each ρ m as a linear combination of X and X.Y as

(d, d.d) t = B -1 quad (ρ m ) (S, T α = 0, R α = 0, S = 0, X.Y ) t , (21)
with X.Y = 0 if X or Y equals T α , R α or S. At this stage, we are finally able to define the 3N linear Cartesian displacements by extracting 3N lines and (3N -7)(3N -4)

2 columns from B -1 quad (ρ m ) to give A inv quad (ρ m ) such that d = A inv quad (ρ m ) (S, S.S) t . ( 22
)
This transformation is represented schematically in Fig. 5 where the red and green colours highlight the extracted parts from B -1 quad (ρ m ). From Eq. ( 22), the link between ρ and ρ can be clearly established by means of geometry we directly obtain the numerical coefficients C i (ρ m ) and C ij (ρ m ) for the function ū(ρ) = cos(ρ). For pyramidal XY 3 molecules, ρ corresponds to the angle between the C 3 axis and the bond X-Y and can be defined by the relation

sin(ρ) = 2 √ 3 sin[(α 1 + α 2 + α 3 )/6)] ( 24
)
where α i is the angle between r j ≡NH j and r k ≡NH k with (ijk)=( 123), ( 213) and ( 312)

α i = cos -1 r j • r k r j r k , (25)
with r j = a j (ρ)+d j -a 4 (ρ)+d 4 . Similarly, we easily obtain the numerical coefficients C i (ρ m ) and C ij (ρ m ) for the function ū(ρ) = sin(ρ). As an illustrative example, Fig. 6 displays the behaviour of these coefficients for NH 3 as a function of u(ρ) = ∆ρ = ρ -π/2.

C. Comparison with existing expressions for 3-atomic molecules

In order to validate our approach, we considered analytical expressions available in the literature and derived by Jensen 33 for 3-atomic molecules. These latter were transformed to symmetry adapted coordinates and applied to H 2 S, SO 2 and O 3 152-154 . We have shown that our numerical coefficients C i and C ij evaluated at different values of ρ were all in excellent agreement (errors < 10 -8 ) with those obtained from the tabulated analytical expressions [START_REF] Jensen | [END_REF] .

- Another source of validation has consisted in comparing our results for the ozone PES 154 to those provided by P. Jensen from the MORBID code 21 . All the coefficients in front of ∆r 1 , ∆r 3 , ∆r 2 1 , ∆r 2 3 and ∆r 1 ∆r 3 obtained from both methods were identical, whatever the value of ρ.

VII. NUMERICAL VALIDATION: COMPUTATION OF VIBRATIONAL ENERGY LEVELS

The final validation is the computation of vibrational energy levels obtained from the extended HBJ Hamiltonian (16) and more generally from the procedure described in Sections V and VI. As a first step, we consider 3-atomic and 4-atomic semirigid molecules (H 2 S, SO 2 , O 3 , PH 3 ), for which numerically exact benchmark results are available using the well established Eckart-Watson formalism and potential energy surfaces from Refs. [152][153][154][155] . This permits checking all algebraic transformations described above with high precision at the

∆E = E W E -E HBJ (cm -1 ). H 2 S SO 2 O 3
Band, Sym E HBJ (cm -1 ) ∆E (cm -1 ) E HBJ (cm -1 ) ∆E (cm -1 ) E HBJ (cm -1 ) ∆E (cm -1 ) level of experimental accuracy. As the second step, we apply our formalism to the nonrigid ammonia molecule possessing large amplitude inversion motion and compare the calculations with the recent results in the literature using the potential energy surface of Refs. 15,77 . In the case of NH 3 , an updated and optimized version of the TENSOR computer code 37,156 initially designed for rigid systems, was used.

ν 2 (A
Hydrogen sulfide, sulfur dioxyde, ozone: For these three molecules, both the kinetic energy and the potential parts were expanded in Taylor series at order 18 in Q k and in the variable u(ρ) = ρ -ρ e . To ensure a convergence better than 0.001 cm -1 in variational cal- As illustrated in Fig. 7, the non-zero term MF 1 µ (0) xρ makes convergence of highly-excited rovibrational energy levels slower, in particular when using small basis sets. In that case, the rotational errors between IAS and MF 1 calculations may reach few wavenumbers. Though these errors are drastically reduced when increasing the number of basis function, we can see that the ro-vibrational energy levels are always better converged in IAS.

Phosphine: Though semirigid, this molecule is also a good candidate because it belongs to the non-Abelian point group C 3v and possesses two doubly degenerate irreducible repre-calculation is estimated to be better than 0.001 cm -1 . Again, the energy levels obtained from the Watson Hamiltonian were taken as reference benchmark values. All vibrational levels up to the octad are given in Tab. III showing the validity of our extended HBJ-type calculations. In Fig. 8, the strong deviation of the energy levels from the full calculation is displayed in case of simplified linear approximation (C ij = 0). This shows the necessity of including all terms in the expansion (9) proving once again that the ρ and ρ coordinates have essentially different nature.

Ammonia: Calculations on semirigid molecules partly proved that the methodology presented in Section VI was correct but a validation on a nonrigid system remains the ultimate proof of the concept. Using the refined NH 3 PES of Refs. 15 and 77 , the vibrational energy levels obtained from the computer codes GENIUSH 16,84 and TROVE 125 for nuclear motion calculations were taken as reference values. In order to properly describe the tunnelling splitting due to the inversion motion, appropriate basis functions had to be introduced in the TENSOR code for a rapid convergence. For the resolution of the GF problem (10), a grid of points ρ m has been generated in the range [0.75, 2.39] by step of π/1000. For both PESs, a reduced (Q k , u(ρ))-kinetic energy operator 16→10 was built using the the simplest linear u(ρ) = ρ -π/2 coordinate before converting to D 3h (M) tensor operators. Concerning the potential part, the S i and ū(ρ) = sinρ e -sinρ representation was considered before performing a 18th-order series expansion in Q k and u(ρ) to ensure a good reproducibility of the 1D double-well curve. The presence of spurious minima in such high order polynomial expansions makes the use of the reduction technique of Ref. 39 absolutely necessary.

After several trial tests, it turned out that only the 18→10 reduced potential had no holes, preventing from some spurious states. Unfortunately, this reduction does not allow to reproduce conveniently the vibrational energy levels for which the inversion mode ν 2 is involved.

We thus followed the strategy explained below Eq. ( 14) that consists in directly fitting the

V (Q k , ρ) potential.
To consider a lower order polynomial expansion, the variable u(ρ) ≡ u (A 1 ) = sinρ e -sinρ was preferred to u(ρ) ≡ u (A 2 ) = ρ -π/2. In order to fit a grid composed of 90000 points up to 7000 cm -1 generated from Eqs. ( 13) and ( 14), the maximum degrees involved in the D 3h (M) potential tensor operator of symmetry A 1 are 6, 6 and 8 in the (Q

(A 1 ) 1 , Q (E ) 3a,b ), u (A 1 ) (ρ) and Q (E )
4a,b coordinates, respectively. Using a nonlinear fitting procedure, we have obtained a rms error of 0.022 cm -1 . For the variational calculation, 6945, 5985 and 12915 functions of TABLE IV. Selected vibrational energy levels for NH 3 using the extended HBJ Hamiltonian model using the PES 15 and comparison with GENIUSH 16 . ∆ lin , ∆ T aylor and ∆ f it stand for the difference between GENIUSH and the linear approximation (see Eq. ( 9)), the Taylor potential and the fitted potential (see text), respectively. All values are in cm -1 . estimated to ∼0.001 cm -1 up to 5000 cm -1 .

Band

A comparison of the first NH 3 energy levels up to 5000 cm -1 between HBJ (this work) and GENIUSH 16 is given in Tab. IV. As expected, the linear approximation (C ij = 0 in Eq. ( 9), fourth column of Tab. IV) gives poor results. The 5th and 6th columns show deviations from GENIUSH using the Taylor and fitted potentials. As mentioned above, the biggest discrepancies arise for nν 2 -type levels when using the Taylor expansion while the agreement with GENIUSH becomes quite good using the potential function fitted on the grid; our precision being limited by the quality of the fit. Finally, we have considered the refined potential of Ref. 77 , but without the diagonal Born-Oppenheimer correction (DBOC)

for convenience of benchmark calculations, because the DBOC part have been expressed in a set of coordinates different from that of the Born-Oppenheimer part. In Tab. V, we compare our results obtained from the fitted potential to those given in Tab. 3.8 of Ref. 157 .

VIII. CONCLUSION

In this paper, we have proposed a general procedure for the series of transformations involving LAM coordinates required for the construction of the generalized Hougen-Bunker-Johns Hamiltonian. First, a combination of algebraic/numerical method is described to build an internal axis system through a ρ-dependent rotation matrix. Illustrative examples are given for HDO, NH 2 D and HNO 3 . Most importantly, we have developed a general procedure for expressing an ab initio curvilinear potential energy surface in the same coordinate system as the HBJ kinetic energy operator. Explicit relations relating the "true" and effective (ρ, ρ)

internal coordinates have been derived. In order to validate the proposed formalism, 3and 4-atomic semirigid and nonrigid molecules were considered. Good agreements with results obtained from the Watson Hamiltonian on the one side, and from the nuclear motion computer programs GENIUSH and TROVE, on the other side, confirm the validity of the proposed approach.

In the near future, an hybrid tensor/numerical approach will be considered to avoid fitting potential functions in the full nuclear space. This work can be thus considered as a first insight into future applications of more complex nonrigid species. This project aims at contributing to the construction of the theoretical spectroscopic databases of molecules possessing a strong astrophysical and planetological interest and to the extension of our TheoReTS Reims-Tomsk information system 44 to ethane-like molecules.

ACKNOWLEDGMENT

FIG. 1 .

 1 FIG. 1. Schematic representation of a 3-atomic nonlinear molecule in its reference configuration (defined by ρ, r e 21 and r e 23 coordinates) and in its instantaneous configuration (defined by ρ, r 21 , r 23 coordinates).

  1 in the case of triatomic molecules). As a simple illustration, Papoušek and Špirko 144 have determined for NH 3 the dependence of ρ on ρ when all instantaneous bond lengths and valence angles are equal, but different from the reference configuration. Here, ρ is the angle subtended by the NH bond of the reference configuration and the C 3 axis (see Fig. 4 of Ref. 144 ). In this particularly case, the following simple relation holds cos(ρ) = r r e cos(ρ),

  FIG. 2. Variation of I αρ (left) and ω α (right) as function of the LAM coordinate ρ for HNO 3 .

  2) and for which no analytical solution exists for the IAS rotation matrix. This nonrigid molecule exhibits a torsional motion (corresponding to the ν 9 band) of the O(3)-H bond (called top) with respect to the much heavier group (called frame) characterized by the ρ angle between the plane containing H-O(3) bond and another one containing O(3)NO(

FIG. 3 . 2 U HNO 3

 323 FIG. 3. HNO 3 in its reference configuration defined in the MF 1 and the definition of ρ.

  while Greek letters α, β, • • • will run over x, y, or z.A. Resolution of the GF problem-strategyOne of the major difficulties in the polyatomic extension of the HBJ approach is to transform the ab initio-based curvilinear potential part into normal and ρ coordinates in a systematic manner. To this end, we define a set of 3N -7 normal coordinates Q k by solving numerically the GF problem 150 on a grid of points ρ m ∈ [ρ min , ρ max ] where the values ρ min and ρ max will depend on the type of LAM. As already mentioned in Introduction, the main drawback of the HBJ approach lies in the fact that the kinematic G matrix, defined from the reference configuration a(ρ), is thus ρ-dependent while the force constant F matrix is determined from the potential function initially expressed in terms of 3N -7 internal coordinates and of ρ. In short, we have G = G(ρ) and F = F(ρ) making it necessary to find the transformation F(ρ) → F(ρ) to express both the potential and KEO parts in the same set of coordinates, namely of (Q k , ρ). In other words, this requires deriving the relation ρ = f (ρ, S) to compute the elements F ij (ρ m ) ≡ ∂ 2 V (ρ m )/∂S i ∂S j where {S} is a set of symmetry-adapted curvilinear coordinates. In more details, these coordinates should be denoted as S (Γr)

FIG. 4 .

 4 FIG. 4. Schematic representation of the B quad matrix defined in Eq. (20) for each value ρ m ∈ [ρ min , ρ max ].

  FIG. 5. Extracted part of B quad and schematic representation of A invquad defined in Eq. (22).

FIG. 7 .FIG. 8 .

 78 FIG. 7.Comparison of the rovibrational energy levels for HDS up to 4000 cm -1 when the reference configuration a(ρ) is computed in MF 1 (see text for the definition) and IAS and when using different basis sets. The HBJ Hamiltonian was expanded at order 14.

  Support from French ANR-16-CE31-0005 project, from Akademic D.Mendeleev program of Tomsk State University and from French-Russian collaborative LIA SAMIA program are acknowledged. The authors acknowledge support from the Romeo computer center of Reims-Champagne-Ardenne. We are grateful for valuable discussions to Prof. P. Jensen who has offered an access to the MORBID computer code during validation tests.

TABLE I .

 I Non-exhaustive list of contributions and computer codes devoted to nonrigid molecules.

	Contributor or computer code ‡	Coordinates Type	Ref.
	Rigid Bender	Normal †	E	29,47,48
	SemiRigid Bender	Normal †	E	48,98-104
	NonRigid Bender	Normal †	E	33,105-113
	NonRigid Invertor	Normal †	E	49-55,114-117
	Contact Transformations	Normal	E
	TENSOR ‡,a	Normal †	V	This work
	† Based on the HBJ or similar idea; E ≡ Effective Hamiltonian, V ≡ Variational calculations.

TABLE II .

 II Selected vibrational energy levels for H 2 S, SO 2 and O 3 computed from the HBJ model (E HBJ , this work) and comparison to those obtained from the Watson-Eckart Hamiltonian. Here,

TABLE V .

 V Selected vibrational energy levels for NH 3 using the extended HBJ Hamiltonian model using the PES 77 without the diagonal Born-Oppenheimer correction (DBOC) and comparison with TROVE given in Tab. 3.8 of Ref. 157 . All values are in cm -1 .

	Band (Sym.) This work (TW)	∆ (TW)-TROVE
	GS (A 2 )	0.775	0.000
	ν 2 (A 1 )	934.561	0.031
	ν 2 (A 2 )	969.566	0.011
	2ν 2 (A 1 )	1601.735	0.073
	ν 4 (E 1 )	1626.524	-0.031
	ν 4 (E 1 )	1627.596	-0.023
	2ν 2 (A 2 )	1884.248	-0.021
	3ν 2 (A 1 )	2386.461	-0.075
	ν 2 + ν 4 (E 1 ) 2542.976	0.062
	ν 2 + ν 4 (E 1 ) 2587.766	0.004
	3ν 2 (A 2 )	2897.189	-0.171
	2ν 2 + ν 4 (E 1 ) 3193.604	0.157
	2ν 4 (A 1 )	3216.397	-0.033
	2ν 4 (A 2 )	3217.976	-0.001
	2ν 4 (E 1 )	3240.644	-0.040
	2ν 4 (E 1 )	3242.015	-0.014
	rms (cm -1 )		0.068

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

sentations of type E. We have used the PES from Ref. 155 expressed as a power series in the S i and ū(ρ) = cosρ -cosρ e coordinates. The HBJ Hamiltonien ( 16) was developed at order 18 in Q k and in u(ρ) = ρ -ρ e before making the reduction at order 10 (see Refs.