Stéphane Richard-Devantoy

Massimiliano Orri

Josie-Anne Bertrand

Kyle T Greenway

Gustavo Turecki

David Gunnell

Chris Power

Ph.D Marie-Claude ; Geoffroy
email: marie-claude.geoffroy@mcgill.ca

Childhood cognitive skills trajectories and suicide by mid-adulthood: an investigation of the 1958 British Birth Cohort

Keywords: Suicide, Cognitive Abilities, Cognition, Academic Performance, Intelligence, Reading, Mathematics, Early-life influences, Longitudinal study, Birth cohort, Cohort Studies

Background: Poor cognitive abilities and low IQ are associated with an increased risk of suicide attempts and suicide mortality. However, knowledge of how this association develops across the life-course is limited. Our study aims to establish whether individuals who died by suicide by mid-adulthood are distinguishable by their child-to-adolescence cognitive trajectories.

Methods: Participants were from the 1958 British Birth Cohort and were assessed for academic performance at ages 7, 11, and 16 and intelligence at 11 years. Suicides occurring by September 2012 were identified from linked national death certificates. We compared mean mathematics and reading abilities and rate of change across 7 to 16 years for individuals who died by suicide versus those still alive, with and without adjustment for potential early-life confounding factors. Analyses were based on 14,505 participants.

Results: Fifty-five participants (48 males) had died by suicide by age 54 years. While males who died by suicide did not differ from participants still alive in reading scores at age 7 (effect size [g]=-0.04, p=0.759), their reading scores had a less steep improvement up to age 16 compared to other particpants. Adjustments for early-life confounding factors explained these differences. A similar pattern was observed for mathematics scores. There was no difference between individuals who died by suicide versus participants still alive on intelligence at 11 years.

Conclusions: While no differences in tests of academic performance and IQ were observed, individuals who died by suicide had a less steep improvement in reading abilities over time compared to same-age peers.

INTRODUCTION

There is strong and consistent evidence from large prospective studies of an inverse association between measures of cognition in childhood and adolescence (ages 7 to 18), [START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF][START_REF] Alaraisanen | Good school performance is a risk factor of suicide in psychoses: a 35-year follow up of the Northern Finland 1966 Birth Cohort[END_REF][START_REF] Osler | Impaired childhood development and suicidal behaviour in a cohort of Danish men born in 1953[END_REF][START_REF] Bjorkenstam | School grades, parental education and suicide--a national register-based cohort study[END_REF][START_REF] Andersson | Association of IQ scores and school achievement with suicide in a 40-year follow-up of a Swedish cohort[END_REF][START_REF] Gunnell | School performance and risk of suicide in early adulthood: follow-up of two national cohorts of Swedish schoolchildren[END_REF]) and/or early adulthood (ages 18 to 20) [START_REF] Gravseth | Suicide in young Norwegians in a life course perspective: population-based cohort study[END_REF]) and later suicide outcomes [START_REF] Alaraisanen | Good school performance is a risk factor of suicide in psychoses: a 35-year follow up of the Northern Finland 1966 Birth Cohort[END_REF][START_REF] Andersson | Association of IQ scores and school achievement with suicide in a 40-year follow-up of a Swedish cohort[END_REF][START_REF] Batty | Psychosis alters association between IQ and future risk of attempted suicide: cohort study of 1,109,475 Swedish men[END_REF][START_REF] Bjorkenstam | School grades, parental education and suicide--a national register-based cohort study[END_REF][START_REF] Gunnell | School performance and risk of suicide in early adulthood: follow-up of two national cohorts of Swedish schoolchildren[END_REF][START_REF] Gunnell | Low intelligence test scores in 18 year old men and risk of suicide: cohort study[END_REF][START_REF] Gravseth | Suicide in young Norwegians in a life course perspective: population-based cohort study[END_REF][START_REF] Osler | Impaired childhood development and suicidal behaviour in a cohort of Danish men born in 1953[END_REF][START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF]. In a large Norwegian cohort including 610,359 participants, poorer academic performance at 18 years predicted increased risk of suicide mortality up to three decades later [START_REF] Gravseth | Suicide in young Norwegians in a life course perspective: population-based cohort study[END_REF]. In another cohort following over 26,315 Swedish individuals an association between poor intellectual quotient (IQ) at 13 years and later suicide attempts in adulthood was fully explained by overall poor academic performance at age 16 years [START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF], suggesting that associations of poor IQ with suicide outcomes operate through its effect on educational attainment [START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF]. The few studies that reported sex differences suggested that an association of impaired cognition with suicide outcomes is seen in males, but less consistently in females [START_REF] Gunnell | School performance and risk of suicide in early adulthood: follow-up of two national cohorts of Swedish schoolchildren[END_REF]).

However, it remains unclear how cognitive deficits develop across the life-course, as most existing prospective studies have been conducted with measures of cognition at only one point in time. As far as we are aware, only one prior study examined associations between general indicators of cognitive development at different ages (e.g. mother reports of age at first speech, doctor report of overall alertness at age 7 and intelligence test at age 15 years) and later suicide mortality [START_REF] Neeleman | Predictors of suicide, accidental death, and premature natural death in a general-population birth cohort[END_REF]. Delayed speech development and low intelligence were risk factors for accidental death, but not suicide, though the number of suicides was low (n=11). To our knowledge, no existing study has documented the development of cognitive abilities in suicidal individuals with prospective repeated measurements of academic performance and IQ from childhood to adolescence (7, 11, 16 years) while accounting for a broad set of confounding factors. This is a key period in the life-course, as neuroimaging research suggests that intellectual capacity changes during adolescence [START_REF] Ramsden | Verbal and non-verbal intelligence changes in the teenage brain[END_REF].

In a previous study based on the 1958 British Birth Cohort, we failed to find evidence of associations between academic performance in reading and mathematics at age 7 and suicide mortality by age 50 years, but statistical power was limited due to the low number of suicides in the cohort [START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF]).

Here we investigate (a) whether associations between academic performance (reading and mathematics) and/or IQ and suicide mortality emerge with measures recorded later in childhood and adolescence and (b) to determine to what extent such associations (if any) are accounted for by pre-existing early-life potential confounding factors (i.e. child's birth order and weight, mother's age, father social class) known for their associations with both cognition [START_REF] Mosing | Associations between birth characteristics and age-related cognitive impairment and dementia: A registry-based cohort study[END_REF]) and suicide [START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF][START_REF] Orri | In-utero and perinatal influences on suicide risk: a systematic review and meta-analysis[END_REF] across the life-course. A better understanding of the development of cognitive impairments in suicidal individuals could provide important information about the etiology of suicide and its prevention.

METHODS

Participants

The 1958 British Birth Cohort is an ongoing longitudinal national birth cohort of 17, 638 individuals born in one week of March 1958 in England, Scotland, and Wales and representing 98% of all births in Britain in that week. Details about the cohort, including study design and response rates can be found elsewhere [START_REF] Power | Cohort profile: 1958 British birth cohort (National Child Development Study)[END_REF].

Participants were flagged for deaths on the National Health Service Central Register (NHSCR). The end of the mortality follow-up was September 2012 when participants were 54 years old. The NHSCR is not notified of deaths of emigrants; hence cohort members who had emigrated permanently from Britain (up to 2009) were excluded from our analyses. Information on academic performance and IQ was collected

Measures

Suicide mortality

As with our prior publication with the cohort, [START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF]) suicides were identified using the International Classification of Diseases, ninth revision (ICD-9) codes E950-59 (suicide) and E980-89 (undetermined intent) or tenth revision (ICD-10) codes X60-84 (suicide) and Y10-34 (undetermined intent). Suicide and death of undetermined intent were combined [START_REF] Gunnell | Time trends in coroners' use of different verdicts for possible suicides and their impact on officially reported incidence of suicide in England: 1990-2005[END_REF]. We have excluded pending verdicts (ICD-9 code 988.88; ICD-10 code Y33.9).

Cognitive ability from childhood to adolescence

At 7, 11 and 16 years, age-appropriate tests for mathematics and reading were administered by the participant's school teacher. The arithmetic test at age 7 comprised 10 problems with graded levels of difficulty; teachers read the questions to poor readers. At age 11, the mathematics test was constructed by the National Foundation for Educational Research in England and Wales. At 16 years, a mathematics comprehension test was constructed at Manchester University. The Southgate test [START_REF] Southgate | Southgate Group Reading Tests: Manual of Instructions[END_REF]) (range 0-30) was used to detect poor readers at age 7: children selected from several words one that corresponded to a picture, and teachers also read out words that the children identified from a list. Reading tests at ages 11 and 16 years were similar to the Watts Vernon comprehension test. At 11 years, a general ability standardized 80-item test approximating general intelligence (e.g. IQ) with verbal and non-verbal scales was also administered [START_REF] Douglas | The Home and the School[END_REF]. At each age, tests were standardized for month and year of assessment for each sex separately. To ease interpretation, we converted all scores to internally standardized Z scores (mean=0 and SD=1).

Early-life confounding factors

Potential confounding factors were identified from our prior studies with this cohort as risk factors associated with suicide deaths by 49 years (P ≤ 0.10) in multivariable analyses [START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF][START_REF] Orri | In-utero and perinatal influences on suicide risk: a systematic review and meta-analysis[END_REF]) and cognition [START_REF] Geoffroy | Child Neglect and Maltreatment and Childhood-to-Adulthood Cognition and Mental Health in a Prospective Birth Cohort[END_REF]). These include: low birthweight (low: <2.5kg vs normal: ≥2.5kg, recorded at birth); child's birth order (reported by mothers when their child was 7 years, including all live and still births and deaths by 7 years; coded 1, 2, 3 or ≥ 4), maternal age at the time of the study member's birth (categorized ≤19, 20-29 and >29 years) and father's social class in 1958 (using the 1951 Registrar General's Classification, categorized as non-manual (I/II/IIINM) and manual (IIIM/IV/V)). We have additionally controlled for fathers' social class at birth as it is strongly associated with cognitive skills [START_REF] Jefferis | Birth weight, childhood socioeconomic environment, and cognitive development in the 1958 British birth cohort study[END_REF].

Statistical analysis

First, using latent growth models we estimated trajectories of reading and mathematics abilities from age 7 to 16 years separately for two vital status groups, i.e., death by suicide versus still alive by September 2012. Models were estimated using maximum likelihood estimator in Mplus version 7.4. A major advantage of growth models is that they rely on the Full Information Maximum Likelihood estimation, which allows each participant with at least one data point to be included in the analysis [START_REF] Enders | The Relative Performance of Full Information Maximum Likelihood Estimation for Missing Data i. Structural Equation Models[END_REF]. Given that the distribution of reading abilities at age 7 was skewed, we performed sensitivity analyses using a robust (Huber-White) estimator for the standard errors, obtaining virtually identical results. Model fit was good according to the chi-square test of model fit (reading, chi-square 0.532; p=0.766; mathematics, chi-square=1.61, p=0.445), the Comparative Fit Indices (reading and mathematics, 1.00), and the Root Mean Square Error of Approximation (reading, 0.006 [90%CI 0.00-0.02]; mathematics, 0.00 [90%CI 0.00-0.03]). We used Wald tests to compare reading and mathematics abilities of individuals who died by suicide to those still alive at each time point (i.e., age 7, 11, and 16 years). Similarly, reading/mathematics rate of change (i.e. slopes) were compared using Wald tests to determine whether individuals who died by suicide versus those still alive differed in their rate of change. These differences were expressed as effect size (Hedge's g); values <0.20 are interpreted as "small", 0.21-0.50 "medium", 0.51-0.80 "large", and >0.80 "very large", similarly to Cohen's d. Given known sex differences in suicide mortality in this cohort [START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF]) and sex differences in the associations between cognitive and suicide outcomes [START_REF] Gunnell | School performance and risk of suicide in early adulthood: follow-up of two national cohorts of Swedish schoolchildren[END_REF], cognitive ability trajectories were estimated in the whole sample and separately for males and females (data shown for males only, as there were few suicides in females (n=7)), however findings must be interpreted with caution as there was no statistical evidence for sex-differences in the associations (P-values for mathematics =.841 and reading =.812).

Second, all models were re-estimated after adjusting for a priori identified potential early-life confounders.

We tested for non-linear associations for mathematics and reading skills at any time point, but none reached significance. In additional analyses, using t-tests we examined associations separately for intelligence and suicide mortality.

Missing data on covariables varied between 1% (father's social class) to 11% (birth order). To minimize further data loss, we imputed missing information on early-life confounders using multiple imputations by chained equations [START_REF] Azur | Multiple imputation by chained equations: what is it and how does it work[END_REF] and conducted analyses across the 20 imputed data sets (primary analysis). Results based on complete cases were comparable to those based on multiple imputations.

RESULTS

Population

Of 16,470 participants in the birth survey of 1958 and who had not emigrated permanently, we excluded 1524 participants who had died from causes other than suicide by age 54 (most deaths had occurred within the first year of life). A further 441 participants were excluded due to missing information on all mathematics and reading measurements at 7, 11, and 16 years, leaving 14,505 participants for analyses.

Of these, 55 participants had died from suicide before age 54 years (48 males and 7 females, including 12 of undetermined intent). The median age of suicide was 41 years (range 18-52 years) for males and 44 years (range 21-49 years) for females.

Reading and mathematics skills trajectories from age 7 to 16 years

Reading and mathematics skills trajectories during childhood for participants who died by suicide compared to those still alive are shown in Figure 1. For reading, participants who died by suicide had similar scores to those still alive at age 7 (effect size [g]=-0.02; p = 0.857) (model 1 in Table 1). At subsequent ages, the mean difference between the reading score of participants who died by suicide increased compared to those still alive (at age 11 years, g=0.09; p=0.555; at age 16 years: g =0.18; p=0.226). The effect size for the overall rate of change over time (slope) was moderate (g=0.21, p=0.099) indicating a less steep improvement in reading abilities for the participants who died by suicide.

Restricting the analyses to males yielded similar results. However, the mean slope difference was larger than that observed in the entire sample (g= 0.30, p=0.048 for males, g=0.10, p=0.830 for females).

For mathematics, a similar pattern was observed (model 1 in Table 2). However, the differences between participants who died by suicide and participants still alive were smaller than those observed for reading, including the slope difference (entire sample, g=0.01, p=0.438; males only, g=0.07, p=0.694).

Role of early-life factors in explaining reading skills trajectory differences between groups

We tested whether the observed difference in rate of change in reading between participants who died by suicide versus those still alive could be confounded by early-life factors. To do so, models were adjusted for key factors including paternal social class, maternal age at child's birth, birth weight, and birth order.

As shown in model 2 in Table 1, the estimated mean difference in reading between the slope of individuals who died by suicide versus those still alive reduced in both the entire sample (g=0.15, p=0.336) and in males only (g=0.18 p=0.483). The strongest associations with reading slopes were seen for birth order and maternal age. This suggests that the observed differences in the slope were explained by other factors in general, and by birth order and maternal age in particular.

We found no difference in intelligence scores at age 11 years between participants who died by suicide and those still alive (whole sample: suicide, mean= 0.00, SD=0.94; those still alive, mean=-0.10, SD=1.07; p=0.434; among males only: suicide, mean=0.01, SD=0.94; those still alive, mean=-0.09, SD=1.12; p = 0.508).

DISCUSSION

Our findings suggest that the participants who died by suicide had similar IQ scores and reading and mathematics scores in childhood and adolescence (as assessed with age-appropriate tests at school) than those still alive, but had a less steep improvement in their reading abilities. However, associations between reading ability trajectories and suicide mortality were accounted for by factors measured at birth that are associated with both suicide risk and cognitive skills (low birth weight, birth order, maternal age at childbirth, father's social class) [START_REF] Geoffroy | Child Neglect and Maltreatment and Childhood-to-Adulthood Cognition and Mental Health in a Prospective Birth Cohort[END_REF][START_REF] Geoffroy | Prenatal and childhood antecedents of suicide: 50-year follow-up of the 1958 British Birth Cohort study[END_REF].

Alterations (typically deficits) in various cognitive functions have been associated with increased suicide risk in various populations [START_REF] Keilp | Neuropsychological dysfunction in depressed suicide attempters[END_REF][START_REF] Keilp | Attention deficit in depressed suicide attempters[END_REF][START_REF] Richard-Devantoy | Deficit of cognitive inhibition in depressed elderly: a neurocognitive marker of suicidal risk[END_REF], Richard-Devantoy et al. 2014c[START_REF] Richard-Devantoy | Decision-making in unipolar or bipolar suicide attempters[END_REF], Richard-Devantoy et al. 2014a[START_REF] Vadini | Neurocognitive impairment and suicide risk among prison inmates[END_REF] and in various stages of life [START_REF] Keilp | Neuropsychological dysfunction in depressed suicide attempters[END_REF][START_REF] Keilp | Attention deficit in depressed suicide attempters[END_REF][START_REF] Richard-Devantoy | Deficit of cognitive inhibition in depressed elderly: a neurocognitive marker of suicidal risk[END_REF], Richard-Devantoy et al. 2014c[START_REF] Richard-Devantoy | Decision-making in unipolar or bipolar suicide attempters[END_REF], Richard-Devantoy et al. 2014a[START_REF] Vadini | Neurocognitive impairment and suicide risk among prison inmates[END_REF].

Our results are not consistent with previous studies showing negative associations between childhood cognitive abilities and suicidal risk such as performance (e.g., low grade point average) [START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF][START_REF] Bjorkenstam | School grades, parental education and suicide--a national register-based cohort study[END_REF][START_REF] Alaraisanen | Good school performance is a risk factor of suicide in psychoses: a 35-year follow up of the Northern Finland 1966 Birth Cohort[END_REF] or IQ [START_REF] Andersson | Association of IQ scores and school achievement with suicide in a 40-year follow-up of a Swedish cohort[END_REF][START_REF] Batty | Psychosis alters association between IQ and future risk of attempted suicide: cohort study of 1,109,475 Swedish men[END_REF], Gunnell et al. 2013[START_REF] Wallin | Suicide attempt predicted by academic performance and childhood IQ: a cohort study of 26 000 children[END_REF]), but the relatively small number of suicides (n=55) compared with 58-553 in previous studies limits power to detect relatively small effects. Our findings emphasise the importance of the developmental period from childhood to adolescence in understanding suicide [START_REF] Turecki | The neurodevelopmental origins of suicidal behavior[END_REF]). To our knowledge, our study is the first to show that differences in reading skills between suicidal and non-suicidal individuals are likely to emerge in the course of childhood and adolescence, rather than being evident in early childhood. Of note a prior study suggested that cognitive function at the age 18 years were more strongly associated to suicide mortality in adulthood than cognitive function at the age of 12 years, suggesting that the reduction in cognitive abilities per se may be a marker of suicidal risk [START_REF] Osler | Impaired childhood development and suicidal behaviour in a cohort of Danish men born in 1953[END_REF]).

In the current study, controlling for early-life influences abolished associations between change in reading abilities over time and suicide risk, suggesting that the observed differences could be due to the impact of pre-existing early-life factors.Early-life environment has also been found to be associated with longlasting cognitive impairments [START_REF] Skogen | The impact of early life factors on cognitive function in old age: The Hordaland Health Study (HUSK)[END_REF][START_REF] Zhang | The Long Arm of Childhood in China: Early-Life Conditions and Cognitive Function Among Middle-Aged and Older Adults[END_REF]. For example, low socioeconomic status in the family may result in decreased access to cognitively-stimulating environments for the child, which in turn might influence the development of his/her cognitive abilities [START_REF] Christensen | Socioeconomic status, child enrichment factors, and cognitive performance among preschool-age children: results from the Follow-Up of Growth and Development Experiences study[END_REF].

Although in our study, father's social class was measured at participants birth, social class is relatively stable over time (Power andMatthews 1997, Hackman et al. 2015), and it is possible that social class later in childhood rather than early-life per se explained the relatively weaker improvement in reading abilities over time. Moreover, young motherhood is associated with psychopathology, lower economic resources and maladaptive parenting practices, which are factors that can affect both suicidal risk and cognitive development of the individual [START_REF] Pinderhughes | Discipline responses: influences of parents' socioeconomic status, ethnicity, beliefs about parenting, stress, and cognitive-emotional processes[END_REF]. Such hypotheses should be formally tested in well-designed studies aiming to unravel the mechanisms linking early-life factors, development of cognitive skills, and suicidal risk. However, it is worth noting that in several prior studies reporting associations between IQ or school performance and suicide risk, these associations remained, even if attenuated, after controlling for early life variables such as socioeconomic status [START_REF] Gunnell | School performance and risk of suicide in early adulthood: follow-up of two national cohorts of Swedish schoolchildren[END_REF][START_REF] Bjorkenstam | School grades, parental education and suicide--a national register-based cohort study[END_REF][START_REF] Batty | Psychosis alters association between IQ and future risk of attempted suicide: cohort study of 1,109,475 Swedish men[END_REF][START_REF] Sorberg | Cognitive ability in early adulthood is associated with later suicide and suicide attempt: the role of risk factors over the life course[END_REF]).

In our study, the less steep improvement in academic abilities in individuals who died by suicide was of 'moderate' effect size in males and of 'small effect size and non-significant for mathematics. As far as we are aware our study is the first study to compare the development of reading and mathematics abilities in individuals who died by suicide versus those still alive, and further studies are needed to clarify the nature of associations between development of reading abilities in childhood and adolescence and later suicide mortality. For example, the less steep improvement in reading abilities observed among individuals who died by suicide in comparison to those still alive could be an early marker for poor executive functions which are in turn associated with suicidal risk [START_REF] Turecki | The neurodevelopmental origins of suicidal behavior[END_REF].

This study has a number of strengths including (1) the use of a large population-based sample followed from birth into adulthood (5 decades), (2) the use of validated tests for the assessment of reading and mathematics skills and objective (coroner) information to determine the cause of death, (3) the use of repeated measures of academic abilities allowing us to model the developmental course of reading and mathematics skills from early childhood to adolescence. Despite those strengths, some limitation must be acknowledged. First, the main limitation is the use of reading and mathematics skills measures as a proxy of cognitive performance or abilities. More specific cognitive measures such as tests of executive functions, memory, and decision-making would be more informative as previous studies have shown that deficits in suicidal patients specifically concerned these functions (Richard-Devantoy et al. 2014a, Richard-Devantoy et al. 2014b). For example, a recent meta-analysis including 25 studies revealed that patients with suicide attempt history perform significantly worse in the Iowa Gambling Task (i.e., a measure of decision making), the Stroop Test (i.e., a measure of executive functions), and a categorical verbal fluency test (Richard-Devantoy et al. 2014a). It would be of interest to have such measures incorporated in future longitudinal studies including data points from early life. However, there is a strong correlation between school performance and cognitive skills such as executive function, [START_REF] Fuhs | Longitudinal associations between executive functioning and academic skills across content areas[END_REF] and fine grained cognitive evaluation is more suited to small scale cross-sectional studies and are difficult to obtain in large representative population-based samples. Although comparison of our results with published literature is hampered by differences in measures (cognitive tests used in the 1958 British Birth Cohort (i.e. IQ and school performance tests) differ from tests used elsewhere (e.g., The Raven's Standard Progressive Matrices [START_REF] Alati | Is IQ in childhood associated with suicidal thoughts and attempts? Findings from the Mater University Study of Pregnancy and its outcomes[END_REF], and National Adult Reading Test (NART) [START_REF] Gunnell | Is low IQ associated with an increased risk of developing suicidal thoughts? A cohort study based on an 18-month follow-up of the national psychiatric morbidity survey[END_REF] we expect all tests to be correlated to some extent as recognized by general intelligence [START_REF] Deary | The neuroscience of human intelligence differences[END_REF].

A second limitation is that, despite the size of the cohort (n~14,000 participants), suicide is thankfully a relatively rare phenomenon, so we had limted statistical power to detect small effects. Further we were not able to examine the separate associations of cognitive skills with suicide mortality among females because of the small numbers of females who died by suicide in this cohort. Previous studies found that presence of psychosis reverses the association between IQ or school performance and suicide risk, i.e., whereby higher IQ or better school performance predict higher risk of suicide in individuals with psychosis [START_REF] Alaraisanen | Good school performance is a risk factor of suicide in psychoses: a 35-year follow up of the Northern Finland 1966 Birth Cohort[END_REF][START_REF] Andersson | Association of IQ scores and school achievement with suicide in a 40-year follow-up of a Swedish cohort[END_REF][START_REF] Batty | Psychosis alters association between IQ and future risk of attempted suicide: cohort study of 1,109,475 Swedish men[END_REF]. If some of our 55 individuals who died by suicide were affected by psychosis, this might have introduced heterogeneity in our sample and reduced the strength of our associations. However, given that 5% of individuals who died by suicide are affected by psychosis [START_REF] Hor | Suicide and schizophrenia: a systematic review of rates and risk factors[END_REF], it is unlikely that our estimates are significantly biased in this manner.

Conclusions

Our findings are consistent with the hypothesis that individuals who died by suicide had a poorer development in terms of their reading abilities in adolescence than those still alive, but no diffenrence in reading abilities was present in the early years of life. Future longitudinal studies are needed to document the development of cognitive abilities of suicidal and control individuals over time. Such a study could help determine if differences in cognitive performance between suicidal and non-suicidal individuals continue to increase over the lifespan.

  in childhood (7 and 11 years) and in adolescence (16 years) via standardized tests. Ethical approval was given (South-East Multicentre Research Ethics Committee ref. 01/01/44), and informed consent was obtained from all participants.

Fig. 1 .

 1 Fig. 1. Trajectories of mathematics and reading abilities from 7 to 16 years in individuals who died by suicide v. still alive, unadjusted associations.

This work was supported by the National Institute for Health Research Biomedical Research Centre at Great Ormond Street Hospital for Children NHS Foundation Trust and University College London. The views expressed in the publication are those of the authors and not necessarily those of the Department of Health. The funders had no input into study design; data collection, analysis, and interpretation; in the writing of the report; and in the decision to submit the article for publication. We thank the Centre for Longitudinal Studies (CLS), UCL Institute of Education, for the use of 1958-NCDS data and the UK Data Service for making them available. Neither CLS nor the UK Data Service bear any responsibility for the analysis or interpretation of these data. DG is supported by the NIHR Biomedical Research Centre at University Hospitals Bristol NHS Foundation Trust and the University of Bristol, England. The views expressed in this publication are those of the author(s) and not necessarily those of the NHS, the National Institute for Health Research or the Department of Health and Social Care. MO is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 793396 and the Canadian Institutes of Health Research. MCG is fellow of the Fonds de Recherche en Santé du Québec (FRSQ) and is supported by the American Foundation for Suicide Prevention (AFSP). Dr Turecki holds a Canada Research Chair (Tier 1) and a NARSAD Distinguished Investigator Award. He is supported by grants from the Canadian Institute of Health Research (CIHR) (FDN148374 and EGM141899). Drs Geoffroy and Turecki are supported by the Fonds de recherche du Québec -Santé (FRQS) through the Quebec Network on Suicide, Mood Disorders and Related Disorders.

Word count: 3340

Word count for abstract: 250

Character count for title: 105


Jean-Antoine Queyroux
email: queyroux@protonmail.ch

Sandra Ninet

Gunnar Weck

Gaston Garbarino

Mohamed Mezouar

Frédéric Datchi
email: frederic.datchi@sorbonne-universite.fr

Structure of liquid ammonia at high pressures and temperatures

Keywords:

 than with more recent equation of state (EoS) models. The experimental structure factor and PDF are well reproduced by AIMD simulations using either the BLYP or the PBE exchange-correlation functional. The shapes of S mol (Q) and g mol (r) vary little over the investigated pressure range and suggest a compact liquid with weak orientational correlations between molecules, which is corroborated by the coordination number varying from 12.7 to ∼14. The simulations are used to study the evolution of the site-site pair distribution functions, which reveals that the number of H-bonds per molecule (between 1.5 and 2) do not evolve with density, and that the distribution of H atoms around N atoms becomes more and more anisotropic with pressure.

I. INTRODUCTION

Ammonia is a major component of the giant icy planets and their moons and has many uses in the chemical industry, which explains the large amount of research devoted to its study. At the fundamental level, a longstanding interest in ammonia arises from its ability to form hydrogen (H) bonds. The H-bond refers to the interaction between a hydrogen atom and a negatively polarized atom (the acceptor A), such that a bridge, noted D-H• • • A, is formed between the atom to which the hydrogen is bonded (the donor D) and the atom A. It is well known that the properties of water are determined by the O-H• • • O H-bond. Without it the melting and boiling points of H 2 O would be more than 200 K lower and life as we know would be impossible. The hydrogen bonds are also responsible for the diversity of polymorphs known for water ice.

The H-bonds in NH 3 have been described as weaker and less directional than those formed by H 2 O, and their relevance to explain the structural motifs formed by NH 3 molecules is disputed. In the gas phase, NH 3 has been found to be a good H-bond acceptor but a poor donor [2]. In the solid phase I, stable below 195 K at ambient pressure, the N atoms adopt a pseudo face-centered cubic (fcc) arrangement and each molecule has 12 neighbors in the first shell [3]. Among the latter, 6 are at a N-N distance of 0.3375 nm and the 6 others are at 0.3940 nm. Those at the shorter distances have been described as H-bonded to the central molecule, 3 as acceptor and 3 as donor. Contrary to water ice H-bonds, these bonds are non-linear (the bond angle is 160 • ) and weaker than in water ice -a single lone electronic pair is shared by 3 H-bonds, and their existence was actually questioned by Ref. [2]. However both high-resolution x-ray diffraction [4] and ab initio molecular dynamics (AIMD) calculations [5] show that the electron density distribution along the bond is consistent with H-bonding. It is more difficult to state on the presence of H-bonds in the plastic solids II and III, which occur at higher pressures along the melting line, due to their large orientational disorder. Above 3.6 GPa at 293 K, ammonia transforms to another proton-ordered phase, phase IV, whose structure was determined by neutron diffraction [6]. In this phase, the N atoms have a pseudo hexagonal close-packed arrangement, and as in phase I, 6 H-bonds may be assigned to each molecule. But in contrast to phase I, the three N-D• • • H bonds are all different in length and angle. This asymmetry increases at higher pressure, in particular above the transition to the isosymmetric phase V, where the shortest and thus strongest H-bond becomes nearly linear at 178 • , while the two others remain bent at 150 • and 160 • [7]. In summary the presence of Hbonds in the ordered solid phases is well established, but unlike H 2 O, they do not have a strong impact on the solid structures -except probably for imposing proton ordering, since all of them are close to compact molecular arrangements.

The situation is less clear-cut in the liquid phase. Evidence for hydrogen bonding have been obtained from neutron diffraction [9][10][11][12], elastic [13,14] and inelastic x-ray scattering [15] experiments, as well as from theoretical simulations [16][17][18][19][20]. The N-N pair distribution function at 277 K obtained by the x-ray diffraction experiment of Ref. [14] presents a maximum at 0.34 nm and a shoulder at 0.37 nm, which, by analogy with the solid phase I, may be interpreted as the respective signature of H-bonded and non-bonded contacts, and thus suggests a well-structured first shell. More recent neutron diffraction experiments and simulations found, however, that the liquid is much less structured. They indicate that an ammonia molecule forms on average ∼2 H-bonds with its neighbors. These H-bonded molecules appear on the N-H(D) pair distribution function (PDF) as a shoulder at ∼0.25 nm of the strong oscillation peaked at ∼0.4 nm. Nevertheless, the spatial distribution of the nitrogen atoms shows that no extended H-bonded network exists in the ammonia liquid since the first solvation shell con-tains 12 to 13 molecules. Again, this is in sharp contrast with the case of liquid water where the solvation structure of a molecule is determined primarily by hydrogen bonding interactions. Except for the recent work of Guthrie et al. [12], all previous experimental studies of the structure of liquid ammonia have been carried out at or near ambient pressure. For water, several neutron [21,22] and x-ray diffraction experiments [23,24] have been reported to pressures of several GPa. They showed that the density increase between ambient and 1 GPa has a large effect on the structure, disrupting the 3D network of H-bonds although preserving the tetrahedral H-bonding scheme around a given molecule. In ammonia, the comparison between the structure at ambient and 1 GPa did not show a strong difference in the short-range structure, however a substantial increase of the spatial correlations at large distance was found [12].

The purpose of the present study is to extend the experimental investigation to higher pressures in order to determine the structural evolution over a larger range of densities. To this end, we used recently developed techniques for the measurement of the structure of low-Z fluids in the diamond anvil cell by synchrotron x-ray diffraction [25], which have been successfully used in the recent past to study liquid hydrogen [26] and carbon dioxide [27] and allowed to determine the melting lines of liquid nitrogen [28] and ammonia [29] up to the megabar range. Here we present experimental data for liquid ammonia obtained at 690 K and 800 K up to 6.3 GPa. The data is analyzed to extract both the x-ray molecular structure factor and pair distribution function. Experimental data are compared to ab initio calculations performed in the same pressure-temperature (P-T) range as experiments.

II. EXPERIMENTAL METHODS

A. Experiments

Ammonia samples (99.99%, Air Liquide) were loaded cryogenically into a membrane diamond anvil cell (DAC) equipped with diamond anvils of Boehler-Almax design [30], 300 µm culet size and 70 • 4θ x-ray aperture. A rhenium foil of 200 µm thickness was used as gasket and pre-indented to 80 µm. To prevent any chemical reaction, the ammonia sample was separated from rhenium by a gold liner of thickness about 5 µm. The latter was also used as the in situ pressure sensor, using the P-V-T equation of state (EoS) of gold from Ref. [31]. The pressures obtained with other EoS [START_REF] Jamieson | High-pressure research in Geophysics[END_REF][START_REF] Shim | [END_REF][34][35][36][37] deviate by less than 0.5 GPa in this pressure range. The DAC was heated with a ring-shaped external resistive heater enclosing the DAC and temperatures were determined with a type-K thermocouple (made of Chromel and Alumel alloys) in contact with the diamond anvil. The estimated pressure and temperature uncertainties are 0.15 GPa and 15 K, respectively.

Angular dispersive x-ray diffraction experiments were carried out at the ID27 beamline of the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The x-ray beam of wavelength 0.3738 Å (33 keV) was focussed to produce a spot size of 3×3 µm on the sample. A MAR345 image plate was used for detection. Diffraction patterns were integrated using the Dioptas software [38].

The main difficulty of present experiments is to extract the very weak coherent signal scattered by the sample from the incoherent (Compton) scattering originating from the 1.7 mm-thick diamond anvils. In order to reduce this background, a multichannel collimator (MCC) was used as described in Ref. [25]. The MCC is composed of two concentric sets of circular slits positioned between the sample and the image plate whose effect is to isolate a diffracting volume centered on the sample. The slits are carefully aligned to make their rotation center coincide with the sample. In order to cover all diffraction angles the MCC is rotated over the angular distance separating two adjacent slits during the collection. To further limit the parasitic noise coming from the experimental hutch, the MCC is enclosed by a lead cover extending up to the detector. The spatial selectivity of the slits increases with the diffraction angle and reduces the anvil contribution by a factor ∼100 at 20 nm -1 and ∼400 at 80 nm -1 . The disadvantage of the MCC is the loss of signal due to the absorption by the slits and to their limited height which only covers 16 • in the vertical direction. The acquisition time is thus increased relative to experiments without the MCC and was set here at 300 s. The background (empty cell) signal was measured at 300 K after pressure was released and the sample escaped from the gasket hole.

B. Data analysis

In order to process the measured intensity I meas (Q) (Q = 4π sin θ/λ is the momentum transfer) and obtain a molecular structure factor S mol (Q) and a molecular pair distribution function g mol (r), we use the methodology presented in details in Eggert et al. [39] and Weck et al. [25], which we briefly recall below.

The measured signal can be written as follows:

I meas (Q) = T DAC (Q)T M CC samp (Q)I samp (Q) + sI bkg (Q) (1) where T DAC (Q) is the DAC transmission, T M CC
samp (Q) is the MCC transmission for the sample scattering. They are calculated analytically using the known geometry of the DAC and MCC, as explained in Weck et al. [25]. I samp (Q) is the total scattering from the sample, I bkg (Q) is the background signal originating mainly from the Compton scattering of the diamond windows, and s is a scaling factor representing the ratio of the incident xray beam intensity used for the measurements of I bkg (Q) and I samp (Q). Since it is not possible to measure s with high enough accuracy, it is taken as a fit parameter in the analysis [39]. The background contribution was obtained by measuring the scattered signal from the empty cell at the end of the experiment. I samp (Q), can be expressed as the sum of the coherent and incoherent scattering, I coh (Q) and I incoh (Q), respectively,

I samp (Q) = 1 N α I coh (Q) + N p I incoh p (Q) ( 2
)
where α is the normalization factor used to express the sample scattering in atomic units and I incoh p (Q) is the incoherent Compton scattering from atoms of type p in the sample computed using the analytic atomic formulas given by Hajdu [40]. The molecular structure factor is defined as:

S mol ≡ I coh (Q) N Z 2 tot f 2 e (Q) (3)
Where f e (Q) is an effective electronic form factor defined as

f e (Q) = (f N (Q) + 3f H (Q))/Z tot , where f N (Q) is
the atomic form factor for N computed using the analytic atomic formulas given by Hajdu [40] while the atomic form factor for H is directly calculated from exact theoretical expressions [41].

S ∞ = (K 2 N + 3K 2 H )/Z 2 tot
, where K X (X = N, H), is an average effective atomic number defined as the average over the integration Q-range of the functions

K X (Q) = f X (Q)/f e (Q).
The molecular pair distribution function g mol (r) is then obtained by the inverse Fourier transform

F (r) ≡ 4πρ[g mol (r) -1] (4) = 2 π ∞ 0 Qi(Q) sin(Qr)dQ (5)
where i(Q) = S mol (Q)-S ∞ and ρ is the average molecular density.

The normalization coefficient α is calculated according the Krogh-Moe [42] and Norman [43] method:

α = Z 2 tot -2π 2 ρ + Qmax 0 [J(Q) + S ∞ ] Q 2 dQ Qmax 0 I samp (Q)/f 2 e (Q) Q 2 dQ ( 6
)
where

J(Q) = p I incoh p (Q)/Z 2 tot f 2 e (Q).
Initially, density ρ is not known and Eq. ( 6) is only exact when Q max = ∞. This makes the determination of α somewhat uncertain. We also need to know the scale factor s with high precision. To determine α and s we therefore use a self-consistent procedure based on a method developed by Kaplow et al. [44] and refined in Ref. [39]. The idea is to remove the nonphysical oscillations in F (r), dominant at low r, coming from errors in α and s. To this end, we use the fact that for r values below the minimum inter molecular distance (noted r min ), F (r) should match the theoretical form:

F theo (r) = -4πρ + F intra (r) for r < r min (7) where F intra (r) corresponds to the intramolecular contribution to F (r). A frozen molecule formulation is chosen for F intra (r), which considers the molecule as a rigid body with a fixed (N,H) distance of 1.1024 Å and (H,H) distance of 1.6242 Å [START_REF] Linstrom | NIST Chemistry webBook; NIST standard reference database No. 69[END_REF]. This approximation should be good enough in the pressure range of present experiments. The difference ∆F (r) between the initial F (r) from Eq. ( 4) and F theo (r) is calculated and Fourier transformed to evaluate the correction to i(Q). The new F (r) is then calculated using Eq. ( 4), and the procedure is iterated until convergence. As noted by Kaplow et al. [44] and Eggert et al. [39], ρ is an independent variable in this analysis with a direct influence on α and ∆F (r) so that ρ can be extracted by minimizing ∆F (r). The optimum value for the scale factor s can similarly be extracted. Thus ρ and s are varied to minimize ∆F (r) using a χ 2 figure of merit defined as χ 2 = rmin 0 [∆F (r)] 2 dr. During this optimization, a damping function of the form exp(-AQ 2 ) is used in the Fourier transform to reduce the termination effects from the finite value of Q max , but the final S mol (Q) and g mol (r) are calculated without the use of damping.

We also note that in order to determine the transmission of the sample scattering through the MCC, the knowledge of the sample thickness is required. This cannot be measured at present so it is taken as an additional fit parameter in the procedure [25].

Finally, in order to compare our S mol (Q) and g mol (r) with other experimental and theoretical data, we use the following relations:

S mol (Q) = (K 2 N S N N (Q) + 9K 2 H S HH (Q) + (8) 6K N K H S N H (Q))/Z 2 tot + S ∞ g mol (r) = K 2 N g N N (r) + 9K 2 H g HH (r) + 6K N K H g N H (r) Z 2 tot (9)
where S αβ and g αβ are the partial structure factors and partial pair distribution functions (PDF), respectively. For Q max = 90 nm -1 , K N = 8.4085 and K H = 0.5305, thus the respective contributions of g N N (r), g N H (r) and g HH (r) to g mol (r) are 70.7%, 26.8% and 2.5%.

C. Ab initio molecular dynamics

We conducted Born-Oppenheimer molecular dynamics simulations at 8 different densities using the PW code of the Quantum ESPRESSO package [START_REF] Giannozzi | [END_REF]. Most simulations used the Perdew-Burke-Ernzerhof (PBE) [47] exchange correlation functional with a plane wave energy cutoff of 800 eV, since PBE functional was found to provide very good agreement with experiment for the equation of state [48] and Raman modes [49] of solid ammonia at high pressure. We also compared the PBE functional results to those obtained with the Becke-Lee-Yang-Parr (BLYP) functional [50,51] at three different densities (29, 35 and 39.4 molecules.nm -3 ), as BLYP functional was found to agree well with the structural data of the ambient pressure liquid [17,18]. Brillouin zone was sampled at the Γ point only. The temperature was equilibrated with a Berendsen thermostat [52].

A supercell containing 108 NH 3 molecules was used. The system was equilibrated at two different temperatures (690 and 800 K) and each pressure was run in the NVT ensemble with a Berendsen thermostat during 10 ps.

III. RESULTS & DISCUSSION

A. Experiments

In total, seven x-ray diffraction patterns of fluid NH 3 were collected from 1.0 to 6.3 GPa at 690 K and 800 K. At these P -T condition, ammonia is in the supercritical state (the critical point of NH 3 is at 11.3 MPa, 405.5 K). To determine whether the data set crosses the Frenkel line delimiting the nonrigid (or gas-like) from the rigid (or liquid-like) states [53], we estimated the latter using the sound-velocity minimum criterion and the temperature dependence of the sound velocity from the Nationnal Institute of Standards and Technology database [START_REF] Linstrom | NIST Chemistry webBook; NIST standard reference database No. 69[END_REF], as in Trachenko et al. [54]. As can be seen in Fig. 1, the presure along the Frenkel line remains below 0.11 GPa for T < 800 K, thus all our data points correspond to liquid-like states.

Fig. 2 illustrates the raw signal obtained from a sample of fluid NH 3 at 4.3 GPa, 690 K obtained with and without the MCC. It can be observed that without the MCC, the signal of liquid ammonia cannot be distinguished from that of the empty cell above 30 nm -1 , while with the MCC, the signal stands above the background up to ∼70 nm -1 .

Nonetheless, the use of the MCC introduces parasitic oscillations due to mechanical imperfections of the MCC, which become more and more visible with increasing Q as the liquid signal becomes very weak. These oscillations are only partly removed when subtracting the empty cell contribution, thus in order to reduce their effects on the data treatment, such as artificial ripples in the pair distribution function, the sample scattering intensity is smoothed using a cubic spline smoothing routine which allows the amount of smoothing to vary as a function of Q [55].

The procedure outlined above to process the raw data requires to fix some parameters (Q max , r min , smoothing and damping factors), while the others are obtained from the optimization process (scale factor s and normalization factor α, sample thickness s th ).

One of the most critical parameter is the maximum Q value, Q max , to which the data is used. Although experimental data were collected up to 100 nm -1 , the rapid decrease of signal-to-noise ratio with Q imposed to set Q max to a lower value. We found here that the best compromise was obtained by setting Q max at 90 nm -1 , except for the lowest density data (1 GPa, 800 K) where Q max was set at 65 nm -1 . The optimum values of ρ and s obtained by the minimization of the low-r oscillations in F (r) are also sensitive to Q max . For our final ρ, we report the average values obtained by varying Q max in a range where S mol (Q) and g mol (r) are well behaved, and use the standard deviations as error bars [39]. The variation of the number density and scale factor when Q max is varied between 64 and 90 nm -1 is shown in Fig. 3 for the pattern at 4.8 GPa and 800 K. The average value at this pressure is 36.3 ± 1.5 molecules.nm -3 .

The number density obtained for the complete data set is plotted against pressure in Fig. 4. In the range of present experiment, the density continuously increases by 57% from 24.6 ± 1.1 to 38.7 ± 1.6 molecules.nm -3 . To our knowledge, there is no other reported experimental data in the literature in this P -T range. We thus compare our experiments to semi-empirical EoS, on one hand, and to the results of our AIMD simulations, on the other hand. The first semi-empirical EoS which we examined is that of Tillner-Roth et al. [1], since this EoS is nowadays considered as the reference EoS for liquid NH 3 . It is based on the parametrization of the free energy using a large set of experimental data for the thermodynamic properties of liquid ammonia up to 1 GPa and 720 K.

The second one is from the recent work of Jahangiri and Behnejad [56] and is based on a model developed by Song and Mason [58] using the statistical-mechanical perturbation theory of hard spheres. The latter EoS was found to fairly agree with the AIMD results of Bethkenhagen et al. [57] up to 20 GPa, 3000 K, and predicts a less compressible liquid than Tillner-Roth et al. [1]'s EoS, the difference between the two reaching 3 molecules.nm -3 (8 %) at 7 GPa. Our experimental densities are in very good agreement with the EoS of Tillner-Roth et al. [1] up to ∼ 4 GPa, then we find higher densities by up to 1 molecules.nm -3 at 6.3 GPa, but this difference remains within the experimental uncertainties. By contrast, the densities given by the EoS of Jahangiri and Behnejad [56] are too low and the difference with our data exceeds the error bars. This is consistent with the fact that the AIMD densities of Bethkenhagen et al. [57] also underestimates the experimental values. By contrast, our AIMD simulations using the PBE functional predict too high densities at low pressures but the agreement with experiment improves with pressure. This is consistent with the good agreement observed between DFT calculations of the EoS of solid NH 3 using the PBE functional [48]. We finally note that for a given density, the AIMD simulations with the BLYP functional give higher pressures than those using the PBE functional, and the resulting P -ρ data closely follows the EoS of Tillner-Roth et al. [1] in the covered P -T range.

The evolution of S mol (Q) with pressure is shown in Fig. 5(a). At all pressures, the molecular structure factor of fluid NH 3 is dominated by a first diffraction peak (FDP) located at 21 nm -1 at 1 GPa, 800 K, followed by broad oscillations. The main effects of increasing pressure is to sharpen the diffraction peaks, increase their intensity and shift them to higher Q. The maximum value of the FDP varies from 1.36 to 2.06 (+51%) and shifts in Q by 13.3% from 1 to 6.3 GPa. Fig. 5(a) also shows the S mol (Q) of liquid ammonia at 273 K close to the saturation pressure (ρ = 22.6 molecules.nm -3 ), which was calculated using FIG. 6. Evolution with density of the coordination number in liquid NH3. Squares and lozenges represent experimental and AIMD results, respectively. Red and blue are for T = 800 K and 690 K, respectively.

Eq. 8 and the partial structure factors reported by Ricci et al. [11]. There is a good correspondance between this S mol (Q) and the one we determined at 1 GPa, 800 K (ρ = 24.6 ± 1.1 molecules.nm -3 ), except that the latter has broader peaks, which is likely due to the effect of the higher temperature. The shape of S mol (Q) is at first glance, very similar to that of a liquid with no or weak orientational correlations such as H 2 [26,59], N 2 [60,61] and O 2 [61,62]. In dense monatomic liquids it is empirically found that the position of the maximum of the FDP, Q M is related to the density ρ through the relation Q M 4.4(4πρ/3) 1/3 [63], and thus (Q M /Q M 0 ) 3 ρ/ρ 0 where Q M 0 and ρ 0 correspond to a reference state. This relation has also been found to be verified in molecular liquids with weak orientational correlations, such as H 2 [26] and N 2 [28]. To test whether it also holds in ammonia, we plot in Fig. 5(c),(Q M /Q M 0 ) 3 as a function of ρ/ρ 0 , using the data of Ricci et al. [11] at 273 K, 0.483 MPa as the reference state. It can be seen that these two quantities are linearly related in liquid ammonia, but the slope of 0.76(1) is lower than 1, thus indicating that liquid NH 3 somewhat differs from a simple liquid. It can also be noticed that the Q 3 M data fall on the same line, irrespective of the temperature, and are thus primarily dependent on density.

Figure 5b shows the evolution of g mol (r) with pressure. The first peak in g mol (r) at ∼ 0.1 nm corresponds to the intramolecualr N-H peak and is fixed by our data analysis. Its broadening is entirely due to the truncation of the Fourier transform at Q max and can thus serve as a measure of the experimental resolution. The H-H intramolecular peak at 0.16 nm cannot be observed, most likely because it is too weak. Similarly as S mol (Q), the shape of g mol (r) is not much affected by pressure in the range of present experiments, the main effects being again a sharpening and intensification of the oscillations as they move to lower r. Consistently, we find a good correspondence between the g mol (r) of Ricci et al. [11], calculated using Eq. 9, with our lowest density data. A noticeable difference can however be observed in the lowr part which is more structured in Ricci et al. [11]'s data, which we believe is due to the lower resolution of our data (see also the comparison with AIMD simulations below). The present x-ray data also confirm the absence of the oscillations at 0.37 nm and 0.46 nm reported by Narten [14] which were not observed in the neutron diffraction studies [9][10][11][12].

The position of the first intermolecular peak approximately gives the distance between first neighbours, and is plotted versus density in Fig. 5(d). This distance varies from 0.345 to 0.315 nm over the experimental pressure range. The ratio between the positions of the first and the second peak are 2 and 1.8 for the lowest and the highest pressure reached in this experiment, respectively. These values characterize a liquid without extended preferential orientation, which is consistent with the observations made above.

Figure 6 shows the evolution of the average coordination number of ammonia molecules with density. This quantity was determined by integrating 4πr 2 ρg mol (r) up to the first minimum of g mol (r). Using our AIMD simulations, we checked that the difference in coordination number obtained by using g mol (r) instead of g N N (r) is below +6%. It is seen that the coordination number slowly increases with density, going from 12.7 to 14.1. These numbers are similar to those found earlier in the neutron diffraction studies (14 in Ricci et al. [11]) and indicate a close packing of molecules at all pressures.

Our experimental data thus suggest that the structure of liquid NH 3 is weakly affected by pressure. Unfortunately, the interpretation of the data in terms of structure is limited because g mol (r) is a weighted sum of three partial PDFs (Eq. 9) which cannot be extracted from a single diffraction pattern. In the following, we compare our experimental data to AIMD simulations and use the latter to obtain greater details on the liquid structure.

We first start by comparing the AIMD predictions for the partial PDFs using the PBE functional, on one hand, and the BLYP functional, on the other hand. This is illustrated in Fig. 7(a) for the highest density of our data set, ρ=39.4 molecules.nm -3 , and for T = 800 K. We find an almost exact correspondence between the two sets of partial PDFs. These two functionals thus give identical results in terms of liquid structure at a given density, however the predicted pressure is higher (8.6 GPa) for BLYP than for PBE functional (6.3 GPa). Next, we compare the experimental g mol (r) with the simulated ones, calculated with the help of Eq. 9. This is shown in Fig. 7(b) for four simulation runs at densities close to exper- imental ones. There is a fair agreement between the experimental and theoretical g mol (r) over the probed range of densities: the peak positions are well reproduced, as well as the amplitude of the second oscillation. The first peak is smaller and broader in the experiment, which we think is due, at least in part, to the low Q max truncation. Indeed the same effects were previously observed in the x-ray study of water [23,39]. A noticeable feature is that the non-zero intensity between ∼ 0.18 and ∼ 0.22 nm, which, as seen later on, is assigned to NH and HH contributions, is observed in both experimental and theoretical g mol (r). Due to the limited size of the simulation box, the theoretical g mol (r) only extends to 0.7 nm, which was found insufficient to determine meaningful S mol (Q) by inverse Fourier transform. This prevents a direct comparison with the experimental S mol (Q) which do not suffer from the "Q max problem". The good agreement between simulations and experiments also extends to the density evolution of the first-neighbor distance and coordination number, as seen in Fig. 5(d) and Fig. 6.

The present AIMD simulations enable to access the evolution of the partial (or site-site) PDF with density which are shown in Fig. 7(c). Like g mol (r), to which it contributes by ∼ 70%, g N N (r) presents a first peak which is unstructured and has a similar shape as encountered in simple compact liquids. There is, however, a noticeable asymmetry on the high-r side at low density, which tends to decrease as the peaks sharpen and move to lower r with pressure. More interesting is the evolution of the NH distribution. At low density, the latter displays a first intermolecular peak with little structure, but on increasing density, several humps appear and the peak seems to be composed of at least three different contributions at r 0.21, 0.29 and 0.37 nm. Turning to g HH (r), the first intermolecular peak presents a bimodal distribution at all densities, and the effect of increasing pressure, apart from the general shift, is to increase the intensity of the low-r contribution. These results suggest that, despite the compact nature of the first-neighbour solvation shell, the H distribution around nitrogen atoms are not isotropic and that orientational correlations exist.

We recall that the presence of H-bonded molecules in the low-pressure, low-temperature liquid was inferred by the presence of a pronounced shoulder in the neutron g N H (r) at r 0.235 nm [11]. In Fig. 8, we compare the experimental partial PDFs of NH 3 at 273 K, 0.483 MPa determined by Ricci et al. [11] to the AIMD simulations of Boese et al. [18] at the same thermodynamic conditions, on one hand, and to our own simulations at the lowest density (24.5 molecules.nm -3 ) at 700 K, on the other hand. The simulations of Boese et al. [18] used the BLYP functional and reproduced fairly well the experimental PDFs, except for slight difference in amplitudes. In particular, the shoulder at ∼ 0.235 nm is present in the NH distribution. Our simulations used the PBE functional, but as seen above, we found no difference between the PDFs produced with the PBE and BLYP functionals. It can be seen that at the higher temperature of our simulations, the shoulder is absent, but the total width of the first intermolecular peak is very similar. This indicates that H atoms which form H-bonds to the central molecules cannot be distinguished from those which do not, which implies that non-bonded molecules can move closer to the central one in the high temperature liquid. FIG. 8. Partial PDFs of liquid ammonia at low density. The green solid, blue dotted and red dashed lines are, respectively, the experimental data of Ricci et al. [11] at 273 K, 0.483 MPa, the AIMD simulations of Boese et al. [18] with the BLYP functional at the same P, T conditions, and our AIMD simulations at 24.5 molecules.nm -3 , 700 K using the PBE functional.

The integration of 4πρr 2 g N H (r) using the same boundaries as Ricci et al. [11], i.e from 0.19 to 0.26 nm, gives an average number of H bonds of ∼1.5 per molecule, instead of ∼2 at 273 K. We note however that the apparent difference with the low-temperature liquid should be taken cautiously as the number of H-bonds are very sensitive to the integration boundaries which are arbitrarily set since there is no resolved H-bond peak.

The shape of the NH distribution found at the highest density (39.4 molecules.nm -3 ) is sensibly different from the ones at low density, as well as from Ricci et al. [11]'s g N H (r) at 273 K. The low-r shoulder assigned to H-bonded H atoms is recovered, but remains less marked than in the liquid at 273 K. The integration of 4πρr 2 g N H (r) up to the inflection point at ∼0.235 nm gives an average of ∼2 H bonds per molecule at this density, which is the same as for the low-T liquid, showing that the number of H-bonds per molecule is not much sensitive to density.

IV. CONCLUSION

In this paper, we have reported the first, to our knowledge, structural data on liquid ammonia above 1 GPa and up to 6.3 GPa at 800 K, using synchrotron x-ray diffraction in a diamond anvil cell. In this P -T range, liquid ammonia is found stable and experience a density increase of 57 %. The present experiments made use of advanced filtering techniques enabling to reject a large part of the background scattering from the diamond anvils. This was found crucial in order to extract highquality data for the x-ray molecular structure factor up to 90 nm -1 , from which the x-ray molecular pair distribution function was obtained. The present study also enabled to estimate the density of liquid NH 3 in a range where no experimental data existed.

The shapes of S mol (Q) and g mol (r) were found to little vary with density in the experimental range, and are found similar to that of liquids with weak orientational correlations between molecules. The coordination number varies from 12.7 to ∼ 14, which is typical of compact liquids. In order to obtain more details of the structure, AIMD simulations using either the PBE or BLYP functional were performed. These two functionals are found to give identical results in terms of pair distribution functions at a given density, but BLYP functional gives higher pressure values than PBE functional. A good agreement is observed between the simulated and experimental g mol (r), first-neighbor distance and coordination number at the same density. The most significant effect of the density increase is seen on the NH distribution, which displays a intermolecular peak with no apparent structure at low density while three separate humps appear high density. The signature of H-bonds is more difficult to observe than for the liquid at low temperature, but the integration of the NH distribution over the range of expected H-bond distances gives on average from 1.5 to 2 H bonds, similar to the low T liquid.

The present data should be useful to benchmark theoretical models of ammonia in the 10 GPa, 1000 K range. Experimentally, the extension of the present work to higher pressures would be very interesting as theoretical studies [57,64] indicate major transformations in the liquid phase above ∼ 20 GPa, where molecules start to dissociate. Although x-ray diffraction only allows accessing the total pair distribution functions, such a radical change in the liquid should be easily observed.

FIG. 1 .

 1 FIG. 1. (Top) Phase diagram of NH3 after Ninet and Datchi [8]. Red circles show the P-T points where experimental data were taken, black star is the critical point and dashed line is the Frenkel line. (Bottom) Representation of the structures of NH3-I and NH3-IV. Colored lines show hydrogen bonds, with different colors for bonds of different lengths and angles.

1 )

 1 FIG. 2. Integrated x-ray diffraction patterns from a NH3 fluid sample at 4.3 GPa, 690 K, with and without the MCC. The red and blue solid lines are patterns obtained respectively with and without the MCC. The red and blue dashed lines are the respective patterns from the empty cell with and without MCC. The intensities for the patterns obtained without MCC have been divided by 20 to scale with those measured with the MCC.

FIG. 3 .

 3 FIG.3. Optimum density and scale factor for ammonia at 4.8 GPa and 800 K as a function of Qmax. The average density is 36.3 ± 1.5 molecules.nm -3 .

FIG. 5 .

 5 FIG.5. Experimental results from the present x-ray diffraction study of liquid NH3. Panels (a) and (b) show respectively the molecular structure factor and pair distribution functions obtained at the indicated pressures. The curves are offset by 0.5 for clarity. The temperature was 800 K except for the patterns labelled by the * symbol, for which temperature was 690 K. The dashed lines are S mol (Q) and g mol (r) calculated using the experimental data from Ricci et al.[11] at 273 K and 0.483 MPa. the Panel (c) shows (QM /QM0) 3 as a function of ρ/ρ0, where QM is the position of the main diffraction peak, ρ is the liquid density, QM0 = 20.65 nm -1 and ρ0 = 22.6 molecules.nm -3 correspond to the values of liquid ammonia at 273 K and 0.483 MPa reported by Ricci et al.[11]. Filled and open squares are experimental data from the present work and Ref.[11], respectively. Blue and red are for T = 690 K and T = 800 K, respectively. The solid line is a linear fit to the data. Panel (d) shows the first-neighbor distance as a function of density, determined as described in the text. Squares are present experimental data, and lozenges are from AIMD simulations with the PBE functional.

FIG. 7 .

 7 FIG.7. Results from AIMD simulations. (a) Comparison between the NN, NH and HH partial PDFs obtained using either the PBE or BLYP functional at the fixed density of 39.4 molecules.nm -3 and 800 K. Offsets of 1.5 and 3 are applied to NH and HH distributions respectively. The arrows designate the three humps at r 0.21, 0.29 and 0.37 nm in the NH partial PDF. (b) Comaparison between the experimental (red solid lines) and theoretical (blue dashed lines) g mol (r) at close densities. The ρ values in molecules/nm 3 are indicated for each curve. The theoretical g mol (r) were obtained with the PBE functional. Curves are offset by 1. (c) Evolution with density of the theoretical partial PDFs. The blue dotted, green dashed, orange dot-dashed and red solid lines correspond to ρ values of 24.5, 29, 33.8 and 39.4 molecules.nm -3 , respectively. Offsets of 1.5 and 2.5 are applied for NH and HH distributions, respectively.

ACKNOWLEDGMENTS

We acknowledge the ESRF for provision of beam time under Long-Term HD-463 and the Agence Nationale de la Recherche for financial support under grant ANR 13-BS04-0015 (MOFLEX).