Some of the Research highlights from our users in 2017:

2017 list of publications


New Models give insight into the heart of the Rosette Nebula

New research, led by the University of Leeds, offers an explanation for the discrepancy between the size and age of the Rosetta Nebula’s central cavity and that of its central stars.

The Rosette Nebula (Credit Nick Wright, Keele University)

The Rosette Nebula in the Milky Way Galaxy, roughly 5,000 light-years from Earth, is known for its rose-like shape and distinctive hole at its centre. The nebula is an interstellar cloud of dust, hydrogen, helium and other ionised gases with several massive stars found in a cluster at its heart.

Stellar winds and ionising radiation from these massive stars affect the shape of the giant molecular cloud. But the size and age of the cavity observed in the centre of Rosette Nebula is too small when compared to the age of its central stars — something that has has puzzled astronomers for decades.

Through computer simulations run in part on DiRAC Resources, astronomers at Leeds and at Keele University have found the formation of the Nebula is likely to be in a thin sheet-like molecular cloud rather than in a spherical or thick disc-like shape, as some photographs may suggest. A thin disc-like structure of the cloud focusing the stellar winds away from the cloud’s centre would account for the comparatively small size of the central cavity.

Study lead author, Dr Christopher Wareing, from the School of Physics and Astronomy said: “The massive stars that make up the Rosette Nebula’s central cluster are a few millions of years old and halfway through their lifecycle. For the length of time their stellar winds would have been flowing, you would expect a central cavity up to ten times bigger.”

A projection of the simulated molecular cloud and nebula along a specific line of sight (and) A slice through the simulation of the Rosette Nebula. (Credit: C. J. Wareing et al.)

“We simulated the stellar wind feedback and formation of the nebula in various molecular cloud models including a clumpy sphere, a thick filamentary disc and a thin disc, all created from the same low density initial atomic cloud.

“It was the thin disc that reproduced the physical appearance – cavity size, shape and magnetic field alignment — of the Nebula, at an age compatible with the central stars and their wind strengths.”

More information on this result can be found here, the STFC Press Release can be found here and Alt-metrics for the paper, indicating it has had a very high level of attention, can be found here here.


Gravitational Waves 2017

2017 saw a series of impressive discoveries during the LIGO-Virgo second observing run. Early in January we made our third binary-black-hole observation (GW170104). Another black-hole signal arrived later in June (GW170608).

In August came two impressive observations. The first followed the inclusion in August of the French-Italian Virgo detector in the observing run, with the first binary-black-hole observation that included three detectors, on August 14th. A third detector made it possible to accurately triangulate the position of the source on the sky, making a dramatic improvement over previous observations.

Only three days later, on August 17th, came the first observation of gravitational waves from two neutron stars spiralling together. This observation was astounding in several ways. The source was only 40 Mpc away, ten times closer than the breakthrough binary-black-hole merger observed in 2015; that lead to the strongest signal we have yet observed. This signal made it possible to provide the strongest support yet for Einstein’s general relativity, an independent measurement of the Hubble constant, the first strong-field constraints on the equation of state of neutron stars.

Beyond that, the results of the merger – a gamma-ray burst and afterglow – were observed by many EM telescope teams. This marks the beginning of multi-messenger astronomy, and the implications for astronomy and astrophysics are still being understood.

Fig. 1: Time-frequency map of two neutron stars spiralling together.

An integral part of all of these GW measurements were the theoretical models of binary-black-hole signals that were developed with the aid of numerical simulations performed on DiRAC. In the weeks following the first binary-neutron-star observation, these models were extended to include neutron-star tidal terms. In preparation for the next observing run, due to begin in early 2019, these models have also been extended to an approximate model of subdominant harmonics, and this model is now being more precisely tuned to numerical simulations.



Shining Light on the Structure of Hadrons

Understanding matter at its deepest level is complicated by the fact that the fundamental constituents known as quarks are hidden from us. Instead of being available for experimental scrutiny, as electrons are, quarks are always bound together into particles called hadrons. We have to infer the properties of quarks, and that of the strong interaction that binds them, from experimental studies of hadrons. A window into the hadron is provided by striking it with a photon, the quantum carrier of the electromagnetic interaction. The photon can penetrate the hadron to interact with the quarks inside, but the bound state nature of the hadron is revealed by the fact that the struck quark must redistribute the photon’s momentum between its fellow constituents for the hadron to change direction as momentum conservation requires (see Fig. 1). The hadron’s electromagnetic form factor, F(Q2), parameterises this bound-state behavior as a function of the squared 3-momentum transferred from the photon. New experiments at Jefferson Lab in the USA will shortly be taking data to determine the electromagnetic form factor of hadrons called pi and K mesons up to much larger values of Q2 than has been possible before. The aim is to reach a regime in Q2 where the strong interaction becomes relatively weak and F(Q2) is expected to have simple behavior from QCD, the theory of the strong interaction.

Fig 1. When a pi meson interacts with a photon (wavy blue line), momentum must be redistributed via gluon exchanged (curly blue line).

However, it is not at all clear where this regime begins. Hadron electromagnetic form factors can also be calculated using the numerical techniques of lattice QCD. The HPQCD collaboration has recently been able to reach a value of Q2 of 7 GeV2 for the determination of F(Q2) for a ‘pseudopion’, a particle like the pi meson but made of strange quarks instead of up and down quarks.

Fig. 2. Lattice QCD results from arXiv: 1701.04250 compared to a simple pole model at low Q2 and leading-order perturbative QCD at high Q2.

This improves numerical speed and statistical accuracy, but still allows us to compare to high-Q2 expectations. We find Q2F(Q2) to be flat in Q2 but in clear disagreement with the simple high Q2 result from perturbative QCD (see Fig. 2).

Calculations of F(Q2) for pi and K meson are underway to make direct predictions for the JLAB experiments. Darwin, and now CSD3, at Cambridge have proved ideal for this work since we can store intermediate quark propagators for re-use at multiple Q2 values.


Stellar Structure and Nucleosyntheis

Movie showing the evolution of the velocity magnitude of a 10243 resolution run of Carbon-burning shell (video credit: Dr Raphael Hirschi). This movie shows the turbulence in the central convective zone and gravity waves in the top and bottom stable regions. In addition, this movie nicely shows the Kelvin-Helmhotz (shear) instability developing at the top boundary due to the fast horizontal motion of the convective flow when it is deflected and “de”celerated by the boundary. Results from Cristini et al 2016, in prep.Very high resolution movie of the C-shell.


COSMOS Consortium: Flashes of light on dark matter

Flashes of light on dark matter

Figure 1: The “cosmic web”, image credit Matteo Veil.

Simulations of the dense cosmic forest in intergalactic space illuminated by distant quasars have revealed new insights into the properties of dark matter that have recently been published in Physical Review Letters. The study and its results further support the theory of Cold Dark Matter, which is composed of particles that move very slowly, while, for the first time, they highlight the incompatibility with another model, the Fuzzy Dark Matter (FDM), for which dark matter particles have larger velocities.In order to study these dark matter properties, COSMOS Consortium investigators analysed the interaction of the “cosmic web” – a network of filaments made up of gas and dark matter present in the whole Universe – with the light coming from very distant quasars and galaxies. Photons interacting with the hydrogen of the cosmic filaments create many absorption lines called the “Lyman-alpha forest”. The atomic interaction of photons with the hydrogen present in the cosmic filaments is used to study the properties of the cosmos and of the dark matter at enormous distances. On the basis of the results obtained, they were able to infer some of the characteristics of the particles that compose the dark matter. In particular, evidence showed for the first time that the mass of the particles, which allegedly compose the dark matter according to the FDM model, is not consistent with the Lyman-alpha forest observed by the Keck telescope (Hawaii, US) and the Very Large Telescope (European Southern Observatory, Chile). The data, instead, support the scenario envisaged by the cold dark matter model.

GRChombo goes public!

In February 2018, GRChombo, a new COSMOS Consortium code for numerical general relativity simulations targeted at theoretical physics applications, was made open source. It is developed primarily by a consortium of physicists from Argonne National Laboratory, Cambridge University, Göttingen University, King’s College London and Queen Mary University of London, in close collaboration with the COSMOS Intel Parallel Computing Centre at DAMTP Cambridge. GRChombo is maintained by a collaboration of numerical relativists and cosmologists with a wide range of research interests, from early universe cosmology to astrophysics and mathematical general relativity. Further information can be found at, including details of how to obtain and use the code. The GRChombo team hopes that this will be a useful resource for other researchers wishing to investigate the effects of strong gravity in fundamental physics. The next release of GRChombo (version 1.1) will shortly incorporate in-situ or “on the fly” visualization capabilities for AMR grids which were developed as part of the Software-Defined Visualization collaboration involving Intel, TACC and COSMOS IPCC (see



A new model of black hole growth and AGN feedback: an analytic model that explains the phenomenology seen in the EAGLE simulations as well as the observed colour bimodality in the galaxy population. In galaxies less massive than 3×1010 M, young stars and supernovae drive a high entropy outflow. Above a halo mass of 1012 M the outflow ceases to be buoyant and star formation is unable to prevent the build-up of gas in the central regions, triggering a strongly non-linear response from the black hole. Its accretion rate rises rapidly, disrupting the incoming supply of cool gas and starving the galaxy of fuel. The host galaxy makes a transition to the red sequence.

Large scale distribution of the hot, X-ray emitting gas (purple) and stars (yellow) in one of the C-EAGLE cluster regions. The filamentary structure surrounding the cluster halo can clearly be seen. The image is approximately 10 Mpc across. Credit: David Barnes and the Virgo Consortium.

A new theory for the nature of the galaxies that reionised the Universe, based on the EAGLE simulations, in which bursts of star formation can create channels through which ionising photons escape. It was shown that carbon-enhanced metal poor (CEMP) stars are the siblings of the stars that reionised the Universe. Today, these may be found in the intracluster light in groups and clusters.

Construction and analysis of synthetic HI rotation curves of dwarf galaxies in the APOSTLE cosmological simulations of Local Group analogues. Using the same techniques employed by observers (e.g. tllted ring modelling), for each individual galaxy a large diversity of rotation curves is obtained, depending on the orientation of the line of sight. These variations arise from non-circular motions in the gas. As a result, the presence of a core in the dark matter distribution is sometimes incorrectly inferred when there is actually a central cusp, as predicted in the simplest Λ-CDM models.

Completion of the C-EAGLE set of 30 galaxy cluster resimulations using the EAGLE code (including examples with self-interacting dark matter). Their X-ray properties, such as the spectroscopic temperature soft-band luminosity, metal content and Sunyaev-Zel’dovich properties, are in agreement with the observed relations.

The present-day cosmic large scale structure on a scale of 25 Mpc. Circles indicate dark matter haloes resimulated in E-MOSAICS; the mages on the right are zooms showing the distribution of gas, stars and globular clusters in one such galaxy. The latter are coloured to indicate whether they formed in-situ within the galaxy or were brought in by accreted satellites. The images at the bottom show the evolution of the galaxy from the epoch of reionisation at z=10 to the present day.

A new implementation of radiation feedback on galactic scales whereby star formation and the effects of their radiation field on the evolution of the gas are calculated self-consistently.

The first simulations that self-consistently follow the formation and evolution of star clusters and their host galaxy in a full cosmological setting. These “E-MOSAIC” simulations graft the MOSAICS semi-analytic model of star cluster formation and evolution into the EAGLE simulations (see Figure 1).

Simulations of galaxy formation using the EAGLE code in models where the dark matter consists of (warm) sterile neutrinos. Large regions of the sterile neutrino parameter space can be ruled because not enough satellite galaxies are produced.

An explanation for the mass discrepancy-radial acceleration relation (MDAR) of disc galaxies in the context of Λ-CDM and its verification in the EAGLE and APOSTLE hydrodynamics simulations, demonstrating that the MDAR does not require MOND.

Figure 1: Projections around one of our haloes at z~6.4, showing (top) gas column density, (middle) Minimum temperature, and (bottom) hydrogen ionization fraction, showing clear self-shielding.

Using artificial neural networks to constrain the halo baryon fraction during reionization. Radiative feedback from stars and galaxies has been proposed as a potential solution to many of the tensions with simplistic galaxy formation models based on Λ-CDM, such as the faint end of the UV luminosity function. The total energy budget of radiation could exceed that of galactic winds and supernovae combined, which has driven the development of sophisticated algorithms that evolve both the radiation field and the hydrodynamical response of gas simultaneously, in a cosmological context. In Sullivan, Iliev and Dixon (2018) we probed self-feedback on galactic scales using the adaptive mesh refinement, radiative transfer, hydrodynamics, and N-body code Ramses-RT. Unlike previous studies which assume a homogeneous UV background, we self-consistently evolved both the radiation field and gas to constrain the halo baryon fraction during cosmic reionization. We demonstrated that the characteristic halo mass with mean baryon fraction half the cosmic mean, Mc(z), shows very little variation as a function of mass-weighted ionization fraction. Furthermore, we found that the inclusion of metal cooling and the ability to resolve scales small enough for self-shielding to become efficient leads to a significant drop in Mc when compared to recent studies. Finally, we developed an Artificial Neural Network that is capable of predicting the baryon fraction of haloes based on recent tidal interactions, gas temperature, and mass-weighted ionization fraction. Such a model can be applied to any reionization history, and trivially incorporated into semi-analytical models of galaxy formation.

Figure 2: Halo baryon fraction, in units of the cosmic mean,The colour bars show (top to bottom) the strongest tidal force exerted over the last dynamical time, mass-weighted hydrogen ionization fraction, mass-weighted gas temperature (in units of the virial temperature) and the virial ratio (ratio of kinetic to potential energies). The dashed lines show the mean baryon fraction in bins of logarithmic halo mass, where the error bars correspond to the standard deviation.We use the existing halo catalogue data, as it is presented in the left-hand column, for our network predictions. The ANN is trained using only our fiducial simulation, and is then tested against our other models, showing excellent agreement.


Breaking resonances between migrating planets

One of the major breakthroughs in astronomy in the last decade has been the discovery that extra-solar planets are ubiquitous. We now have a variety of ways of detecting exoplanets, and our telescopes have now found thousands of planetary systems. Perhaps the most surprising discovery was that most known extra-solar planetary systems look nothing like the Solar System: we see giant planets orbiting their stars in just few days; planets on highly eccentric orbits; and even systems where multiple rocky or icy planets are packed into very closely-spaced orbits. The latter type of system, discovered by the Kepler mission, are particularly interesting, as these planets orbit too close to their host stars to have formed in their current locations. The conventional picture of planet formation suggests that such planets formed much further from their stars, and then “migrated” inwards due to gravitational interactions with their parent protoplanetary gas disc. Planets of different mass migrate at different rates, however, and in multi-planet systems this process invariably leads to neighbouring planets becoming trapped in so-called mean-motion resonances. These resonances are special locations in planetary systems, where adjacent planets’ orbital periods are integer multiples of one another (i.e., a 2:1 resonance occurs when the outer planet has exactly twice the period of the inner planet). When in such a resonance the planets’ feel one another’s gravity much more strongly than they would do otherwise. Planets which become “trapped” in these locations are therefore expected to remain in resonance indefinitely, but this is at odds with what we see in real systems: most exoplanets in multiple systems are not in resonance with their neighbours. This stark discrepancy between theory and reality has led astronomers to question the very basics of our understanding of how planetary systems form.

Figure 1: Numerical simulation of a pair of “super-Earth” planets migrating though a protoplanetary disc, performed on the DiRAC-2 Data Analytic cluster. The simulation followed the planets’ migration into, and eventually out of, the 2:1 mean-motion resonance, showing that such resonances can be readily broken. (Figure adapted from Hands & Alexander 2018.)

In the last few years astronomers have suggested a wide variety of ways around this problem, but as yet no solution has gained acceptance. The Theoretical Astrophysics Group in Leicester has recently used high-resolution numerical simulations to study how planets in multiple systems migrate through their parent protoplanetary discs. The computational power provided by DiRAC allowed us to follow the evolution of multi-planet systems for many thousands of orbits, capturing the planet-disc and planet-planet interactions in exquisite detail (see Figure 1). These simulations show that under certain conditions (conditions which are prevalent in many planet-forming discs) the gravity between the disc and the planet can overcome the gravity between the planets, breaking the system out of the resonance. This mechanism appears to operate efficiently, with planets escaping resonance after as little as a thousand years. This potentially allows multi-planet systems to migrate into their observed tightly-packed orbits without becoming “trapped” in resonances, and offers an elegant answer to the question of how to assemble the tightly-packed planetary systems observed by Kepler.


Extreme QCD: Quantifying the QCD Phase Diagram

Our collaboration uses the DiRAC supercomputers to simulate the interaction of quarks, the fundamental particles which make up protons, neutrons and other hadrons. The force which holds quarks together is Quantum ChromoDynamics, “QCD”. We are particularly interested in the behavior of QCD as the temperature increases to billions and even trillions of Celsius. These conditions existed in the first moments after the Big Bang, and are recreated (on a much smaller scale!) in heavy ion collision experiments in CERN (near Geneva) and the Brookhaven laboratory (in New York State).

The intriguing thing about QCD at these temperatures is that it undergoes a substantial change in nature. At low temperatures, QCD is an extremely strong, attractive force and so it’s effectively impossible to prise quarks apart, whereas at temperatures above the “confining” temperature, Tc, it is much weaker and the quarks are virtually free and the hadrons they once formed “melt”.

We study this effect by calculating the masses of protons and other hadrons and their “parity partners”, which are like their mirror-image siblings. At low temperatures, hadrons have lower masses than their partners. Above Tc, due to the fact that QCD acquires more symmetry, hadrons and their partners may become equal in mass.

Figure 1: Plot of the “R”-value as function of temperature.

Our results are summarized in the plot where we show the “R” value which measures how close the masses of the partners are: R=1 corresponds to very different masses, and R=0 means that they have identical masses. For a variety of hadrons including the proton (N), delta (Δ) and omega (Ω) as the temperature increases above Tc, we see that R is close to zero implying that the parity partners have near equal masses. This confirms the long-held, but unconfirmed suspicions regarding the nature of QCD. Drilling deeper into this data, it’s clear from the figure, that the Ω hadron partners are the least close in mass. This can be understood because the constituent quarks in Ω are the “strange” quarks which are significantly heavier than the virtually massless “up” and “down” quarks inside the N and Δ. This extra mass means that the symmetry mechanism mentioned does not work as well.


Hadron structure

Hadrons fall into two classes: baryons, such as the familiar proton, and mesons, for example the pion. Deep Inelastic Scattering was the first direct evidence that baryons really were made up of much smaller constituents, “partons”, i.e. quarks and gluons. Investigating the structure of hadrons has become a vital component of many experiments at (for example) JLAB and the LHC.

Figure 1: First results of the method leading to a reconstruction of a parton distribution function (left). Ratio of the electric and magnetic form factors as function of momentum transfer (right).

In [1] we have proposed an alternative method of measuring the parton distribution in hadrons. Parton distribution functions tell us how the momentum of a nucleon is shared out between its partons. The traditional lattice method of measuring parton distributions proceeds through the operator product expansion. In that approach lattice operators are constructed which correspond to moments of the distribution – each moment needs a different operator, with its individual renormalisation constant. Our new approach directly involves simulating scattering on hadrons – in essence we are repeating a deep inelastic scattering experiment on the lattice. First results of the method leading to a reconstruction of a parton distribution function (left panel) are shown above.

Other investigations are designed to determine the distribution of electromagnetic currents in hadrons due to the presence of quarks. These can be parameterised by two functions (electric and magnetic form factors). For most of the second half of the 20th century, it was believed that they both depended on momentum in the same way, so that their ratio was constant. However recent experiments have shown that at higher momentum transfer the ratio decreases and perhaps passes through zero. To investigate this theoretically is difficult (due to the larger momentum transfer). We have developed a method, [2], that allows this range to be considerably extended, approaching the region where the sign might change, as shown in the right panel.

Bibliography: [1] A. J. Chambers et al., Phys. Rev. Lett. 118 (2017) 242001. [2] A. J. Chambers et al., Phys. Rev. D96 (2017) 114509.


Lattice QCD tests Standard Model at high precision

Figure 1: Fermilab g-2 Experiment (Photo: Wikipedia/Glukicov)

We have recently studied the properties of the electron’s heavier sibling, the muon. Its anomalous magnetic moment, aμ, quantifies how vacuum fluctuations contribute when it interacts with an external magnetic field. With accuracy of the order of 1ppm, aμ is one of the most precisely determined quantities in experimental as well as theoretical particle physics. Since its value is sensitive to contributions from yet unknown new physics, the persistent 3−4σ tension between experiment and theory is generating a lot of interest within the particle physics community. Is this really a first sign of new and yet unknown physics? Two new experiments at Fermilab (USA, see Figure 1) and J-PARC (Japan), respectively, are trying to address this question. Both efforts are expected to reduce the uncertainty on the experimental value by a factor of around four, down to 140ppb. This puts much pressure on the theory community to match this precision – only the combination of theory and experiment will be able to answer this important question.

Figure 2: Summary of results for the HVP [4] – RBC/UKQCD 2018 are the results addressed in the text – note the tension with the “No new physics data point”. Green points are other lattice determinations, magenta points are phenome-nological determinations.

On the theory side the error budget is dominated by hadronic uncertainties, in particular uncertainties on predictions for vacuum polarisation effects due to the strong interaction. DiRAC has enabled the first-ever computation of these effects starting from first principles, i.e. from the Standard Model Lagrangian by allowing simulations of Lattice Quantum Chromodynamics including effects of Quantum Electrodynamics. In a series of papers [1-4] researchers from the Universities of Edinburgh and Southampton, together with their collaborators in the US, have developed the required theoretical understanding and the numerical and algorithmic techniques needed to complete their ambitious and state-of-the-art program. In Lattice Quantum Chromodynamics one performs the calculation by constructing a discrete four-dimensional space-time grid (the lattice) on which one numerically solves the QCD equations of motion. Such lattice QCD simulations are the only known method to address the above questions without having to rely on ad hoc assumptions. This type of computation crucially relies on access to the world’s fastest parallel supercomputers like DiRAC’s Extreme Scaling system.

There is still some way to go until a full theoretical computation will match experimental precision. The DiRAC researchers have however developed a method that allows to get very close to the required precision by combining complementary experimental data and their numerical simulations. It will be exciting to compare to new experimental results expected by end of 2018!

Bibliography:[1] Calculation of the hadronic vacuum polarization disconnected contribution to the muon anomalous magnetic moment, Blum et al., Phys.Rev.Lett. 116 (2016) no.23, 232002 [2] Lattice calculation of the leading strange quark-connected contribution to the muon g–2, RBC/UKQCD Collaboration, Blum et al., JHEP 1604 (2016) 063 [3] Isospin breaking corrections to meson masses and the hadronic vacuum polarization: a comparative study, Boyle et al., JHEP 1709 (2017) 153 [4] Calculation of the hadronic vacuum polarization contribution to the muon anomalous magnetic moment, Blum et al., arXiv:1801.07224


Black Hole highlight
Modelling black hole spins and jets in the new era of gravitational wave astronomy

Active galactic nuclei (AGN), with powerful nuclear emission that can outshine the entire host galaxy, represent some of the most extreme physical environments in the Universe. At the heart of these systems reside rapidly growing supermassive black holes (SMBH). They act as efficient gravity engines releasing vast amounts of energy through accretion, over most of cosmic time. It is believed that these energy outbursts may change not only the properties of gas in the vicinity of the AGN, but transform the whole host galaxy, expelling material all the way into intergalactic space.

Figure 1: Column density rendering of a self-gravitating gaseous disc surrounding an equal mass, 106 M sol SMBH binary. Gas streams excited by the time-varying gravitational potential of the binary develop from the cavity edge to the SMBHs. These streams provide fresh material for accretion and affect the binary dynamics through gravitational torques, extracting angular momentum and reducing the SMBH separation. Simulations have been performed on DiRAC supercomputers at the Universities of Cambridge, Leicester and Durham.

There are, none the less, many outstanding questions on the connection between SMBHs and their host galaxies, and numerical simulations have a major role in advancing our knowledge and providing testable theories that future observations may prove or rule out. For example, it is well established that accretion discs play a fundamental role in modulating mass accretion onto SMBHs and driving their spin evolution, with profound implications for relativistic jets and gravitational waves from coalescing SMBH binaries. Their contribution is, however, customarily neglected in large-scale galaxy formation simulations due to the inherent complexity.

Researchers at the Institute of Astronomy (IoA) and at the Kavli Institute for Cosmology (KICC) at the University of Cambridge, have taken a major step forward in solving this issue by developing a novel numerical model within the moving-mesh code AREPO, which self-consistently describes mass feeding through an accretion disc, governing both the evolution of mass and spin of a SMBH. In a recently submitted paper (Fiacconi, Sijacki, Pringle, MNRAS submitted, arXiv:1712.00023), they showed that SMBHs in spiral galaxies, like in our own Galaxy, are expected to be highly spinning, while SMBHs in more massive elliptical galaxies likely have lower spins. IoA/KICC researchers are now exploring the pairing of SMBH spins in binaries (see Fig. 1), which will represent the main cosmological source of gravitational waves for the LISA space observatory (launch 2034).

SMBHs spins are thought to be closely related to the jet phenomena, which are often observed to inflate giant lobes of relativistic plasma in galaxy clusters. Understanding the propagation and interaction of jets with the surroundings by means of hydrodynamical simulations is challenging because of the large dynamic range involved, spanning over 8 orders of magnitude in spatialscales. Building on an earlier implementation of a super-Lagrangian refinement scheme that allows unprecedented spatial resolution close to SMBHs (Curtis & Sijacki, 2015, MNRAS, 454, 3445), IoA/KICC researchers have developed a new jet model for AREPO capable of injecting realistic, highly beamed jets (see Fig. 2; Bourne & Sijacki, 2017, MNRAS, 472, 4707). They found that the inflation of lobes can heat the intra-cluster medium (ICM) through a number of processes, although jets seem to be unable to drive significant turbulence within the ICM in agreement with the gas kinematics observed by the Hitomi satellite in the centre of the Perseus

Figure 2: Simulated AGN-driven jet propagating through the ICM. Each panel shows a different gas property in a slice through the jet. Also shown in the middle panel is the Voronoi cell structure, highlighting the wide dynamic range afforded by the super-Lagrangian refinement scheme: smallest cell size is ~10pc while the size of the computational domain is 20Mpc. Simulations have been performed on DiRAC supercomputers at the Universities of Cambridge, Leicester and Durham.

The next immediate aim is to combine the new models of accretion and jet feedback, as jet power and direction are expected to be directly related to the BH spin. With the upcoming Hitomi replacement, X-ray Astronomy Recovery Mission (XARM, launch 2021), and future Advanced Telescope for High Energy Astrophysics (Athena, launch 2028) on the horizon, IoA/KICC researchers’ simulations on DiRAC supercomputers will provide testable predictions and help interpret these new observations.


Exeter: The formation of massive stars

Figure 1: Density cut through the simulation. The colour scale shows density (scaled between 10−19 and 10−14 g/cc). The arrows show the direction of the velocity field, and their colour the speed (between 0 and 100 km/s).

Stars form as molecular clouds gravitationally collapse and fragment. The formation of the most massive stars (generally understood to be those with masses greater than twenty times the mass of the Sun) is complicated by the presence of the intense protostellar radiation field, which produces outward forces which may hinder accretion and even lead to outflows. Incorporating the effects of the radiation field on the gas is a complex problem, involving a combination of microphysical effects (dust heating, scattering, ionisation) that require a polychromatic treatment of a radiation field.

Figure 2: The line profile (upper left), integrated intensity map (upper right in Jy km/s/pixel), moment 1 map (lower left) and position velocity diagram (lower right in mJy/pixel), for the CH3CN J = 13 → 12, K = 3 line.

At Exeter we have developed a novel radiation hydrodynamic (RHD) code that treats the problem using a Monte Carlo method, splitting the radiation field into many millions of photon packets that propagate through the computational domain, undergoing absorptions, scatterings and re-emissions. We have run this code on DiRAC’s Complexity system to simulate the growth of a massive star. We find that the star forms via stochastic disc accretion and produces fast, radiation-driven bipolar cavities. The evolution of the envelope infall rate and the accretion rate on to the protostar are broadly consistent with observational constraints. After 35 kyr the star has a mass of 25 solar masses and is surrounded by a 7 solar-mass Keplerian disc of 1500 au radius (see Figure 1). Once again these results are consistent with those from recent high-resolution studies of discs around forming massive stars.

We are able to construct so called-synthetic observations from our RHD models. For example, molecular line simulations of a methyl cyanide (CH3CN) transition compare well with observations in terms of surface brightness and line width. These results indicate that it should be possible to reliably extract the protostellar mass from such observations (see Figure 2).


Radiative Transfer Simulations for Type I Supernovae

Figure 1: Comparison between synthetic nebular spectra of W7 calculated with ARTIS (coloured by emitting ion) and the model of Fransson & Jerkstrand (2015; black line).

Between tens and hundreds of days after explosion, the material ejected from Type Ia supernovae become optically thin at most wavelengths. At this time, the observed spectra are dominated by emission from forbidden radiative transitions of singly- and doubly-ionised Fe, Co, and Ni. With almost all of the 56Ni having decayed prior to the nebular phase, most of the energy injection is from the decay of 56Co. The decays of 56Co release gamma rays and positrons, which Compton and Coulomb scatter with thermal electrons, resulting in a population of high-energy (~MeV) non-thermal electrons. These non-thermal electrons are especially important for the ionisation balance, as the rate of non-thermal ionisation is generally higher than photoionisation for the most important ions at this time.

Figure 2: Comparison between normalised ARTIS synthetic spectra for W7 and sub-Chandrasekhar models (normal and high metallicity), and the observed nebular spectrum of SN2011fe from Mazzali et al. (2015).

The ARTIS code (Sim 2007, Kromer & Sim 2009) is a multi-dimensional radiative transfer code that uses the Monte Carlo method with indivisible energy packets (Lucy 2002) to compute spectra for modern supernova explosion simulations. To extend the range of validity to late times, we have (2016/17) developed a new non-LTE population/ionisation solver that includes non-thermal ionisation, excitation, and heating processes, and a comprehensive new atomic database with forbidden transitions and modern photoionisation cross sections. This allows us, for the first time, to calculate spectra out to one year post-explosion (early studies with ARTIS were typically limited to less than one month). Extensive testing has verified that our code can reproduce previous one-dimensional work (see Figure 1).

Taking advantage of DiRAC HPC resources, we have calculated nebular spectra for a range of explosion scenarios, including both Chandrasekhar- and sub-Chandrasekhar-mass explosions (Nomoto et al. 1984; Iwamoto et al. 1999; Sim et al. 2010). Results substantiate the potential of the [Ni II] 7378 feature to distinguish between explosion scenarios (Shingles et al.).


Understanding our Milky Way

The most important new result in 2017 was our (Debattista et al. 2017) demonstration that all the trends observed in the Milky Way’s bulge are produced by the evolution of the bar. Trends such as the X-shape of the metal-rich bulge stars, which is absent in the metal-poor stars, the vertical metallicity gradient, the old age of bulge stars, the weak bar of the oldest stars compared with the strong bar in slightly younger stars, and the varying kinematics of stars with metallicity, have all, in the past, been interpreted as being due to the bulge having formed partly in situ and partly accreted. Our simulations show that all these trends result purely from secular evolution of the bar itself (i.e. the bulge is part of the bar), with no hints of an accreted component other than the stellar halo, at about 1% of the Milky Way’s total stellar mass. Our results have resolved 20 years of conflicting interpretations of the origin of the Milky Way’s bulge and represent an important step in understanding it. We explained the physics of this evolution via a new concept, kinematic fractionation: the morphological separation of stellar populations on the basis of their kinematics at bar formation. The simulations predict strongly X-shaped metallicity maps in external edge-on boxy/peanut-shaped bulges, a prediction which has since been confirmed using MUSE observations (Gonzalez et al. 2017), establishing a novel link between the Milky Way and galaxies in general.

Figure 1: Metallicity map.


Protoplanetary disks

Figure 1: The CO emissivity profile of the extensive disc around the young star IM Lupi, demonstrating that the observed profile is best fit by the photoevaporative flow driven by photoevaporation with a low ambient ultraviolet radiation field.

We have completed a study (Rosotti et al 2017) involving the longest term integration of a planet + protoplanetary disc system that has ever been conducted (totaling 3 x 105 planetary orbits) with a view to understanding whether discs can drive significant planetary eccentricity. This is a fundamental issue for interpreting data on exoplanetary eccentricity. We have shown that there is a complicated exchange of eccentricity between the disc and the planet on these timescales which can be understood in terms of two secular modes that are subject to long term damping and pumping. While there is still much to be understood about the processes responsible for mode damping and pumping, our study has shown that eccentricity growth cannot be determined via the sort of short timescale (< 1000 orbit) simulations that have been conducted previously. These GPU accelerated calculations pave the way for a more comprehensive exploration of long term eccentricity driving of planets in discs.

In Haworth et al 2017 we have presented the first calculations that demonstrate that discs may lose significant mass from their outer regions as a result of photoevaporation by ambient far ultraviolet radiation, even in regions where the ambient radiation levels are very low. This study involved a novel code, TORUS-3DPDR which solves for radiative transfer in dust, hydrodynamics and photochemistry within the outflowing wind and which demonstrates that the extended CO emission observed over many 100s of AU from the disc IM Lupi is well explained by a photoevaporative wind.


Excited Charmonia and Charmonium Resonances from Lattice QCD

Figure 1:  Examples of the results: left is the spectrum in the isospin-1 hidden-charm channel with JP=1+ (J is spin, P is parity), right is the isospin-0 doubly-charmed channel with JP=1+. Boxes are the computed masses and one-sigma statistical uncertainties; horizontal lines show the meson-meson energy levels expected if there were no meson-meson interactions with an adjacent number indicating the degeneracy if it is larger than one (lattice QCD calculations are performed in a finite volume and so the spectrum of meson-meson energy levels is discrete). The number of energy levels in the computed spectrum is equal to the number of expected non-interacting meson-meson levels and they all lie close to the non-interacting levels. This suggests that there are only weak meson-meson interactions and no strong indications of a bound state or narrow resonance in these channels.

Experiments have observed a number of puzzling charmonium-like states in recent years. Of particular interest are states with ‘exotic flavour’, i.e. those with flavour quantum numbers which cannot arise from solely a quark-antiquark pair and so signal physics beyond a conventional meson. One proposed explanation for such states is that they are tetraquarks, compact states containing two quarks and two antiquarks. In Cheung et al [JHEP 1711 (2017) 033, arXiv:1709.01417] we developed methods to investigate the existence of tetraquarks using first-principles lattice QCD computations. The spectrum of states was calculated in a selection of exotic-flavour channels (isospin-1 hidden-charm and doubly-charmed channels), including a number where tetraquarks are expected in phenomenological models. There were no strong indications for any bound states or narrow resonances in the channels studied and so no sign of tetraquark candidates.


Thousands of exoplanets have been discovered in the recent years, most of them are gas giant and hundreds appear to be rocky. Many of these exoplanets have very short orbital period, hence hot atmospheres. The hot temperatures are expected to turn these atmospheres into a high pressure steam, thus having very different spectroscopic signatures than cooler objects, which will influence interpretation of the atmospheric observations. At present interpretation is severely impacted by the lack of corresponding spectroscopic data. The massive number of absorbers in the atmospheres of exoplanets also affects the cooling and therefore the evolution of the young hot objects; comprehensive molecular data is also crucial to model these processes.

Figure 1: A giant planet that orbits close to its sun (Illustration: NASA/ESA/G. Bacon (STScI)).

Our group ExoMol specialises in and leads the field of producing theoretical spectral data (so-called line lists) for molecules of atmospheric, astrophysical, and industrial importance applicable for high temperature applications. The ExoMol database already contains over 300 billion transitions for over than 50 molecular species. In order to accomplish this data intensive task, we use some of the UK’s most advanced supercomputers, provided by the Distributed Research utilising Advanced Computing (DiRAC) project and run by the University of Cambridge. Our calculations require tens millions CPU hours, hundreds thousands GPU hours, and up to 6 TB of RAM, the processing power only accessible to us through the DiRAC project. We are in a unique position to provide line lists for more complex molecules (such as large hydrocarbons), comprising up to 1012 lines.

Methane is a major absorber in hot Jupiter exoplanets, brown dwarfs and cool carbon stars. There is a huge demand for a reliable hot methane line list and there are several groups working towards this objective. As part of our project we have generated a new 34 billion methane line list (known as 34to10), which extends our original 10to10 line list to higher temperatures (up to 2000 K). Even larger line lists are now being completed (for silane, ethylene, acetylene, propene, methyl chloride etc).

Billions of transitions make their direct (line-by-line usage) application in radiative transfer calculations impractical. To mitigate this big data issue we have developed a new, hybrid line list format, based on the idea of temperature-dependent absorption continuum. The new hybrid scheme has being now applied to larger hydrocarbon’s with line lists containing hundreds billions of transitions.


Galactic Archaeology in the Gaia Era

Impacts of a flaring star-forming disc and stellar radial mixing on the vertical metallicity gradient. We ran an N-body simulation of the evolution of a barred disk galaxy similar in size to the Milky Way, using our original Tree N-body code, GCD+. We first assumed the thin disk has a constant scale-height, zd, independent of the radii, ”constant zd model”. We then assigned the metallicity of the thin disk particles to follow the radial metallicity gradient of [Fe/H] = 0.5−0.08× R with a dispersion of 0.05 dex with no vertical metallicity gradient at t = 0. Left two panels of Fig. 1 show the initial (at t = 0) and final (t = 8 Gyr) vertical metallicity distribution for the thin disk component of this model (constant zd model). Interestingly, after 8 Gyr of evolution, the vertical metallicity gradient becomes positive, which is inconsistent with the negative vertical metallicity gradients observed in the mono-age population of the thin disk (e.g. Ciucă, Kawata et al. 2018, MNRAS, 475, 1203).

Figure 1: Vertical metallicity distributions of the thin disk component for constant scale-height, zd, model (left two panels) and flaring model (right two panels) in the radial ranges of 6 < R < 10 kpc at t = 0 and 8 Gyr. The white dots and the error bars at t = 8 show the mean and dispersion of [Fe/H] at different height, |z|, and the solid line shows the best fit straight line. The dotted line shows the vertical metallicity gradient assumed at R = 8 kpc at t = 0.

We presented that a flaring thin disk is one possibility to remedy this issue. We made ”flaring model” with a 2nd thin disk component whose vertical scale height increases with radius. This flaring thin disk component mimics a mono-age population that is born in a flaring star-forming region. We assigned the metallicity following the same radial metallicity gradient as our constant zd model to the flaring disk particles. Right two panels of Fig. 1 show the initial (at t = 0) and final (t = 8 Gyr) vertical metallicity distribution for the flaring thin disk component of this model. In this case, the metal poor stars initially in the outer disk become dominant at high vertical height at every radii, which can drive a negative vertical metallicity gradient. Therefore, if mono-age populations of the Milky Way thin disk stars are formed within a flaring star forming disk, each mono-age thin disk population can have a negative vertical metallicity gradient. Then, older population has a steeper negative vertical metallicity gradient, because they have more time to experience more radial mixing. This can explain the steeper negative vertical metallicity gradient observed in the older population of the thin disk with the combined data of Gaia DR1 and RAVE DR5 (Ciucă, Kawata et al. 2018).


Linking the nuclear interaction to the structure of the heavy elements

Figure 1: Successful comparison between the experimental charge distribution of 36S and our simulations, along with the predicted charge profile of 34Si.

Ab initio calculation of the potential bubble nucleus 34Si. The possibility that an unconventional depletion (referred to as a “bubble”) occurs in the center of the charge density distribution of certain nuclei due to a purely quantum mechanical effect has attracted theoretical and experimental attention in recent years. Based on a mean-field rationale, a correlation between the occurrence of such a semibubble and an anomalously weak splitting between low angular-momentum spin-orbit partners has been further conjectured. Energy density functional and valence-space shell model calculations have been performed to identify and characterize the best candidates, among which 34Si appears as a particularly interesting case. While the experimental determination of the charge density distribution of the unstable 34Si is currently out of reach, we studied the potential bubble structure of 34Si on the basis of the state-of-the-art ab initio self-consistent Green’s function many-body method and using HPC resources at CEA Saclay (France) and the DiRAC Complexity cluster (UK). Demanding that the experimental charge density distribution and the root-mean-square radius of 36S be well reproduced, along with 34Si and 36S binding energies, only leaves one model specific model of the Hamiltonian as a serious candidate to perform this prediction. In this context, a bubble structure, whose fingerprint should be visible in an electron scattering experiment of 34Si, is predicted. The occurrence of a bubble structure in the charge distribution of 34Si is convincingly established on the basis of state-of-the-art ab initio calculations.


Understanding the origin of strong galactic outflows

Understanding the origin of strong galactic outflows is one of the key issues in galaxy formation theory. However, due to our incomplete picture of stellar feedback processes and the lack of numerical resolution, much remains uncertain about how gas is blown away from star-forming regions in galaxies. Using high-resolution radiation-hydrodynamic simulations of a gas-rich isolated disk, we investigated the impact of radiation feedback and supernova explosions on the star formation histories and gas outflow rates. In particular, we have introduced radiation pressure from resonantly scattered Lyman α photons for the first time, and found that it efficiently destroys giant molecular clouds in the early stage of bursty star formation episodes in a metal-poor dwarf galaxy. As a result, supernovae explode in low-density environments, carrying a large amount of gas that is comparable to the star formation rate. We also demonstrated that this early feedback self-controls the gas outflow rate by regulating a bursty star formation history, as opposed to some recent claims that multi-scattering photons can drive extremely massive outflows. As Lyman α photons scatter very efficiently in a dust-poor medium, it may also have played a crucial role in suppressing the formation of extremely metal-poor globular clusters in the early Universe.

Figure 1: radiation-hydrodynamic simulations of a metal-poor dwarf galaxy. The left panels display the projected distribution of gas (top) and stellar density (bottom). The right panels show that the majority of Lyman α radiation is generated (top) and scattering (bottom) inside star-forming regions. Lyman α feedback suppresses the formation of a massive stellar bulge by blowing gas out from the galaxy centre. In order to cover a wide dynamic range from several hundred kpc to a few pc, the simulations were carried out on the high-performance computing facility, DiRAC, employing ~10 million adaptively refined grids.

Categories: Highlights