Some of the Research Highlights from our users in 2020: –

Massive Discs around low-mass stars

Planets form within flattened discs of material around young stars. How these discs evolve and how they might form planets remains a mystery. In recent years exoplanet detections have shown planets around low-mass stars (stars lower in mass than about half the sun’s) are extremely common, even more common than around sun-like stars. How these discs’ retain sufficient mass to form these planetary systems is an open questions. 

The problem of the requirement for massive discs around low-mass stars is personified by the famous Trappist-1 planetary system that boasts 7 rocky planets that all orbit their 0.1 solar mass star in periods less than a month. Standard models of the disc that formed the Trappist-1 planetary system suggest they should be so massive it was unstable and would have destroyed itself on a short-timescale. However, young stars are much brighter when they are first born, heating this surrounding disc, increasing its temperature and providing pressure support against the gravity that wants to destroy it. 

TJ Haworth led a collaboration as part of the DiRAC project dp100, to use Smooth Particle Hydrodynamics simulations to show that the luminosity emitted by bright stars was strong enough to stabilise these discs, allowing them to retain the mass necessary to form the observed planetary systems.  

Disc structures around low mass stars determined using SPH simulations. The left panel shows a smooth stable disc heated by a young low-mass star. The middle panel shows the disc structure if the low-mass star’s main sequence properties were used. This disc shows spiral structures indicative of gravitational instabilities that rapidly drain it of material. The right panel shows a disc around a low-mass star that is not emitting, this disc has destroyed itself within a few orbital periods. 


Haworth et al. 2020 MNRAS 494 4130

HPQCD: LUV is in the air!

Figure 1. When the Bc emits a W boson to become a J/psi meson, different helicities of W boson affect the rate of decay in different kinematic regions of angles and momenta.

Nuclear beta decay is the most well-known version of a process in which quarks emit weak interaction W bosons and change type or ‘flavour’ whilst bound by the strong interactions inside a hadron. In nuclear beta decay a neutron becomes a proton as one of its down quarks changes into an up quark. Other quarks undergo beta decay too and in the process provide insights into possible new physics beyond the Standard Model (SM). The W boson here is seen through its decay to a lepton and antineutrino. The SM W boson does not distinguish between leptons, except for effects that arise from their masses. This is known as Lepton (flavour) Universality (LU). Hints of lepton universality violation (LUV) have been seen in some beta decays – it is important to study others to pin this down. Accurate lattice QCD calculations of strong interaction effects are critical here for comparison to experiment.

Figure 2. W Helicity amplitudes for Bc to J/psi decay vs squared momentum transfer. Ht appears with the lepton mass and so only contributes when the W decays to the heavy tau.

The beta decay of a Bc meson to a J/psi meson is a promising such process, currently being studied by the LHCb experiment at CERN. A Bc meson has no net spin but the J/psi charm and anti-charm quarks have their spins aligned for a total spin of 1. This means that the different spin alignments (helicities) of the W boson contribute to the angular distribution of the decay products in ways that may allow tests for new physics (Figure 1). HPQCD has calculated the amplitude for this process via different W helicities and covering the full range of possible J/psi momenta for the first time, using lattice QCD (Figure 2; arXiv:2007.06956, 2007.06957).  This allowed us to determine the ratio of rates of production in the SM of tau leptons in the final state compared to mu (from the W) as 0.2582(38) with a 1.5% uncert-ainty. Initial LHCb results hint at disagreement, which might indicate LUV from new physics, but their large uncertainty masks any significance in the comparison. Using our results in the LHCb analysis should enable reduced experimental systematic errors.  In another lattice QCD ‘first’ HPQCD also calculated the amplitudes for Bc to Bs and Bc to Bd beta decay (arXiv:2003.00914). Here a charm quark transforms to an s quark or d quark; these are on LHCb’s ‘to-do’ list to compare to c to s/d processes in D mesons. 

The HPQCD calculations were done on the DiRAC Data Analytic system in Cambridge; this is ideal for the numerically efficient methods we have developed to study b and c quarks using precision lattice QCD on very fine lattices.  

Spontaneous Symmetry Breaking in the 3d Thirring Model

The Thirring Model describes relativistic fermions moving in a two-dimensional plane and interacting via a contact term between covariantly conserved currents. The physical system it most resembles is that of low-energy electronic excitations in graphene. For free electrons at half-filling on a honeycomb lattice, conduction and valance bands form cones just touching at their vertices at two “Dirac points” lying within the first Brillouin zone. Since the density of states vanishes, and the effective fine structure constant is boosted by a factor vF/c≈1/300, where the pitch of the cones vF is the Fermi velocity, the resulting physics is described by a strongly-interacting relativistic quantum field theory, with equal likelihood of exciting electrons or holes. 

Besides possible applications in layered condensed matter systems, the Thirring model is an interesting theory in its own right, and possibly the simplest theory of fermions requiring a computational solution. A major question is whether or not for sufficiently strong self-interaction a bilinear condensate <ψψ> forms, restructuring the ground state and causing a gap to form between the conduction and valence bands resulting in a phase transition from semimetal to insulator – this process is precisely analogous to spontaneous chiral symmetry breaking in QCD. It is anticipated that this can occur if the number of fermion species N lies below some critical Nc: determination of Nc  requires non-perturbative quantum field theory.

In lattice field theory, the problem seems to require a precise rendering of the correct U(2N) fermion global symmetries. We have been using the Domain Wall Fermion formulation, in which U(2N) GInsparg-Wilson symmetries are recovered in the limit that the domain wall separation Ls∞. Simulations have been performed using the RHMC algorithm on 12and 163 systems with Ls ranging from 8 to 48. It turns out that at strong couplings g2 and light fermion masses m recovery of GW symmetry is slow – nonetheless at weak-to-moderate couplings our results for <ψψ> as a function of g2, m obtained at Ls=48 are compatible with an extrapolation to the large Ls limit based on the entire dataset. They are well-fitted by an equation of state which assumes a continuous symmetry-breaking transition with non-mean field theory critical exponents at a critical 1/g2≈0.28 (see figure). The fit is less convincing at the strongest couplings (lying within the broken phase), which we attribute to the  Ls∞ limit not yet being reached.

A major outcome of the project to date is the confirmation that the critical number of flavors for symmetry breaking in the Thirring model Nc>1. This is a significant step towards understanding the formulation of fermion quantum fields in strongly-interacting theories.

Amplification of superkicks in black-hole binaries through orbital eccentricity

One of the most spectacular discoveries obtained from numerical simulations of black holes in the framework of Einstein’s general relativity are the so-called superkicks. These arise from the net emission of linear momentum in gravitational waves during the inspiral and coalescence of two black holes. Just like firing a bullet imparts a recoil on the shooter, the single black hole resulting in from the coalescence will recoil from the preferential emission of gravitational waves in a specific direction. Superkicks are the most dramatic version of this effect and are generated when the black holes are rapidly rotating but in opposite directions. The recoil can then reach magnitudes of up to 3700 km/s. This is particularly remarkable, since this velocity is large enough to kick the black hole out of even the most massive galaxies. Superkicks can thus explain astrophysical observations of black holes displaced relative to the centre of their host galaxy or quasars (luminous objects sourced by black holes) exhibiting significant blue or redhift relative to their environment. Superkicks also affect the assembly of black holes throughout cosmological history.

In our work, we have shown that moderate values of the orbital eccentricity can even further enhance the superkick magnitudes by about 25%, leading to kicks well above 4000 km/s. Our results have consequences for the ongoing programs of observing sources of gravitational waves. Larger kicks are able to eject a larger fraction of black holes from their hosts, stellar clusters or galaxies, and can reduce the expected stochastic gravitational-wave background targeted in pulsar timing observations. The superkick effect is illustrated by the figure below, taken from one of our simulations at the time shortly after the two black holes have merged into one.

Figure: Snapshot of the gravitational-wave signal generated by an eccentric black-hole binary in the so-called superkick configuration. This snapshot shows the cross section of the wave signal around the time of the black-hole merger inside the plane perpendicular to the orbital plane. Note that the waves are stronger in the upper half than the lower half. This asymmetry corresponds to a net emission of linear momentum and a recoil of the black hole resulting from the merger. This is illustrated in the three panels on the right which display a zoom-in of the central region at different times: the time Tm around merger (as in the large left panel) as well as 60 and 120 time units later (measured in units of the black-hole mass). The (nearly circular) black line in these panels represents the black-hole horizon. Note that the black hole moves downwards as time progresses; this motion is the kick effect.

Hydrodynamic simulations demonstrate how spiral arms in protoplanetary discs probe vertical temperature structure.

HD100453 is a notable protoplanetary disc where a binary companion drives a pair of prominent spiral arms, detectable using a number of dust and gas tracers. In Rosotti et al 2020 1 we demonstrated how the different pitch angles of the spirals traced in the submillimetre continuum and in scattered light can be used to probe the vertical temperature structure of the disc. The idea is simply that the pitch angle (angle between the spiral arm and the local tangential direction) measured using any particular diagnostic depends on the temperature of the region that produces this emission and is expected to be larger in warmer regions of the disc. Scattered light emission derives from small dust grains several pressure scale heights above the disc mid-plane while sub-millimetre continuum emission derives from large grains settled towards the disc mid- plane. Thus it is possible to use forward modeling to simulate discs with a variety of vertical temperature structures in order to reproduce the observed spiral morphologies.

The lefthand panel shows the spirals in HD100453 at 0.88 mm (blue contours) and in scattered light at 0.6 μm, demonstrating that the latter spirals are more open and thus likely to be associated with warmer conditions (the pitch angles are respectively 6° and 19°). We performed 3D hydrodynamical simulations with the FARGO3D code of a disc subject to the dynamical effects of the companion star (shown as asterisk in figure) and then post-processed these simulations to generate simulated observations using the Monte Carlo radiative transfer simulation, RADMC3D. The polar emission plots (right panel) clearly showed that it is impossible to reproduce the spiral morphologies in both tracers using a disc that is vertically isothermal and that instead the surface layers (from which the scattered light emission originates) has to be around four times hotter than the mid-plane regions.

This result is plausible given that the upper layers are exposed to direct irradiation from the central star. However, our simulations provide the first hydrodynamical demonstration for vertical temperature stratification in protoplanetary discs.

Exploring fundamental fields with strong gravity

Two of the most pressing questions in cosmology are the nature of dark matter, and how the Universe began. When we say we do not know what dark matter is, what we mean is that we cannot characterise it as a particle – the mass, spin, self interaction and couplings to standard model matter are unknown. When we say we do not know how the Universe began, we mean that accepted paradigms like inflation still leave many questions unanswered at the earliest moments after the Big Bang. This project has pushed forward our knowledge in these areas, highlighting the model dependence of inflationary models and clarifying the behaviour of dark matter environments around black holes.

A key highlight this year has been the extension of our studies of light dark matter candidates to the binary case, where we have found density enhancements between the black holes that impact on the merger part of the gravitational wave signal. Such effects are small and need to be distinguished from other parameters of the binary, but have the potential to shed light on the particle nature of dark matter.

Dark matter density around binary back hole mergers for bosonic candidates.

Discovery+ Universe Unravelled and advanced in-situ visualization

The series, Universe Unravelled with the Stephen Hawking Centre, premiered in November 2020 on Discovery+, coinciding with the UK launch of this new digital platform.  In 25 short episodes the series explores what we know about the Universe, including how massive objects warp the fabric of spacetime and create gravitational waves. The series employs stunning graphics from Cosmos Consortium simulations performed on STFC DiRAC systems, which were produced in collaboration with Intel’s Advanced Visualization and Rendering team.

An extensive reel of the visualization work created with Intel can be viewed on the CTC YouTube channel:

This animation compilation includes visualizations of gravitational waves from black hole mergers and radiation from cosmic string networks, as shown below:

The Cosmos Consortium has a longstanding collaboration with Intel on in-situ visualisation, a new paradigm that analyses data “on the fly” as it is generated, rather than post hoc as was traditionally the case. It is an essential step in the approach to exascale computing and these new capabilities, developed collaboratively also with Kitware, are now available in the scientific data visualization application Paraview, which operates across multiple platforms from laptops to supercomputers. The underlying Intel OSPRay ray-tracing libraries were awarded a Technical Achievement Academy Award on 13th February 2021 because of their widespread use in movie animation, including the Avengers movies.    

The Universe Unravelled series is available on the new Discovery+ service here:

Convection influenced by rotation: simulations and theory

All stars, and many planets, carry energy by convection in some parts of their interiors. The heat transport associated with this convection, which must be calculated in order to construct a consistent stellar or planetary structure model, is influenced by rotation or magnetism.  But a quantitative understanding of the effects these have on the convection, and so ultimately on stellar/planetary evolution, is still lacking.  A central problem is that the convective dynamics occur on timescales that are much shorter than the evolutionary ones.  Hence, numerical simulations that can adequately resolve the convective flows cannot be evolved over evolutionary timescales in order to investigate long-term effects on stellar structure.

Example snapshot of temperature and vertical flows from a rotating 3D simulation at a latitude of 45°.  The convective plumes tend to align with the rotation axis.  [Currie et al. 2020]

In recent work, we have circumvented this problem by using a hierarchy of approaches, ranging from 3D simulations to 1D models. With the aid of DiRAC resources, we conducted a large survey of 3D simulations in localised, Cartesian domains tilted at some angle with respect to the rotation vector (described in Currie, Barker, Lithwick, Browning 2020, MNRAS, 493, 5233).

These box simulations served as idealised representations of a small part of a rotating star or planet, situated at various latitudes. We used these to analyse the rich variety of phenomena that occur as the rotation rate and latitude are changed, and to compute how the temperature gradient established by the convection varies with these parameters. Broadly, rotation makes the heat transport less efficient, leading to steeper temperature gradients. We compared the results of these calculations to expectations from semi-analytical theory: in particular, we showed that a new multi-mode “rotating mixing length” theory (based on earlier work by Stevenson 1979) provides a reasonably good description of the dynamics at most latitudes and rotation rates. This semi-analytical model could be incorporated into 1D stellar evolution codes.  Separately, we have carried out 1D stellar evolutionary models (with the open-source code MESA), informed by our 3D simulations, to begin studying what effects the stabilising influence of rotation (or magnetism) has on stellar structure (e.g., Ireland & Browning 2018). These results have, for example, provided constraints on the mechanisms responsible for the apparent “inflation” of some low-mass stars, which seem to be systematically larger than standard (non-rotating, non-magnetic) models predict.    

A new mechanism for “dark matter heating” in ultra-faint dwarf galaxies

For many years, we have known that the amount of dark matter at the centres of nearby dwarf galaxies appears to be less than models predicted. This has generated significant interest because it might point to dark matter being more complex than the “cold” thermal relic particle assumed so far.

Figure 1: Ultra-faint dwarf galaxies simulated as part of our DiRAC EDGE project. This plot shows the dark matter density of a patch of a simulated Universe from which these dwarfs are selected (centre), and their individual dark matter densities at the end of the simulation (inset panels).

Recent observations appear to point to a particular solution to the above puzzle: “dark matter heating”. Repeated inflow and outflow of gas during star formation causes the inner gravitational potential to fluctuate with time, causing the dark matter particle orbits to migrate outwards. In this model, we expect an anticorrelation between the inner dark matter density in dwarf galaxies and the amount of star formation they have experienced. This indeed appears to be the case.

However, a puzzle remains. There is mounting evidence for a low dark matter density also at the centres of at least some of the very smallest dwarf galaxies in the Universe: so-called “ultra-faint” dwarfs. These have formed so few stars that “dark matter heating” models do not predict any appreciable lowering of their inner density. In our EDGE project, we set out to run some of the highest resolution simulations of ultra-faints dwarf galaxies to date to address this problem (see Figure). Run on DiRAC, these simulations reach a mass resolution of ~100 solar masses per simulation particle, and a spatial resolution of ~10 light years. A key novelty is the use of “genetically engineered” initial conditions. These allow us to forensically change and explore the impact of late/early assembly on these dwarfs in their full cosmological context. Using this, we found a rather surprising result: late forming ultra-faints were continuously lowering their inner dark matter densities even long after star formation ceased. We tracked this down to a whole new mechanism that can drive dark matter heating in dwarfs: minor mergers. Late assembling dwarfs experience many mergers throughout their lifetimes. As these pass close to the centre on plunging orbits, they can also fluctuate the potential and excite dark matter heating. This could, then, explain the puzzlingly low dark matter density at the centres of some ultra-faint dwarfs.

Project P.I. : Prof. Justin I. Read, University of Surrey 

Authors (lead author in bold): Orkney, Matthew D. A. (Surrey); Read, Justin I. (Surrey); Rey, Martin P. (Lund); Nasim, Imran (Surrey); Pontzen, Andrew (UCL); Agertz, Oscar (Lund); Kim, Stacy Y. (Surrey); Delorme, Maxime (CEA, France); Dehnen, Walter (ZAH, Heidelberg)


Modelling the Epoch of Reionisation in the Simba Simulations 

The Epoch of Reionisation is the period in the early Universe between about 200 million and 1 billion years after the Big Bang. During this epoch, the very Birst stars and galaxies formed and started to eat away at the fog of neutral gas that pervaded the cosmos. Detecting and characterizing these Birst galaxies in one of the main goals of upcoming facilities such as NASA/ESA’s James Webb Space Telescope, and the European Large Telescope, and is already being probed with the Hubble Space Telescope and ESO’s Atacama Large Millimetre Array (ALMA). Simulations such as Simba, run on cosma-7 at Durham, are a critical tool for situating such data within a cosmological context and exploring the underlying physical drivers.


Harvard graduate student Xiaohan Wu, with project PI Romeel Davé at Edinbugh, used the Simba simulation suite to explore the physical and observable properties of EoR galaxies as will be seen by Hubble and Webb. Wu et al. (2020) showed that Simba produced galaxies with sizes and photometric colours in very good agreement with existing Hubble observations, and then made predictions for the sort of improvements that will be seen using Webb (to be launched in 2022). The image below shows redshift 6 simulated galaxies “observed” with Hubble in 3 bands (top row), and with Webb (bottom row). The improvement is stark, and enables detailed studies of the assembly history of early galaxies.


In a separate work, Cornell student Daisy Leung worked a group including PI Davé to use the same Simba suite to study the far-infrared properties of EoR galaxies as would be observed with ALMA (Leung et al. 2020). To do so, she applied a new far-IR line radiative transfer code called SÍGAME to Simba outputs, creating the Birst ensemble of [CII] line emission images from a cosmological simulation in the EoR, and the Birst predictions of the [CII] luminosity function that might be detected with a large ALMA program, or future proposed far-IR observatories such as NASA’s Origins. These results illustrate Simba’s capabilities to make making multi-wavelength predictions for EoR galaxies that can help guide and interpret present and future observations of the Birst billion years.

Small-scale magnetic reconnection from an anisotropic turbulent cascade

The solar wind is an expanding plasma that is formed in the solar corona and spreads across the solar system. As the plasma travels away from the Sun, its electromagnetic and kinetic power spectra exhibit fluctuations over a wide range of temporal and spatial scales. The energy stored in the larger fluctuations cascades to smaller fluctuations. At the end of this cascade, particle heating occurs. It is not clear which mechanisms are responsible for the ultimate heating of the particles. The reason for this lack of knowledge is the general possibility for multiple phenomena to exchange energy between particles and electromagnetic fields depending on the plasma conditions and the nature of the fluctuations. In this context, the outstanding question is whether and how the plasma fluctuations reach the smallest scales to fulfil the conditions for dissipation.

Magnetic reconnection is a fundamental plasma process that modifies the topology of the magnetic field at different scales and increases the kinetic and thermal energy of the particles. The level of turbulence in a plasma before and during magnetic reconnection affects the effectiveness of this process as an energy exchange mechanism. In addition, magnetic reconnection (like the turbulent cascade) transports energy to small scales. Therefore, we ask the science question: “How does magnetic reconnection relate to small-scale plasma phenomena and the turbulent cascade?”

Volume rendering of the electric current density J in one of our DiRAC simulations of anisotropic Alfvénic turbulence. The turbulence forms elongated magnetic structures, which are called “magnetic flux ropes”. The flux ropes are unstable, and magnetic reconnection results from their non-linear interaction. Particles are accelerated during the reconnection process. Artificial spacecraft trajectories through the simulation domain can be compared with in-situ observations in the solar wind. 

To address this problem, the use of peta-scale simulations has proven to be indispensable. The computational power provided by the high-performance computing facilities at DiRAC enables us to perform rigorous studies that resolve the small-scale dynamics, also known as kinetic scales, while preserving the history of the large-scale cascade. In particular, particle-in-cell (PIC) simulations solve the evolution of the plasma based on first principles. This type of simulation accounts for phenomena that only reveal themselves in kinetic theory. With DiRAC’s capabilities, we have come a big step forward to understanding turbulence and magnetic reconnection in plasmas.

We initialise our turbulence through the collision of counter-propagating plasma waves in our simulation box. These waves interact with each other and create turbulence, which is consistent with the observed turbulence in the solar wind. This interaction forms, after some time, magnetic and electric current structures. Some of the magnetic structures are so-called “flux ropes” which undergo magnetic reconnection and accelerate particles.  We define a set of criteria to identify where in the simulated plasma reconnection occurs. With these criteria, we find several reconnection sites in our simulation domain. We select one extensive reconnection site and study its properties in great detail. We find that this event is complex and asymmetric unlike the idealised Harris current-sheet system often used to study magnetic reconnection.  One great advantage of our simulations is that we can “fly” an artificial spacecraft through our simulation box. This spacecraft records data in the same way as a real spacecraft would record data in the solar wind. The measurements from our artificial spacecraft serve as predictions for measurements of reconnection in the solar wind by the latest solar-system spacecraft Solar Orbiter and Parker Solar Probe which have the appropriate instruments onboard to test our predictions.

Categories: Highlights