2015 Highlights

A vacuum is never completely empty space but teems with particles that are created fleetingly by quantum fluctuations in energy and then disappear. The heavier the particle, the less time its brief existence can last. In that time, however, the particle can interact with the much lighter particles in our everyday world and thereby reveal its existence through tiny discrepancies in the properties of these particles from that expected in the Standard Model. The magnetic moment of the muon shows such a discrepancy, a tantalizing 25±9 parts in 10¹⁰. To pin this down more accurately a new experiment, E989 at Fermilab near Chicago, will reduce the experimental uncertainty by a factor of 4. An improved theoretical uncertainty from the Standard Model to match this needs lattice QCD calculations. This year HPQCD has developed a new method using DiRAC that will enable the accuracy needed to be achieved.

Figure 1. Determination of meson decay constants, compared to experimental results, where available. HPQCD, arXiv: 1503.05762, 1408.5768.

The muon is a heavier cousin of the electron with the same electric charge and spin, so that it carries a magnetic moment. The ‘gyromagnetic ratio’, g, measures the ratio of magnetic moment to spin (in units of e/2m). Naively, g=2 but this ignores the interactions of the muon with the cloud of virtual particles discussed above. The anomalous magnetic moment, a_μ = (g – 2)/2, can be measured very accurately from the difference between the circulation and spin precession frequencies as muons fly round a ring transverse to a magnetic field.

The biggest effect in the Standard Model result for a_μ comes from Quantum Electrodynamics where very high order calculations can be done. The largest theoretical uncertainty is from Quantum Chromodynamics and the process of Fig. 1 in which a muon radiates a photon that briefly becomes a quark-antiquark pair. This ‘hadronic vacuum polarisation’ (HVP) contribution needs to be calculated to better than 0.5% in lattice QCD to reduce the theoretical uncertainty. Our new method (arXiv:1403.1778, Phys. Rev D89:114501) makes this possible, having achieved 1% for the first time for the strange quark HVP contribution. Our method has low systematic errors and DiRAC enabled high statistics on realistic gluon field configurations. The up/down quark calculation is now underway. Errors will be reduced by a further factor of 3 as we increase statistics by a factor of 10 in 2015.

Type Ia Supernovae (SNe Ia) are among the most energetic phenomena in the Universe and are important for many different astrophysical research fields (e.g. theories of binary stellar evolution, chemical enrichment of galaxies, and cosmology). These spectacular events are believed to be associated with thermonuclear explosions of white dwarf stars but answers to the questions of when, why and how the explosions are triggered remain unclear.

Figure 1. Flux and polarization spectra (black) predicted 7 days before maximum for the merger model of Pakmor et al. (2012) viewed from a specific orientation. For comparison, grey lines show observed spectra of SN 2004dt at the same epoch.

In the past few decades spectropolarimetry has 
played a leading role in the search for a comprehensive picture of SNe Ia (Wang & Wheeler 2008). 
The particular value of this technique is its ability to probe the geometry of the explosion: linear polarization, attributed to electron scattering, has been detected in SNe Ia, suggesting that asymmetries 
in the progenitor system and/or explosion mechanism are 
present. However, interpreting polarimetric data is challenging and careful modelling is required to make the connections between observed spectropolarimetry and specific ejecta geometries. Our current research aims to address this challenge by making polarization prediction from modern multi-dimensional explosion models.

We recently developed a new Monte Carlo technique for use in our radiative transfer code (Kromer & Sim 2009) that allows us to calculate polarization signatures for hydrodynamical explosion models (Bulla et al. 2015). Compared to simpler techniques, this method leads to a substantial reduction in the Monte Carlo noise and thus in the required computing time (a factor of ~ 50). The combination of this new technique and access to the super-computing facilities provided by DiRAC has now made it feasible to compute polarization spectra for fully three-dimensional explosion models using simulations that run for only a few days.

In a recent paper (Bulla et al. 2016), we presented polarization calculations for the violent merger of two carbon-oxygen white dwarfs. The violent merger model has received considerable attention in the past few years, but our study is the first to investigate its spectropolarimetric signatures in detail. We demonstrated that the particular merger model of Pakmor et al. (2012) provides a remarkably good match with some highly-polarized SNe Ia (see Fig. 1), but is too asymmetric to reproduce the low levels of polarization commonly seen in many other SNe Ia. We are now extending our studies to alternative explosion scenarios: this will allow us to determine how spectropolarimetry can be most effectively used to distinguish and test explosion theories.

Bibliography: Bulla et al. 2015, MNRAS, 450, 967; Bulla et al. 2016, MNRAS, 455, 1060; Kromer & Sim 2009, MNRAS, 398, 1809; Pakmor et al. 2012, ApJL, 747, L10; Wang & Wheeler 2008, ARA&A, 46, 433.

Galaxies come in a variety of shapes. Nevertheless, in a simplified view of galaxy formation, one can classify them either as spirals or ellipticals. However, the orientation and ellipticity of these spirals and ellipticals is not random; it depends on their formation history and their environment. In particular, tidal forces from the large-scale structure of the Universe (see figure below) play a role in determining correlations between ellipticity, orientation and the underlying (dark) matter density field. Such correlations are referred to as ‘intrinsic alignments’.

On the other hand, future imaging surveys from the Euclid satellite and the Large Synoptic Survey Telescope, aim to measure per cent level perturbations on galaxy ellipticities caused by a different effect: the deviation of photons from their original path due to the gravity of these same large-scale structures. This ‘gravitational lensing’ effect constitutes a powerful probe of the nature of the mysterious dark energy that drives the accelerated expansion of our Universe.

Figure 1. Left panel: 28 Mpc slice through the Horizon-AGN simulation cube (140 Mpc on a side), at redshift z=0.5. Right panel: 10x zoom on a region located in the bottom left part of the cube. The different colours indicate gas density (in green), gas metallicity (in blue) and stars (in red). Galaxies are identified as clusters of stars. Their shapes are measured to investigate intrinsic ellipticity correlations that contaminate the weak gravitational lensing signal used by astronomers to constrain the nature of dark energy.

Therefore, the quest for improving our understanding of intrinsic alignments to avoid them contaminating dark energy observables is paramount. Thanks to the DiRAC facility, our group (Horizon-UK) lead an analysis of unprecedented statistical scope and detail of galaxy shape and spin alignments, mining the huge dataset from the state-of-the-art hydrodynamic cosmological simulation Horizon-AGN. These pioneering results are reported in Chisari et al, 2015, MNRAS, 454, 2736.

We have used the DiRAC-2 Complexity cluster to obtain a mass estimate of the Carina dwarf spheroidal (dSph) satellite of the Milky Way at the time when it first fell into the gravitational field of our Galaxy (about 6 Gyr ago). This is the first primordial mass estimate that has been obtained directly from observed data without recourse to matching with cosmological simulations.

Figure 1. Comparison of our mass estimate for Carina (black and red data points with error bars) with predictions for the relationship between stellar mass (y-axis) and dynamical mass (x-axis) obtained via abundance matching of observed galaxies and simulated halos. (Figure reproduced from Ural et al., 2015.)

The dSphs which surround the Milky Way are vital laboratories in which to test theories of galaxy formation. Not only are they are the lowest mass stellar systems known to contain dark matter but they also populate the faint end of the galaxy luminosity function. Over the past 15 years, several issues have been identified which appear to place cosmological galaxy formation simulations in tension with observations: (1) simulations which contain only dark matter over-predict the numbers of objects around a Milky Way-like galaxy by at least an order of magnitude (the “missing satellites” problem; Moore et al., 1998); (2) simulated haloes exhibit central density cusps in their dark matter profiles, while observations are often claimed to be more consistent with uniform density central cores (the “cusp-core problem); (3) simulations predict the existence of dark haloes whose masses should have allowed them to form stars but which are apparently missing from the census of observed galaxies (the “too big to fail problem”; Boylan-Kolchin et al., 2012). The resolution of all three issues requires knowledge of the masses of the dwarf satellites both at the present epoch and at the time of their formation.

Despite the availability of large observed data sets containing measured velocities for 1000s of stars in several of the brighter Milky Way satellites, previous analyses of their mass distributions have made many simplifying assumptions. In particular, most models assume that the dwarf galaxy is in dynamical equilibrium and isolated from external perturbations. For small Milky Way satellite galaxies, both these assumptions are of questionable validity. N-body simulations provide a natural way to relax these assumptions and explore the impact of the Milky Way on the evolution of a dwarf galaxy. We have combined a Markov-Chain-Monte-Carlo algorithm with an N-body simulation code to explore the parameter space of models of the Carina dwarf galaxy using 20,000 medium-resolution simulations.

We identified the Carina dSph as an interesting target for this initial study due to the large amount of available kinematic data and the fact that there were already hints that its outer regions are being disturbed by the Milky Way (e.g. Muñoz et al., 2006). Our algorithm was able to constrain the present-day mass of Carina to be (7±3)×10⁷ solar masses. More interestingly, we were also able to constrain Carina’s mass at the time it first fell into the Milky Way to be 3.9 (−2.4;+3.9)×10⁸ solar masses. Our mass estimate is significantly lower than the halo mass generally associated with dSphs of comparable luminosity to Carina. In Ural et al. (2015), we suggested that this may be evidence that galaxy formation becomes stochastic in low-mass haloes.

The calculations for this project consumed approximately 1M core hours on the Complexity cluster of DiRAC-2, as a very significant number of simulations were required to verify the efficacy of the method using artificial data sets before it could be applied to the observed data with confidence. By using the Complexity cluster we were able to run multiple MCMC chains in parallel which greatly reduced the overall time required to complete the project. The work is reported in Ural, Wilkinson, Read & Walker, 2015 (Nature Communications, 6, 7599).

Quasars are extremely bright objects that reside in the centres of galaxies. They are believed to consist of a supermassive black hole (some as large as a billion times the mass of our own Sun) that is surrounded by gas. These black holes grow in size by swallowing their reservoir of gas and, in doing so, release an enormous amount of energy into their host galaxies. This energy affects the growth of each galaxy, changing where and when stars form and playing an important role in the different appearances of galaxies that we observe in the Universe today.

Figure 1.

Such quasars are very rare in the present day Universe, but were much more active in the past. Because of their brightness, some of the most distant objects we can see are some of the most massive quasars, existing barely 1 billion years after the Big Bang. The presence of these objects raises many questions: how were these objects able to grow to such a large mass so quickly? What do the galaxies that host such objects look like? How does the growth of such a massive black hole affect the surrounding environment?

Using the DiRAC supercomputers at Universities of Cambridge, Leicester and Durham, researchers at the Institute of Astronomy and Kavli Institute for Cosmology, Cambridge simulated the growth of a galaxy containing a quasar similar to those we observe. These results are summarised in a recent Letter to MNRAS by Curtis and Sijacki, accepted on December 10, 2015. Vital to the work was the use of a new computational technique to significantly increase the resolution of simulations around the central black hole, which Curtis and Sijacki developed over the last year (2015, MNRAS, 454, 3445) thanks to unique computational resources provided by DiRAC. This novel technique enabled the capture of the flow of gas around black holes with unprecedented accuracy both in term of mass and spatial resolution.

The work by Curtis and Sijacki provides new insights into what we might hope to observe in the near future around these distant quasars with the James Webb space telescope and the Atacama Large Millimeter/submillimeter Array. Long thin structures, called filaments, transport gas from the outskirts of the galaxy all the way to the central region. These filaments grow a disc of gas that orbits around the black hole. The gas in this disc is relatively cold, which allows it to become dense enough to form stars at a dramatic rate. A fraction of the gas escapes the disc and falls onto the black hole, releasing large quantities of energy. This drives a powerful flow of super-heated gas out of the central galaxy, oriented in such a way that it does not destroy the surrounding disc, which is the key novel result of the model.

In the coming years there will be many observations of quasars in the early Universe. The ways in which they compare with the predictions of Curtis and Sijacki will provide clues to the physical processes that we do not fully understand, allowing deeper insight into the role black holes play in shaping our Universe.

A hot topic in cosmology is to test gravity on scales well beyond our Solar system. 2015 marks the centenary of Einstein’s publication of his General theory of Relativity (GR), which is still our standard gravity theory. Although GR has been tested rigorously in the Solar system, its application in cosmology is still largely an extrapolation. Previous cosmological observations did not have the precision to test it accurately, but the situation is changing quickly in light of the future galaxy and cluster surveys such as Euclid, LSST, DESI and eROSITA, all of which aim to test GR as a primary science objective. Over the next decade, these observations will determine key cosmological parameters with percent-level precision, and bring the cosmological tests of GR to a similar precision for the first time. They will also address another critical question in modern cosmology – whether the mysterious accelerated expansion of our Universe is due to some unknown exotic matter (dark energy) or a different behaviour of gravity on cosmological scales.

Figure 1. The velocity field in a slice of the simulation box of the standard (left) and non-standard (right) gravity model. Colours are used to illustrate the amplitudes of the velocity field. Filamentary structures, towards which matter fall, are clearly visible, as well as the model differences around these filaments. It is by differences such as shown here that different gravity models can be tested against observational data.

To fully exploit those billion-dollar surveys and to achieve their scientific potential of testing gravity, it is crucial to be able to accurately predict the observable signatures of non-standard gravity theories. This is a numerically intensive problem, which depends critically on powerful supercomputing resources such as DiRAC-2. Indeed, DiRAC has enabled us to make significant progresses over the past few years [1-11]. These efforts have led to developments of state-of-the-art algorithms and simulation codes such as the first ever parallelised N-body code for such theories, ECOSMOG [1-3], as well as a recent code [8] implementing a novel algorithm to trace light rays across cosmological simulations on the fly. These codes have been used widely in studies of this field, and ruled out popular non-standard models such as Galileon gravity [7]. More recently, tests on DiRAC led an improvement of the performance of ECOSMOG by over 10 times for some theories [9], making it possible to study such theories by large and high-resolution simulations in the future.

The numerical studies can also have implications to other fields. For example, the Galileon model features nonlinear differential equations similar to the famous Monge-Ampere equation in the optimal transportation problem, and its study has led to a new stable algorithm to solve this type of equation. These previous studies have fully equipped us [8-10] with efficient tools to make accurate predictions for a number of cosmological observables, including weak gravitational lensing, redshift space distortions and clusters of galaxies, for which we expect precise observational data in the next decade. They have prepared us for accurate cosmological tests of gravity and will help to shed light on the understanding of the origin of the accelerated expansion of our Universe. More exciting progresses are expected for the next few years.

Bibliography: [1] Li B., Zhao G., Teyssier R., Koyama K., JCAP, 01, 051 (2012), [2] Li B., Zhao G., Koyama K., JCAP, 05, 023 (2013), [3] Li B., et al., JCAP, 11, 012 (2013), [4] Li B., et al., MNRAS, 428, 743 (2013) [5] Jennings E., et al., MNRAS, 425, 2128 (2012); Hellwing W. A., et al., Phys. Rev. Lett., 112, 221102 (2014) [6] Shi D., et al., MNRAS, 452, 3179 (2015); Cai et al., MNRAS, 451, 1036 (2015) [7] Barreira A., et al., JCAP, 08, 028 (2015); JCAP, 08, 059 (2014) [8] Barreira A.,et al., JCAP, submitted, [9] Barreira A., Bose S., Li B., JCAP, submitted, [10] He J., Li B., Phys. Rev. Lett., submitted

Accretion is an important process relevant to various fields in astrophysics, from star formation to the study of compact binaries and Type I supernovae. Continuous theoretical and numerical efforts have been devoted to understanding the interaction between the inner edge of an accretion disk and the stellar surface. Complementary to these works, there are also many studies devoted to the effect of accretion on the structure and evolution of objects, using either simple semi-analytic approaches or simulations based on one dimensional (1D) stellar evolution codes. But no multi-dimensional numerical study has ever been devoted to the effect of accretion on the structure of accreting objects. Because of this lack of sophisticated modelling, great uncertainties still remain regarding the fraction of accretion luminosity imparted to an accreting object or radiated away, and how the accreted material is redistributed over the stellar surface and the depth it reaches. Early accretion history may impact the properties of stellar objects (luminosity, radius, and effective temperature) in star formation regions and young clusters even after several Myr when the accretion process has ceased. These effects could provide an explanation for several observed features in star-forming regions and young clusters, having important consequences on our general understanding of star formation.

Figure 1. Distribution of tracer particles accreted on a young sun for cold (e1) and hot accretion (e9) cases respectively. Color scale is the radial velocity in units of cm/s and gray scale points show the position of the tracer particles. The cold accretion case shows rapid mixing of the accreted particles within the convective envelope. The hot accretion case shows the formation of a hot surface layer where particles accumulate and do not sink.

The DiRAC Resource, Complexity, has been used to perform the first hydrodynamics simulations describing the multi-dimensional structure of accreting, convective low mass stars. The hydrodynamical simulations were computed using our fully compressible time implicit code MUSIC (Multidimensional Stellar Implicit Code) recently improved with the implementation of a Jacobian-free Newton-Krylov time integration method, which significantly improves the performance of the code. These developments and the benchmarking of the code were made using Complexity. These new accretion calculations have enabled us to understand the effect of so-called “hot” accretion, when the accreted material is characterised by a larger entropy than that of the bulk of the stellar material.

The multi-D simulations show an accumulation of hot material at the surface that does not sink and instead produces a hot surface layer, which tends to suppress convection in the stellar envelope. This enables us to derive an improved treatment of accretion in 1D stellar evolution models, based on an accretion boundary condition. Such treatment of accretion in a 1D stellar evolution code is more physical than the standard assumption of redistribution of the accreted energy within the stellar interior. These new results will be published in Geroux, Baraffe, Viallet et al. (2016, A&A, submitted). The impact of this more physical accretion boundary condition on the evolution of young low mass stars and brown dwarfs and on observations in young clusters will be explored in a follow-up studies. These pioneering simulations performed on the DiRAC facility shows the potential of using time implicit multi-D simulations to improve our understanding of stellar physics and evolution.