Highlights

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

• Massive Discs around low-mass stars
• HPQCD: LUV is in the air!
• Spontaneous Symmetry Breaking in the 3d Thirring Model
• Amplification of superkicks in black-hole binaries through orbital eccentricity
• Hydrodynamic simulations demonstrate how spiral arms in protoplanetary discs probe vertical temperature structure
• Exploring fundamental fields with strong gravity
• Discovery+ Universe Unravelled and advanced in-situ visualization
• Convection influenced by rotation: simulations and theory
• A new mechanism for “dark matter heating” in ultra-faint dwarf galaxies
• Modelling the Epoch of Reionisation in the Simba Simulations
• Small-scale magnetic reconnection from an anisotropic turbulent cascade
• Compton amplitude from the lattice
• Silicate grain growth due to ion trapping in supernova remnants
• Kinematic signatures of gravitational instability
• Orbital evolution of binaries with circumbinary discs
• Rare particle decays on DiRAC super computers
• Guiding new 1D stellar evolution models with detailed 3D hydrodynamic simulations
• Composite Dynamics signatures at the LHC
• The Milky Way total mass profile as inferred from GAIA DR2
• Monte Carlo radiative transfer for white dwarf deflagrations
• Cosmic accretion through massive haloes
• DK I = 0, DK ̅ I = 0, 1 scattering and the D_s0^* (2317) from lattice QCD
• Blandford-Znajek jets in galaxy formation simulations
• Hunting for the sources that reionized the Universe

2020 list of publications

---

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.  

References

Haworth et al. 2020 MNRAS 494 4130

---

HPQCD: LUV is in the air!

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.

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 v_F/c≈1/300, where the pitch of the cones v_F 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 N_c: determination of N_c  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 L_s∞. Simulations have been performed using the RHMC algorithm on 12³ and 16³ systems with L_s ranging from 8 to 48. It turns out that at strong couplings g² and light fermion masses m recovery of GW symmetry is slow – nonetheless at weak-to-moderate couplings our results for <ψψ> as a function of g², m obtained at L_s=48 are compatible with an extrapolation to the large L_s 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/g²≈0.28 (see figure). The fit is less convincing at the strongest couplings (lying within the broken phase), which we attribute to the  L_s∞ 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 N_c>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.

---

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.

---

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 13^th February 2021 because of their widespread use in movie animation, including the Avengers movies.    

---

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.

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.

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)

Paper: https://ui.adsabs.harvard.edu/abs/2021arXiv210102688O/abstract

---

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?”

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.

---

Compton amplitude from the lattice

Understanding the internal structure of hadrons from first principles remains one of the foremost tasks in particle and nuclear physics. It is an active field of research with important phenomenological implications in high-energy, nuclear and astroparticle physics. The static properties of hadrons, from the hybrid structure of quark and meson degrees of freedom at low energies down to the partonic structure at short distances, are encoded in structure functions. The forward Compton amplitude describes the process of virtual photon scattering from a hadron (here a nucleon) and via the optical theorem provides an essential ingredient in the determination of the structure functions.

Lattice QCD, however, has really only been able to probe some limited kinematic corners of the all-encompassing structure functions – primarily being limited to elastic form factors and low(est) moments of leading-twist parton distributions. In a breakthrough paper published in Physical Review Letters [1], we demonstrated how it is possible to extract the forward Compton amplitude from lattice QCD. This is achieved by taking advantage of the efficiency of the Feynman-Hellmann approach and avoids the need to compute 4-point functions. This advance opens a new avenue to study a range of features of hadronic states, including non-trivial correlations in their partonic structure that can provide novel insights into parton dynamics.

In a recent paper [2], we extended this investigation and demonstrated for the first time that the Q²-dependence of the moments of structure functions can be probed directly within a lattice simulation (seen in the left figure below). This breakthrough calculation clearly shows that power-corrections (curvature at small- Q²) can be significant and will have important consequences for attempts to extract the leading-twist parton distributions from experimental data at low- Q². One possible source of power-corrections is due higher-twist effects which go beyond the simple parton picture. Such contributions are suppressed by extra powers of Q², however their size in the small-Q² domain has until now been poorly known. The right figure shows a first attempt [3] at isolating one such higher-twist term for three values of  Q² which clearly shows an enhancement at small Q².

[1] A. Chambers et al., Phys. Rev. Lett. 118 (2017) 242001.

[2] K. U. Can et al., Phys. Rev. D 102 (2020) 114505.

[3] A. Hannaford-Gunn et al., PoS LATTICE2019 (2020), 278; R. Horsley et al., PoS LATTICE2019    (2020), 137.

---

Silicate grain growth due to ion trapping in supernova remnants

Dust material in the clumpy ejecta of supernova remnants is exposed to high gas temperatures and shock velocities. On the one hand, the energetic conditions can cause a significant destruction of the dust grains due to sputtering or grain-grain collisions. On the other hand, the energetic gas ions can penetrate deep into the dust grains. For grain temperatures below ~500 K, the diffusion rate of oxygen  and other heavy ions in silicates is very low and they are trapped once they have intruded into the grain. This process, named as ion trapping, has not been considered so far as a measure to counteract grain destruction by sputtering.

To study this phenomenon we considered a clump impacted by a shock-wave in the oxygen-rich supernova remnant Cassiopeia A. DiRAC HPC Facilities were used to conduct hydrodynamics simulations which were further post-processed with the code Paperboats (Kirchschlager et al. 2019). This allowed us to follow the dust mass and grain size evolution in the shocked ejecta material.

In our new study (Kirchschlager et al. 2020) we showed that penetration and trapping within silicate grains of oxygen, silicon, and magnesium, the same impinging ions  that are responsible for grain surface sputtering, can significantly reduce the net loss of grain material (Fig. 1). We found for a pre-shock gas density contrast between clump and ambient medium of χ = 100 that ion trapping increases the surviving masses of silicate dust by factors of up to two to four, compared to cases where the effect is neglected (Fig. 2). The formation of grains larger than those that had originally condensed is facilitated and allows the presence of micrometre-sized grains in the post-shock medium. For higher density contrasts (χ ≥ 180), we found that the effect of gas accretion on the surface of dust grains surpasses ion trapping, and the survival rate increases to ∼55% of the initial dust mass for χ = 256.

References:
Kirchschlager, Barlow, M. J., &

Schmidt, F. D. 2020, ApJ, 893, p.70-78

---

Kinematic signatures of gravitational instability

Young stars are born surrounded by cold discs or dust and gas. These discs act as a conduit for accretion of mass on to the central star, and also provide the raw material for planet formation. At early times these discs are thought to be sufficiently massive that they can be unstable to their own gravity. This gravitational instability initially leads to the formation of spiral density waves in the disc, and these spirals transport angular momentum very efficiently, increasing the disc accretion rate dramatically. In some cases the instability can also lead to fragmentation of the disc, leading directly to the formation of giant planets or brown dwarfs. Understanding gravitational instabilities is therefore crucial to understanding both star and planet formation. However, in the early stages of star formation discs are typically still enshrouded in their parent molecular clouds, making them difficult to observe, and measurements of the disc mass are plagued by systematic uncertainties (because the molecular hydrogen that comprises the bulk of the disc mass is largely invisible to our telescopes). Moreover, there are numerous other processes (such as disc-planet interactions) which can give rise to spiral structures in discs, and distinguishing these from the spirals caused by gravitational instabilities is challenging. It is therefore still not clear when, where, or even if, real discs become gravitationally unstable.

In recent years new observational facilities, especially the Atacama Large Millimetre/sub-millimetre Array (ALMA), have provided us with astonishing new high-resolution observations of discs around young stars. These discs have uncovered an unexpected wealth of structure in discs – gaps, rings, spirals, shadows – and have given us a huge array of new insights into the processes of star and planet formation. However, in the case of gravitational instabilities we have not been able to exploit these observations fully, because it was not clear what our telescopes should be looking for. By exploiting the power of DiRAC’s Data Intensive at Leicester (DIaL) cluster, we were able to perform a suite of 3-D Smoothed Particle Hydrodynamics simulations of gravitationally unstable discs, and then couple those simulations to radiative transfer calculations to generate synthetic observations of unstable discs. We found that the critical, unambiguous signature of gravitational instability lies not just in the existence of spirals, but rather in their kinematics, and that this signature is readily detectable in ALMA observations. Emission lines from carbon monoxide molecules provide the best probe of disc kinematics. In a “normal” (non-self-gravitating) disc the disc’s rotation results in a well-known “butterfly” pattern in the velocity structure, which is readily seen in ALMA’s so-called channel maps (essentially a series images of the CO emission at a particular velocities). Our DiRAC simulations showed that this is radically altered in self-gravitating discs, with the spiral density waves resulting in a characteristic “wiggle” in the channel maps (see Fig.1). This wiggle is a unique signature of the spirals generated by self-gravity, and represents the smoking gun of gravitational instability in discs around young stars. The race is now on to detect this signal in real observations, and these results will provide crucial new insight into the formation of stars and planets.

These results were published as Predicting the Kinematic Evidence of Gravitational Instability, C.Hall et al., The Astrophysical Journal, vol. 904, article 148 (2020).

---

Orbital evolution of binaries with circumbinary discs

Binaries are common in astrophysics. Stars typically form in binary or multiple systems, and when two galaxies merge their supermassive black holes form a binary system in the merged galaxy’s centre. In each case it is possible for gas to accumulate in circumbinary orbits, and the interaction between this disc of gas and the binary is a recurring theme in astrophysics. Such discs have been recently observed directly in protostellar systems. And they are thought to facilitate the mergers of black holes which creates gravitational wave emission.

Following analytical investigations and pioneering numerical simulation work in the 20th century, the standard picture of circumbinary discs is that they drive the binary to shrink over time. This occurs because there are discrete locations in the disc where the binary orbit resonates with the disc orbits, and these resonances transfer energy and angular momentum from the rapidly rotating object (the binary) to the slower rotating object (the disc). As a result the disc expands and the binary contracts. Circumbinary discs are typically expected to be truncated at a radius of a few times the binary semi-major axis with, in some cases, a small amount of accretion on to the binary components via time-dependent streams from the disc inner edge.

In the last few years, this standard picture has been challenged by numerical simulation results from several groups that report that the binary expands with time upon interacting with a circumbinary disc. These simulations were typically performed with large disc viscosity to reduce computational cost and often in 2D. In this work (Heath & Nixon, 2021) we explore 3D hydrodynamical simulations with more realistic values of the disc viscosity to test whether the binary expands or contracts when interacting with a circumbinary disc. We find that when the disc viscosity is sufficiently large as to overwhelm the binary-disc resonances, the binary expands, and we show that this is due to the accretion of material that comes from orbits with a larger specific angular momentum than the binary orbit. However, we also show that when the simulations are performed with realistic disc parameters, the binary shrinks as the binary-disc resonances dominate the evolution over the accretion of matter in this case. We show that the structure and evolution of the system is dependent on several disc and binary parameters (see Fig. 1), and conclude that in most systems we expect the binary orbit to shrink over time.

---

Rare particle decays on DiRAC super computers

Quarks are the fundamental particles that make up most of ordinary matter, such as protons and neutrons in atomic nuclei. They are bound together by the strong nuclear force, mediated by the exchange of gluons, as described by the theory of strong interactions, Quantum Chromodynamics (QCD). They come in different flavours and are found to change into one another. The study of such flavour-changing transitions could well be the arena in which the first hints of physics beyond the Standard Model are discovered. Increased experimental precision and novel theory calculations are required to better understand them. In this context we note our last year’s Dirac Science Highlight on hadronic K→ππ decays which has in the meantime appeared in Physical Review Letter, where it was highlighted as Editor’s Suggestion.

Within the RBC/UKQCD collaboration with researchers from the UK and USA, we are studying flavour-changing processes that are found extremely rarely in nature and dominated by so-called long-distance effects. Over the past years our collaboration has developed novel theoretical and numerical tools and is now in a strong and unique position to compute the underlying physics.The particular class of processes we are focussing on at the moment are the rare kaon decays K → πl+l− and K → πνν ̄. These Flavour Changing Neutral Current (FCNC) decays are in parallel being studied experimentally at NA62 (CERN) and KOTO (J-PARC). Our calculations rely crucially on the cancellation of power divergences through the GIM mechanism for which inclusion of charm valence quarks on fine lattices is crucial. Because of its longstanding experience in kaon physics and the use of chiral fermions, RBC/UKQCD has a significant advantage to carry out these calculations in a timely fashion that coincides with the NA62 timeline. By repeating our exploratory studies on the most advanced data set generated with the help of DiRAC resources we will make predictions for this decay, which complement NA62’s experimental program.

The problem to be solved is to evaluate quantum field theory correlation functions defined by path integrals, from which particle masses, decay rates and other physical quantities can subsequently be determined. By discretising a finite volume of space time, the path integral becomes an extremely high- but finite-dimensional ordinary integral. A transformation to Euclidean space allows us to treat the action term in the integral as a probability weight and the integral can then be evaluated by importance sampling. Lattice QCD simulations naturally fall into two steps: The first step is to use a Markov chain Monte Carlo to generate an ensemble of field configurations. In the second step correlation functions are evaluated using the previously generated ensemble as input.

Working with a fine enough lattice spacing on a large enough volume and with light quarks at or near their physical mass poses a computational challenge which needs a large-scale HPC facility at the petaflops scale if simulations are to be completed on a competitive time scale. The DiRAC Extreme Scaling Service is ideally suited for this task.

Due to the strong suppression in the Standard Model these processes are highly sensitive to new physics and could very well allow to find first signs of it.

DiRAC resources are now enabling a detailed study which will allow direct comparison to experiment. RBC/UKQCD’s calculation to date is the only one of its kind and therefore highly anticipated.

Bibliography:

• “Lattice QCD study of the rare kaon decay K+→π+ννˉ at a near-physical pion mass”, Christ et al., Phys.Rev.D 100 (2019) 11, 114506, arXiv:1910.10644
• “K+→π+ννˉ decay amplitude from lattice QCD”, Bai et al., Phys.Rev.D 98 (2018) 7, 074509, arXiv:1806.11520
• “Exploratory Lattice QCD Study of the Rare Kaon Decay K+→π+ννˉ”, Bai et al., Phys.Rev.Lett. 118 (2017) 25, 252001, arXiv:1701.02858
• “First exploratory calculation of the long-distance contributions to the rare kaon decays K→πℓ+ℓ−“, Christ et al. , Phys. Rev. D94 (2016) 11, 114516 arXiv:1608.07585
• “Prospects for a lattice computation of rare kaon decay amplitudes II K→πννˉK→πννˉ decays:, Christ et al., Phys.Rev.D 93 (2016) 11, 114517, arXiv:1605.04442
• “Prospects for a lattice computation of rare kaon decay amplitudes: K→πℓ+ℓ−decays”, Christ et al., Phys.Rev.D 92 (2015) 9, 094512, arXiv:1507.03094
• “Standard Model Prediction for Direct CP Violation in K→ππ Decay, Bai et al. Phys. Rev. Lett. 115, 212001

---

Guiding new 1D stellar evolution models with detailed 3D hydrodynamic simulations

We have used the results of 3D simulations of stellar carbon burning (Cristini et al. 2017, 2019) simulated on DiRAC systems to guide the development of a new prescription for convective boundary mixing in stellar evolution models. We applied the phenomenon of turbulent entrainment at convective boundaries, seen in the 3D simulations of convective shells, to the hydrogen-burning convective core of evolutionary models. Shown in the figure (Fig. 6 of Scott et al. MNRAS, 2021) are the evolutionary tracks, in terms of luminosity and temperature, of a selection of models from our grid. These include the previous standard of convective boundary mixing in our evolutionary code (blue curves) and new models with entrainment (red and purple dashed), with dots, plusses and crosses marking 90%, 95% and 99% of the hydrogen-burning lifetime respectively. We compared the coolest edge (rightmost extent) of these curves to both previous model grids using different amounts of mixing (Ekstrom et al. 2012, Brott et al. 2011, dotted lines) and observational constraints (Castro et al. 2014, dot-dashed line). The crosses on the entrainment model tracks, occurring at or near the cool edge of the tracks, follow a trend with mass which more closely follows the shape of the Castro et al. line than the previous model grids. This suggests that the physics included in the entrainment prescription could be used to reproduce observed features which are not seen in standard models. In future work, evolutionary models with improved convective boundary mixing physics will be used to initialise new 3D simulations of convection. For more details on this work, follow the weblink to Scott et al. 2021 article given below

Weblinks: 
MNRAS: https://academic.oup.com/mnras/advance-article-abstract/doi/10.1093/mnras/stab752/6169727
ArXiV: https://arxiv.org/abs/2103.06196

---

Composite Dynamics signatures at the LHC

Understanding which novel interactions can give rise to the emergence of the Standard Model at a low energy scale of about 1 TeV without conflicting with the high precision tests of the latter theory is one of the biggest challenges of current Particle Physics. This problem requires robust input from Quantum Field Theory and phenomenology that can guide current and future experimental searches.  Composite Higgs models are based on a new fundamental interaction which features an emergent Higgs boson which properties are similar to the Standard Model Higgs boson, a key requirement motivated by the current experimental results provided by the LHC experiments.  Among the restricted class of gauge theories that feature all the key requirements, particularly relevant are those based on the symplectic group Sp(2N).

Our project has provided insights in the dynamics of Sp(2N) gauge theories relevant for experimental searches of composite dynamics.  We finalised the first calculation of the low energy parameter of the Goldstone Boson in Sp(2)=SU(2) gauge theory that controls the production of new vector resonances and which is constrained by  vector boson scattering processes at the LHC [1].  Our predictions is currently being used to constrain further the model and to study the constraints provided by future colliders. In addition, our project has performed a benchmark calculation of the scattering properties of the Goldstone boson in the singlet channel[2]. Ultimately, such a quantity will allow to constrain from first principle the Higgs properties in a Composite Higgs scenario.  Furthermore a detailed study of the spectrum of Sp(2N) purely gluonic bound states in the absence of matter was performed[3,4]. 

The study shed light on the dependence of the spectrum on the underlying gauge symmetry, addressing a fundamental question to inform model building. A smooth variation of the spectrum with the number of colours N has been found, showing that the effect of the finite value of N can be seen as a correction to well-defined infinite-N masses of the various states.  

[1] V. Drach et al.,  Scattering of Goldstone Bosons and resonance production in a Composite Higgs model on the lattice,  accepted for publication in JHEP [arxiv:2012.09761]
[2] V. Drach et al., Scattering of Goldstone Bosons in  SU(2) model with N_f=2 fundamental  fermions, Asia Pacific Symposium for Lattice Field Theory, August 2020 (paper in preparation).
[3] E. Bennett et al., Color dependence of tensor and scalar glueball masses in Yang-Mills theories, Phys Rev D.102.011501 [arXiv:2004.11063]
[4] E. Bennett et al.,  Glueballs and Strings in Sp(2N) Yang-Mills theories, accepted for publication in Phys. Rev. D  [arxiv:2010.15781]

---

The Milky Way total mass profile as inferred from GAIA DR2

Over the past few years, there has been a breakthrough in the ability to simulate realistic galaxies and galaxy populations through cosmological hydrodynamical simulations. Some of the main programmes in this area have been conducted using DiRAC facilities.

This kind of simulations not only inform us about the physical processes at work during the formation of galaxies but also allow us also to connect cosmological theory directly to astronomical observations. In this work we use simulations from three of the major cosmological hydrodynamics simulations of recent years, the EAGLE, APOSTLE and AURIGA projects with data from the GAIA satellite which has transformed our view of the Milky Way. We use the simulations to construct a detailed physical model of the Milky Way.

The model galaxy contains 6 baryonic subcomponents: a bulge, a thin and a thick stellar disc, an HI and a molecular gas disc and a circumgalactic medium (CGM), all embedded in a dark matter halo. The simulations have shown that the presence of baryons induces a contraction of the dark matter distribution in the inner regions of the galaxy which we include in our Milky Way model.

We fit our model to the Gaia DR2 Galactic rotation curve and other data and find the best-fit model for the Milky Way. This has a dark matter halo mass, M_DM = 0.97 ± 0.20 × 10¹² M_o, and concentration before baryon contraction of 9.4 ± 1.3, which lie close to the median halo mass–concentration relation predicted in ΛCDM. We provide values for all the baryonic components of the Milky Way.

---

Monte Carlo radiative transfer for white dwarf deflagrations

Type Ia supernovae (SNe Ia) play a number of key roles in astrophysics. These include synthesising heavy elements such as iron, injecting kinetic energy in galaxy evolution and acting as cosmological distance indicators. Despite this, however, the way in which SNe Ia occur is still poorly understood and there is currently no definitive explosion scenario that explains SNe Ia. Additionally, thanks to data from modern transient surveys, it has become clear that while “normal” SNe Ia make up a relatively homogenous sequence of objects, there exist a number of “peculiar” sub-classes within the SNe Ia population. Among the largest of these sub-groups is the set of Type Iax supernovae (SNe Iax), which were estimated by Foley et al. (2013) to make up 31_-11⁺¹⁷% of the SNe Ia population. SNe Iax are a diverse class of objects that are sub-luminous and span a much wider range in brightness compared to normal SNe Ia.

Understanding which explosion mechanism(s) produce different SNe Ia is a key research question. The Chandrasekhar-mass, pure-deflagration model has long been suggested as a possible explosion mechanism to explain SNe Ia. In this scenario, a carbon-oxygen white dwarf accretes material from a companion star until it approaches the Chandrasekhar mass limit. Nuclear reactions are then ignited in the white dwarf’s core, which leads to a thermonuclear runaway. The resulting explosive energy release will gravitationally unbind either part or all of the white dwarf. When it is assumed that the propagation of the nuclear burning front is a sub-sonic deflagration, this model struggles to account for the high luminosities of “normal” SNe Ia. However, studies have shown that pure-deflagration models with variations to their initial conditions can explain a wide variety of the observed properties of SNe Iax.

Our DiRAC project involves using computer simulations of radiative transfer to predict the spectrum of light that would be emitted by theoretical models of supernovae. We use the comparison between these simulations and observed supernovae to better understand the underlying explosion models and gain insight into which explosion models are favoured and disfavoured by the data for each specific sub-class of supernovae observed. The work we highlight here specifically focusses on the Chandrasekhar-mass, pure-deflagration scenario as a possible explosion scenario to explain the SNe Iax class.  

Previous deflagration models of Fink et al. (2014) varied the number of ignition sparks used to initiate the deflagration as a way of varying the strength of the explosion and found that a variety of SNe Iax properties could be explained in this way. However, the previous models did not account for the full range of observed luminosities and, moreover, the multi-spark ignition approach has been shown by Nonaka et al. (2012) and Zingale et al. (2011) to be less probable than single-spark ignition. Therefore, in collaboration with our colleagues at the Heidelberg Institute for Theoretical Studies, we have carried out a new set of deflagration simulations to systematically explore the degree of diversity that can be accounted for with single-spark ignition models by varying the ignition position, central density, and composition. Sample results are shown in the accompanying figure. Our results indicate the deflagration scenario can cover a significant region of the parameter space occupied by observed SNe Iax. In particular, the models may be almost a magnitude fainter than the models in the Fink et al. (2014) sequence, and thus give luminosities extending down to the faintest SNe Iax. This work (to be published in a paper led by our collaborators) has therefore confirmed the full observed diversity in brightness of the SNe Iax class can be explained by pure deflagrations of Chandrasekhar-mass, carbon-oxygen white dwarfs. However, systematic differences between models and data remain and require further investigation. We are therefore currently working on refining simulations with more advanced microphysics (Shingles et al. 2020) to better understand these differences. We intend to publish this follow up work later this year.

This work also used the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

References:

• Fink et al., 2014, MNRAS, 438, 1762
• Foley et al., 2013, ApJ, 767, 57
• Nonaka et al., 2012, ApJ, 745, 73
• Shingles et al., 2020, MNRAS, 492, 2029
• Zingale et al., 2009, ApJ, 704, 196

---

Cosmic accretion through massive haloes

During the hierarchical gravitational collapse of primordial dark matter overdensities, the baryonic content of haloes grows by accreting gas and stars from the intergalactic medium (IGM). This gas is what then feeds the formation of a galaxy at the halo’s centre. The mechanism by which this feeding occurs is not fully understood, as it involves a complex interplay between cooling, shock heating, and feedback from a variety of physical processes in the galaxy itself. All of these processes act in concert in the circumgalactic medium (CGM), which thus plays host to the baryon cycle that regulates the process of galaxy formation.

In the `classical’ analytic model of accretion onto massive dark matter haloes gas shock heats as it infalls, forming a hot atmosphere. Recently however, a significant amount of observational evidence points to a more complex, multi-phase CGM with a substantial cold component across a wide variety of galaxy masses.

In massive haloes, cold gas is thought to bypass accretion shocks and penetrate deep into haloes via filaments that deliver cold gas through the CGM. Numerically this is difficult to model, as the CGM is primarily made up of low-density, multi-phase gas that cosmological simulations do not focus resolution on. Accretion shocks themselves are also numerically broadened in simulations which may affect gas shock heating and energy dissipation at the shocks. These effects can potentially lead to inaccurately modelled properties of CGM gas.

Within this project we introduce a novel computational technique that automatically and on-the-fly boosts numerical resolution around shocks during the simulation. This allows us to better resolve the accretion shocks at the boundary of the CGM and IGM in massive haloes. This is the region that determines if primordial gas filaments can penetrate into the hot halo, and how turbulent motions are generated in the wake of curved shocks.

To do this, we resimulate a massive, high-redshift galaxy cluster progenitor with the successful physical model employed in the FABLE suite of simulations using the moving-mesh code AREPO. Figure 1 shows how our ‘shock refinement’ scheme boosts resolution significantly in and around cosmic filaments and at the accretion shocks, compared with AREPO‘s ‘base refinement’ scheme which is commonly adopted in the literature.

Figure 2 demonstrates clearly the physical effect of our novel method in maps of  Mach number, temperature, metallicity and turbulent velocity. The top panels show that with the ‘shock refinement’ scheme active we see sharper and better defined accretion shocks around the halo. The second row illustrates how the amount of cold gas penetrating into the hot halo is increased with ‘shock refinement’, leading to a much more multi-phase CGM more in line with observations. These cold gas filaments are shown to also have a low metallicity in the third row, indicating how our ‘shock refinement’ scheme allows primordial gas to more easily penetrate deep within the virial radius of the halo. Finally we see that higher numerical resolution and better defined curved accretion shocks lead to a significant boost in turbulent velocities in the halo, which is important when comparing models to the next-generation of X-ray observations from eROSITA and Athena.

Having a well-resolved accretion shock also allows us to create a 3D shock surface, as shown in Figure 3. Here we colour the surface with the Mach number, indicating that a significant amount of the surface area has no associated accretion shock at all. To explain this we also plot the positions of inflowing gas filaments in grey, whose points of intersection with the halo correspond very well with the unshocked regions of the shock surface. This visually demonstrates the filamentary feeding of galaxy growth, and how even in massive dark matter haloes with a hot gaseous atmosphere, cold, primordial gas can be delivered directly to the central galaxy.

Improved physical models and a better understanding of the CGM and the feeding of galaxies are important for a wide variety of current and upcoming observations, including those by MUSE and the Keck Cosmic Web Imager and enhanced Sunyaev-Zel’dovich observations from SPT and ACTpol.

The simulations for this project were carried out on COSMA7 in Durham and Peta4 in Cambridge as part of the DiRAC project dp012. 

---

DK I = 0, D  I = 0, 1 scattering and the from lattice QCD

Published in JHEP 02 (2021) 100 [arXiv:2008.06432]

Much of observable strongly interacting matter can be understood in terms of the quark model whereby hadrons are composed of a quark and an antiquark or three quarks. In recent years significant deviations have been observed from quark model expectations, including unexpected experimental results for hadrons composed predominantly of a charm quark and a strange antiquark. In particular, the  was found much lower in mass than was expected.

Quantum chromodynamics (QCD) is the theory of strongly interacting matter, and in this work, we applied state-of-the-art numerical lattice QCD methods to help understand the enigmatic , and other closely-related systems. We found a bound state in DK scattering in S­­-wave (i.e. with no relative angular momentum) that resembles the experimental state and is strongly influenced by the nearby DK kinematic threshold. One useful feature of lattice QCD is the ability to vary the quark masses, and this can help us better understand why a state has a certain mass and what its composition is. We used two different light-quark masses – the D and K mesons both contain a light quark and so are strongly dependent on this – and found that the state becomes less bound as the light-quark masses are reduced from a heavy value to a value closer to that found in experiment.

The figure shows the computed DK scattering amplitudes for the lighter light-quarks (left) and the heavier light-quarks (right). The more rapid rise in the S-wave amplitude at DK threshold in the left plot is due to the bound state being closer to threshold and so having a stronger influence. The observed strong coupling of this state to S-wave DK is suggestive of a large molecular DK component – this is not included in the quark model and may account for a significant portion of difference between the observed mass and the expected mass. This study also laid vital groundwork for a more recent preprint [arXiv:2102.04973] investigating the closely related  resonance in S-wave D  scattering.

We also investigated manifestly exotic D  scattering, mapping out the scattering amplitudes for the first time – any states in this channel cannot be composed of solely a quark and an antiquark, and require (at least) two quarks and two antiquarks. While this is difficult to produce experimentally, in the calculation we found attraction in isospin-0 D  and evidence for a virtual bound state (attraction but not strong enough to form a bound state). Further calculations to investigate if an exotic bound state is present in analogous systems with different quark masses will be very interesting.

This work was made possible in part through the DiRAC Data Intensive Service hosted at the University of Cambridge.

---

Blandford-Znajek jets in galaxy formation simulations

Jets launched by active galactic nuclei (AGN) are believed to play a significant role in shaping the properties of galaxies and provide an energetically viable mechanism through which star formation can become 
quenched.

There are numerous open questions surrounding the nature of this jet-mode feedback, however, and due to the complex interactions between the different processes at play in these environments, simulations prove a vital tool to explore these questions. In our work we focus on simulations that can resolve galactic-scale processes but, due to the vast dynamic range that would be required, these simulations cannot resolve the jet launching scales. This means that ab-initio AGN jets will not form and the jet launching processes are necessarily encoded in sub-grid models.

Many such sub-grid models assume that the jets have fixed powers and directions but we know that, in reality, the jet power is modulated by the gas accreted by the black hole and that the amount of gas available depends on the heat that the jet itself is able to impart to its surroundings, linking the properties of the jet with those of the galaxy and forming a feedback loop. One of the most promising processes by which these jets are launched is that of the Blandford-Znajek mechanism. In this model, the jet direction is tied to that of the black hole spin and, thus, is not guaranteed to be fixed. It is therefore of particular importance that self-regulation of AGN jets is fully investigated in the context of galaxy-scale simulations.

Motivated by this, we have developed a self-consistent sub-grid model for AGN feedback in the form of a Blandford-Znajek jet. Our model assumes that the black hole is surrounded by a sub-grid, thin accretion 
disc which modulates the accretion flow onto the black hole, allowing us to accurately follow both the evolution of the mass and angular momentum of the black hole.

Our model has been designed such that it can be employed in galaxy scale simulations with the ultimate aim of assessing whether self-consistent Blandford-Znajek jets can reproduce large-scale galaxy and cluster observables. First, however, we apply our model to the central region of a typical radio-loud Seyfert galaxy embedded in a hot circumgalactic medium (CGM) which allows us to resolve the parsec-scale interactions of jets with the cold dense circumnuclear gas found close to the black hole as well as the interactions of the jet with the surrounding hot CGM.

These simulations were carried out on COSMA7 in Durham and Peta4 in Cambridge as part of the DiRAC project dp012.

We found that jets launched into high pressure environments thermalise efficiently due to the formation of recollimation shocks and the vigorous instabilities that these shocks excite increase the efficiency of the mixing of CGM and jet material (see Fig. 1). The beams of more overpressured jets, however, are not as readily disrupted by instabilities so the majority of the momentum flux at the jet base is retained out to the head, where the jet terminates in a reverse shock. All jets entrain a significant amount of cold circumnuclear disc material which, while energetically insignificant, dominates the lobe mass together with the hot, entrained CGM material. The jet power evolves significantly due to effective self-regulation by the black hole, fed by secularly-driven, intermittent mass flows. The direction of jets launched directly into the circumnuclear disc changes considerably due to effective Bardeen-Petterson torquing. Interestingly, these jets obliterate the innermost regions of the disc and drive large-scale, multi-phase, turbulent, bipolar outflows (see Fig. 2).

---

Hunting for the sources that reionized the Universe

When the first stars emitted light around 100 million years after the Big Bang, their photons began to reionize the Universe bringing it out of the Dark Ages. How long this reionization period lasted, and importantly the main sources providing the photons for reionization are widely debated. Candidates for these sources include dwarf galaxies, massive galaxies, active galactic nuclei (AGN), accretion shocks, stellar mass black holes, globular clusters, and dark matter annihilation and decay.

Ideally one would want to observe these objects directly at the epoch of reionization and simply measure the Lyman Continuum (LyC) radiation that escapes and ionizes the intergalactic medium. However this is near impossible because the photons are absorbed by intervening neutral Hydrogen.  Hence it becomes necessary to turn to numerical simulations to gain insight into how many ionizing photons can emerge from these sources and to see whether reionization epoch galaxies that emit LyC, called LyC leakers, share observational signatures with low redshift analogues that can more easily be observed.

As photons escape by moving through channels cleared out by supernovae and/or AGN, their success at escaping a galaxy in simulations depends on resolving the interstellar medium. In Katz et al. 2020, MNRAS, 498, 164 we report on our work using DiRAC to run a suite of high-resolution Adaptive Mesh Refinement cosmological radiation hydrodynamics resimulations (the Aspen suite) of a massive galaxy (M_DM = 10^11.8 M_sun) at redshift 6.  The simulations include a seven species non-equilibrium chemistry (including H₂) coupled to radiative transfer. We reach resolutions of 1000 M_sun for the stars, 4 X 10⁴ M_sun/h for the dark matter and 9.1pc/h for the gas. We post-processed the simulations with a photoionization code to predict emission luminosities for IR lines and nebular continuum for ~1000 galaxies up to redshift 9.1.

One way that observations try to identify LyC leakers is through observing an excess of [OIII] to [OII]. But our simulations found that such an excess is not a sufficient condition for high escape fractions as only ~11% of our simulated galaxies with an excess of [OIII] to [OII] had LyC escape fractions greater than 10%.  We found that viewing angle (figure 1) as well as ISM physics (figure 2) could cause galaxies with similar [OIII] to [OII] ratios to have very different photon escape fractions.

Therefore, in search of a better diagnostic to identify LyC leaker analogues at low redshift, we used machine learning and found that [CII]_158mm and [OIII]_5007A are the most predictive indicators. As photon escape fractions appear to be regulated by feedback, it is not surprising that lines that probe neutral ([CII]_158mm) and ionized ([OIII]_5007A) gas can pick out LyC leakers.  But since those two lines are unlikely to be observed by the same instrument as one is in the optical whereas the other is in the infrared, we instead tested the power of replacing [OIII]_5007A with the infrared line [OIII]_88mm. We applied our new diagnostics on a sample of local dwarf galaxies to identify LyC leakers and then on epoch of reionization galaxies. The most promising candidate to be actively contributing to the reionization of the Universe was found to be MACS1149_JD1 at z = 9.1.  

These simulations were carried out on DiAL in Leicester as part of DiRAC project dp016.

Some of the Research Highlights from our users in 2019:

• HPQCD: B meson oscillations.
• The external photo-evaporation of planet-forming discs.
• Understanding the fundamental physics of the Earth’s radiation belts.
• Exploring fundamental fields with strong gravity.
• UKMHD: How is the solar corona heated?
• The Simba Cosmological Galaxy Formation Simulations.
• Understanding the Milky Way’s metallicity distribution.
• Supermassive black hole jets stir old cosmic puzzles.
• The effects of ionising star formation on clouds in galactic spiral arms.
• How stars form in the smallest galaxies.
• A steeply-inclined trajectory for the Chicxulub impact.
• Identifying and Quantifying the Role of Magnetic Reconnection in Space Plasma Turbulence.
• Extreme Gravity and Gravitational Waves.
• Spontaneous Symmetry Breaking in the 3d Thirring Model.
• The universal structure of dark matter halos over a mass range of 10²⁰.
• How primordial magnetic fields shrink galaxies.
• Stellar Hydrodynamics, Evolution and Nucleosynthesis.
• The first observational evidence that planets migrate?
• Computational spectroscopy of exoplanets.
• Exploring Jupiter’s Magnetic Field.
• ‘Fast and furious’ planets around tiny stars.
• Monte Carlo radiative transfer for exotic thermonuclear supernovae.
• COSMOS: Studying cosmic strings.
• The evolution of high mass-ratio circumbinary discs.
• Astrophysical Maser Flares from Shock Impact.
• Dust destruction in supernova remnants.
• Gaseous Dynamical Friction in the Early Universe.
• The b1 resonance in coupled-channel scattering from Lattice QCD.
• Electron acceleration and transport in turbulence.
• Extreme QCD: Quantifying the QCD Phase Diagram.
• Deep inelastic structure of the proton.
• New LIGO-Virgo gravitational waves.
• UNITY: Multi-wavelength simulations of relativistic cosmology.
• Impact of Reionization on the Intergalactic Medium.
• Tidal disruption by supermassive black holes.
• DiRAC enables prediction for matter-anti-matter asymmetry in the Standard Model.
• Ion-neutral decoupling in the nonlinear Kelvin—Helmholtz instability: Case of field-aligned flow.
• Simulation of Sp(2N) gauge theories for Composite Higgs models.

2019 list of publications

---

B meson oscillations

Watching the oscillations of two coupled pendulums (for example, attached to the same not-too-rigid support) provides insights into simple but rather counter-intuitive physics. Set one pendulum oscillating and, after a while, the other one will start up as the amplitude of the first dies down. Later still the amplitude of the second pendulum will die down and the first, rather astonishingly, start up again. The energy of the motion sloshes back and forth between the two pendulums because the swinging of either single pendulum is not a simple (‘normal’) mode of oscillation of the coupled system. Setting one pendulum in motion in fact sets both normal modes oscillating and these combine the motion of both pendulums in different ways. As time progresses the interference between the normal modes, which have different frequencies, results in this puzzling (at first sight) behaviour.

An analogous phenomenon occurs in the physics of particles that are bound states of a bottom quark and either an anti-strange or anti-down quark, known as a Bs or Bd meson respectively. We can distinguish these particles from their anti-particles by their decay processes in our particle detectors. However the B and anti-B are coupled together through the weak interaction process shown in Figure 1. The particles with simple time behaviour (the ‘normal modes’) are then mixtures of the B and anti-B. A particle produced as a B may oscillate back and forth into an anti-B (and decay as one in our detector) at a later time with a frequency that can be measured very accurately in experiments at the Large Hadron Collider at CERN. Calculation of the coupling induced by Figure 1 requires the numerical solution of the theory of the strong force (that holds the B and anti-B together) using the method known as lattice QCD. HPQCD has recently done the most accurate calculations of this to date (arXIv:1907.01025) using DiRAC. Comparison to the experimental results enables us to determine the coupling between W bosons, top quarks and strange or down quarks (Vts and Vtd) as shown in Figure 2. Our results agree well with the Standard Model expectation that the CKM matrix of W-boson couplings is unitary, and limit the room for new physics. We also predict the annihilation rate of Bs and Bd mesons to a µ+µ- pair, a key mode for new physics searches at LHC. HPQCD has also done the first lattice QCD calculation of the difference of decay rates of the two normal modes of the Bs (arXiv:1910. 00970).

---

The external photo-evaporation of planet-forming discs

Planets form within flattened discs of material around young stars. These young stars themselves form in clustered groups of up to hundreds of thousands of stars. Most planet-forming discs therefore find themselves in a situation where they are being shone upon by ultraviolet (UV) radiation from these nearby young stars. This has the effect of stripping material from the disc in what is called a photoevaporative wind. This could play an important role in influencing how the discs evolve and the planets that they form.

There are two serious issues that limit our understanding of how external photoevaporation of discs. The first issue is that disc evaporation is very hard to model. To get the temperature correct, which determines the structure of the wind and the rate that mass is lost, one has to solve for the chemistry of the gas and how that is affected by the incoming UV radiation. Furthermore, each cell tries to cool by having “line photons” escape and carry away energy. Each region of the calculation therefore depends on every other region, since if the surroundings are very dense these photons won’t escape and vice versa. We have therefore only had 1D models of the process in the past, except in the strongest UV environments where thee modelling is easier.

T. J. Haworth, as part of the DiRAC project dp100, has developed the TORUS-3DPDR code (Bisbas et al. 2015, Harries et al. 2019), to tackle this problem in 2D models for the first time (Haworth & Clarke 2019). This has vastly improved our understanding of how discs are evaporated, including the rate at which they lose mass and where the mass is actually lost from the disc.

The other key challenge to understanding of how external photoevaporation is observing it in action. 1D models are insufficient for this purpose, but with these new 2D models, predictions can be made as to what should be observed (Haworth & Owen 2020)

References

• Bisbas et al. (2015), MNRAS, 454, 2828
• Harries et al. (2019), A&C, 27, 63
• Haworth & Clarke (2019), MNRAS, 485, 3895
• Haworth & Owen (2020), MNRAS, 492, 5030

---

Understanding the fundamental physics of the Earth’s radiation belts

Electromagnetic waves interact strongly with charged particles in the Earth’s inner magnetosphere. It is important to be able to model the evolution of these particles, since we rely on the many satellites that orbit within this hazardous radiation environment. Our current approaches to modelling the effect of wave-particle interactions in the outer radiation belt, over macroscopic time and length-scales, rely heavily upon quasilinear wave theory (e.g. see [1]). This theory describes wave-particle interactions as diffusive processes in the plasma, flattening out gradients and slowly moving a small fraction of electrons to higher energies (leading to acceleration), or different pitch-angles (leading to loss). However, the quasilinear formalism relies upon a number of assumptions that may not hold for all important wave types in the radiation belts. Use of the DiRAC facilities has enabled us to carry out important studies into the efficacy of the quasilinear formalism.

In [2], we used the EPOCH particle-in-cell code (see [3]) to conduct numerical experiments to investigate electron interactions with an incoherent spectrum of whistler-mode waves. Our novel approach directly extracts diffusive characteristics across all energy and pitch angle space. This benchmarking work establishes a framework for future investigations on the nature of diffusion due to whistler-mode wave-particle interactions, using particle-in-cell numerical codes with driven waves as boundary value problems.

In [4], we use the techniques developed in [2] to investigate the electron response to whistler-mode waves as a function of increasing wave amplitude. We find that whistler-mode waves with amplitudes of order (dB/B)² ~ O(10^-10) – O(10^-6) in a uniform B give diffusive and advective dynamics. Over appropriately short timescales the diffusive component of the dynamics agrees with quasilinear theory even for the highest amplitude waves. These timescales range from thousands to tens of gyroperiods for the lowest and highest amplitude waves respectively. This provides very helpful insight into the fundamental nature of electron dynamics in energy and pitch-angle space due to the wave-particle interaction under the uniform background-field assumption, used in the basic quasilinear theory. This also motivates a variety of future numerical experiments and theoretical investigations on the applicability (or otherwise) of quasilinear diffusion theory to electron dynamics in the Earth’s radiation belts, using the particle-in-cell method established in [2].

References

• [1] Glauert, S. A., and Horne, R. B. (2005), doi:10.1029/2004JA010851.
• [2] Allanson, O. et al. (2019), doi:10.1029/2019JA027088
• [3] Arber, T. D. et al. (2015), doi:10.1088/0741-3335/57/11/113001
• [4] Allanson, O. et al (2020), doi:10.1029/2020JA027949

---

Exploring fundamental fields with strong gravity

Einstein’s theory of general relativity is both complex and beautiful; even after 100 years, there is still much it can reveal to us about the nature of gravity. It may also help us answer other fundamental questions, such the nature of the dark matter component of the universe. So far, this component is only known to interact gravitationally with other matter, and so its behaviour in strong gravity regimes may be the only way to probe key properties, such as its mass.

This project used numerical simulations of fields in gravitational backgrounds, to explore the behaviour of light dark matter particles around black holes. We demonstrated how accretion of light dark matter onto a black hole leaves a distinctive pattern in the density (above left) and began to study how its behaviour around spinning black holes might probe its self interactions (below right). Future work will extend these studies to the binary case, where such dark matter imprints may be observable by future detectors.

---

UKMHD: How is the solar corona heated?

The only viable mechanism to heat the outer solar atmosphere must involve the magnetic field. Two main approaches have been adopted. One is to model the heating along individual magnetic field lines and, by specifying the spatial and temporal form of the heating term in the thermal energy equation, try to compare predicted temperatures and densities with the observed values. In this way, constraints on the form of the heating are determined. However, the heating used is not determined from physical processes. The advantage of this approach is that the rapid spatial variations in the transition region can be adequately in 1D. The other approach is to solve the full MHD equations, without any assumptions on the form of the heating and see how heating occurs due to Ohmic and viscous effects. This second approach, while physically more realistic, requires the numerical resolution in 3D of large but narrow current sheets in which magnetic reconnection can occur. In the nanoflare heating model of Parker, a large number of intense current layers are required and this requires extremely large numerical grids to resolve the sharp current gradients. Hence, Dirac is the only UK computing facility that can cope with the computational requirements. However, it remains extremely difficult to resolve the transition region fully when the thermodynamics is included.

In a unique set of numerical experiments, the evolution of the solar coronal magnetic field is modelling by slowly stressing it, through slow photospheric motions. Once a significant amount of free magnetic energy has built up, and instability triggers a large release of energy. The important point is that the magnetic field does not release all the stored energy but tries to relax towards an non-potential state. After this large event, the imposed photospheric motions continue and the magnetic field releases more energy in the form of heat and kinetic energy in a series of smaller but more frequent events. The energy release is spread across the whole region covered by the magnetic field and this is known as an MHD avalanche process. For the first time, Reid et al, A&A 633, 2020, have used Dirac to investigate in detail how the heating process occurs and where the large number of heating are spatially located. Hence, the form of the coronal heating function is determined, without any prior assumptions, and can be used with the field aligned modelling to determine the plasma temperature and density for comparison with observations.

Figure 1 (left) shows a contour plot of the heating along one field line as a function of distances along the field and time. Red indicates strong heating events. They are localised in space and time. Figure 1 (middle) shows the time averaged heating along the same field line. There is one large heating event about 0.4 of the distance along the field line but this is due to the very first event. After that there is no preferred location for the heating. Figure 1 (right) shows the spatially averaged heating as a function of time. Again, there is one large event that start off the avalanche process and then the heating comes in regular bursts. For more details see Reid et al Astron. & Astrophys., 633, id.A158.

---

The Simba Cosmological Galaxy Formation Simulations

Creating accurate models of our Universe is one of the great quests of humankind, going back to the ancient Greeks, through to Copernicus, Galileo, and Einstein. In recent years, huge progress has been made in developing physical models that connect the early conditions we see shortly after the Big Bang with the present-day population of galaxies, clusters, and intergalactic gas that we see around us today. Owing to the complex and highly nonlinear nature of gravitational collapse and radiative processes, supercomputer simulations play an increasingly critical role in developing models to explain why the Universe looks the way it does.

The Simba simulations are a next-generation suite of galaxy formation models, run on Cosma-7 at Durham. Simba includes unique physical modules that address some of the most pressing questions in modern galaxy formation theory, such as the growth of supermassive black holes and their impact on their host galaxy, and the formation and evolution of dust grains within a cosmological context. Simba predictions are in remarkably good agreement with a wide range of galaxy demographic properties, and have elucidated new ways in which galaxies and black holes interact with their surrounding cosmic ecosystems to shape the observable properties of both. Simba represents our most comprehensive and faithful view yet of the so-called baryon cycle that drives the formation and evolution of galaxies like our own Milky Way over cosmic time.

Simba has had a major impact on the community since its release in July 2019, with 18 publications as of June 2020 based primarily on Simba data, and over 30 active projects currently listed on the Simba pbworks page. The first annual Simba team meeting was held in February 2020, with over 30 attendees representing 17 institutions from 5 countries, demonstrating the global impact of Simba. Simba is gearing up for a public data release in the Fall, which will further widen the user base and establish it among the pre-eminent cosmological galaxy formation simulations today.

---

Understanding the Milky Way’s metallicity distribution

An important aspect of understanding the formation of bulges of galaxies, including that of the Milky Way, is understanding their metallicity trends. The most notable feature of the metallicity distribution is that when the bulge is viewed edge-on, it appears more vertically pinched than the density. In the Milky Way this manifests as an X-shape when viewed from the Solar perspective. Simulations that include gas, star formation and the chemistry of stars have managed to reproduce the pinched appearance of the metallicity (Debattista et al. 2017). However these simulations are expensive and do not permit an easy detailed matching to the Milky Way. In the past this has forced researchers to paint metallicity on particles in N-body simulations (i.e. without any gas or star formation). This type of painting has assigns a metallicity based on the position of a star (both its distance from the galactic centre and its height above the mid-plane). However this is equivalent to assuming that all stars at a particular position are born at the same time with the same metallicity, which is highly implausible. Using a suite of N-body simulations, we have shown that where a star ends up on a bulge depends strongly on its initial dynamical actions, which measure its angular momentum as well as its radial and vertical random motion. This allowed us to develop a novel technique for painting metallicity onto star particles using the dynamical actions. We showed that this technique produces metallicity maps comparable to those of star forming simulations. This more realistic metallicity assignment allows us to much more rapidly model the formation of the bulge of the Milky Way. The figure shows an example of such a metallicity mapping: the metallicity distribution (colours and black contours) is more pinched along the minor axis than is the density (white contours), as observed in the Milky Way and in external galaxies.

Moreover, using the actions to understand bulge formation has lead us to a number of important new insights. We showed for instance that the vertical gradient of the vertical action is largely erased by the formation of the bulge, whereas the very weak gradient of the radial action is vastly enhanced during the same time. This means that the vertical gradients observed in galaxies are not the result of the superposition of a metal-rich thin disc and a metal-poor thick disc. Another, quite surprising, result of our use of dynamical actions is that there are populations that vertically cool during bulge formation. We are investigating the properties of these populations and hope to be able to use them to constrain what the Milky Way looked like before the bulge formed.

---

Supermassive black hole jets stir old cosmic puzzles

Galaxy clusters are the most massive gravitationally bound systems and encode unique information on the composition of our Universe. They comprise of a massive dark matter halo and up to a few thousand rapidly moving galaxies enveloped by a very hot X-ray emitting plasma, known as the intracluster medium (ICM). Of particular importance for regulating the heating and cooling of the ICM, the decades long “cooling flow” problem, are the relativistic jets produced by supermassive black holes with masses in excess of a billion solar masses. Notwithstanding massive theoretical and observational effort, a major puzzle left to crack is how the energy from these powerful jets is transferred to the ICM to stop the catastrophic cooling.

Using the moving-mesh code AREPO, we performed state-of-the-art simulations with the highest-resolution jets to-date within a large scale, fully cosmological cluster environment to shed light on this issue. The simulations make use of a novel refinement technique (see Figure 1) that not only allows high resolution close to the central black hole but also within the jet lobes themselves. We found that the mock X-ray observations of the simulated cluster revealed the so-called “X-ray cavities” and “X-ray bright rims” generated by supermassive black hole-driven jet, which itself is distorted by motions in the cluster remarkably resemble those found in observations of real galaxy clusters (see Figure 2). While it is well accepted that the vast amount of energy injected into the jet lobes (which would be enough to meet the Earth’s power usage for more than 10³¹ years) is sufficient to offset the cooling losses within the ICM and prevent the formation of a cooling flow, how this energy is effectively and isotopically communicated to the ICM is still an open debate. Our simulations showed that the ICM motions or “weather” are crucial to solving this problem. Through the cluster weather action, the jet lobes can be significantly moved and deformed, stirring them around the cluster, which ultimately leads to jet lobe disruption and effective energy transfer to the ICM which could be the long sought-after solution to combat the “cooling flow” puzzle.

This work has been published in the Monthly Notices of the Royal Astronomical Society: “AGN jet feedback on a moving mesh: lobe energetics and X-ray properties in a realistic cluster environment” by Martin A. Bourne, Debora Sijacki and Ewald Puchwein. Royal Astronomical Society press release: “Stormy cluster weather could unleash black hole power and explain lack of cosmic cooling”, 14 October, 2019.

---

The effects of ionising star formation on clouds in galactic spiral arms

The formation and life-cycle of giant molecular clouds is a crucial part of understanding both the evolution of galaxies and the varying environments in which the star formation process occurs. The origins of the physical phenomena that govern giant molecular clouds and the star formation within them include external effects caused by the galactic environment (such as galactic potentials and external radiation fields) as well as more local effects such as stellar feedback and turbulence. Recently, we performed simulations of a galactic spiral arm consisting of complexes of giant molecular clouds that included ionising radiation from young stars. Unlike previous work in this area, which typically simulated isolated clouds, our models include a number of clouds along a spiral arm, and thus ionisation fronts propagate in a much more realistic environment. The initial conditions were achieved by extracting a 500 pc² section of spiral arm from a galaxy scale simulation by Dobbs & Pringle (2013) and increasing the resolution.

Using the simulations we examine the effect of photoionising feedback caused by the most massive stars forming within the clouds. Lyman continuum photons emitted from these massive stars ionise the surrounding gas creating bubbles – HII regions – which are bounded by expanding shocks. Photoionisation has the highest energy budget of the early onset stellar feedback mechanisms. Comparison of simulations with and without photoionisation shows that the timing and location at which stars form is significantly affected (Fig. 1). While the final fraction of gas converted to stars remains relatively unchanged, the time taken is reduced by a factor of nearly 2. This is due to star formation in nearby dense gas being compressed by shocks from the surrounding ionised regions. The result is an acceleration of the star formation rate and a much wider distribution of star formation across the region.

---

How stars form in the smallest galaxies

It has been a long-standing puzzle how the smallest galaxies in the Universe have managed to continuously form stars at such a remarkably low rate – tiny galaxies like Leo P (pictured above) form approximately one star every 100,000 years! Making use of the DIRAC DIaL supercomputer, our research team has established how such low-levels of star formation arise because small dwarfs can stay quiescent for several billions of years, slowly accumulating gas in their centres. When this gas eventually collapses, star formation resumes at rates consistent with observations.

Through high-resolution computer simulations, we demonstrate that star formation in the very smallest dwarf galaxies is shut down in the early Universe as a result of heating and ionisation from the strong light of new-born stars, exploding stars and stellar winds. However, we find that this quiescent state can be temporary. Some of our simulated dwarf galaxies are able to slowly accumulate fresh fuel in the form of hydrogen gas, over many billions of years. Once a sufficient amount of gas has piled-up to overcome any residual heating, star-formation re-ignites and proceeds continuously until today. This scenario explains, for the first time, the puzzlingly low star formation rates observed in faint dwarf galaxies such as Leo P, Leo T, and others. The computational power provided by DIRAC DIaL was crucial to establish this result, allowing us to afford detailed modelling of the star formation and heating processes in these tiny galaxies, over the full lifetime of the Universe.

By deepening our understanding of dwarf galaxies, we gain new insights into how astrophysical processes, such as stellar explosions and the heating and cooling of cosmic gas, regulate galaxy formation. This work also opens the door to harnessing these galaxies as a new test of the “Big Bang” cosmological model. Further work is currently underway to predict how many such star-forming dwarfs exist in our Universe and how many will be discovered by the next generation of astronomical telescopes.

---

A steeply-inclined trajectory for the Chicxulub impact

Simulations performed on the DiRAC High Performance Computing Facility have revealed the asteroid that doomed the dinosaurs struck Earth at the ‘deadliest possible’ angle.

The simulations reveal that the asteroid hit Earth at an angle of about 60 degrees, leaving a signature of its trajectory in the rocks beneath the Chicxulub Crater in Mexico.

Such a strike would have maximised the amount of climate-changing gases thrust into the upper atmosphere, including billions of tonnes of sulphur that blocked the sun and triggered the KPg mass extinction. It is thought that up to 75% of plant and animal species on Earth became extinct, including dinosaurs.

The new simulations are the first ever fully 3D numerical calculations to reproduce the whole impact event, from the moment of impact to the moment the final Chicxulub crater was formed. Results of the simulations were compared with geophysical data from the site of the impact crater to diagnose the direction and angle of impact.

The results from these simulations were published in Nature Communications (Collins et al., 2020), and attracted widespread international media attention [2, 3]. The article was also featured on the STFC website [4].

• [1] Collins, G.S., Patel, N., Davison, T.M. et al. A steeply-inclined trajectory for the Chicxulub impact. Nat Commun 11, 1480 (2020). https://doi.org/10.1038/s41467-020-15269-x
• [2] https://www.altmetric.com/details/82771548
• [3] https://www.imperial.ac.uk/news/197694/dinosaur-dooming-asteroid-struck-earth-deadliest-possible/
• [4] https://stfc.ukri.org/news/uk-supercomputers-simulate-impact-of-the-asteroid-thought-to-have-wiped-out-dinosaurs/

---

Identifying and Quantifying the Role of Magnetic Reconnection in Space Plasma Turbulence

Plasma is the fourth state of matter after the solid, liquid, and gaseous states. Almost all of the visible matter in the Universe is in the plasma state. This includes all the stars, the material between the stars, the intergalactic medium, and also the gas in the Earth’s cosmic neighbourhood: the solar wind and the geospace environment in Earth’s magnetosphere. Very often, these plasmas have a low viscosity. Therefore, we expect that these plasmas are highly turbulent. The goal of our project is to understand the role of a process called “magnetic reconnection” in the omnipresent plasma turbulence in space.

Magnetic reconnection is a multiscale plasma phenomenon in which the magnetic field re-organises and transfers energy into the plasma particles. It occurs in laboratory plasmas, such as fusion or confinement experiments, and in astrophysical plasmas, such as the solar wind, coronal mass ejections, solar flares, explosive events in planetary magnetospheres, accretion discs, and star-formation regions. Although reconnection has been studied for over 60 years, there is still no consensus about a complete theory to describe magnetic reconnection at all scales involved.

Current theories and numerical simulations suggest that turbulence and magnetic reconnection are closely related. Turbulence transfers energy across the scales of the plasma. Energy is injected at large spatial scales and then transported to the very small scales where it dissipates as heat. This transfer is known as the turbulent energy cascade. It is proposed that the ultimate dissipation of energy at small scales strongly depends on a link between turbulence and reconnection.

We study this link between turbulence and reconnection with large computer simulations of plasma turbulence. The computational power provided by high-performance computing facilities at DiRAC enables us to study plasma processes from first principles using three-dimensional cutting-edge particle-in-cell (PIC) simulations. This type of simulation resolves the smallest plasma scales and accounts for phenomena that only reveal themselves in kinetic theory. Even with the enormous power of DiRAC, full PIC simulations are unable to cover the whole range of scales involved in plasma turbulence and reconnection in a real plasma, since these simulations are very expensive in terms of computing memory and time. With DiRAC’s capabilities, however, we have taken a big step forward to understanding turbulence and magnetic reconnection in plasmas.

We initialise our computations 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. We find that, after some time, the turbulence itself generates current sheets and forms regions of magnetic reconnection. We define a set of criteria to identify where in the simulated plasma reconnection occurs. With these criteria, we find one extensive reconnection site in our simulation and study its properties in great detail.

One huge 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 space plasmas like the solar wind. The measurements from our artificial spacecraft will serve as predictions for measurements of reconnection in the turbulent solar wind by ESA’s and NASA’s latest solar-system spacecraft Solar Orbiter and Parker Solar Probe which have the appropriate instruments onboard to test our predictions.

---

Extreme Gravity and Gravitational Waves

When two black holes inspiral and merge, they emit gravitational waves that have famously been detected for the first time by LIGO in 2015 [1]. These gravitational waves carry energy as well as momentum away from the black-hole binary system. The energy loss is responsible for the inspiral and eventual merger. The linear momentum radiated away, on the other hand, can impart a recoil on the black-hole system, just like the firing of a bullet imparts a recoil on the shooter. This effect follows from conservation of momentum; for black holes, the gravitational waves carry net momentum in some direction, so the black-hole emitter has to respond by moving accordingly in the opposite direction.

One of the most dramatic results of numerical relativity has been the discovery of “superkicks”, very large recoil velocities of up to 3,700 km/s that occur when the black holes start their inspiral with specific spin directions [2,3]. Such large kicks are sufficient to eject black holes from even the most massive host galaxies; the escape velocities from giant elliptic galaxies are about 1,000 km/s. This has sparked a lot of interest among astrophysical observers to look for ejected black holes, because this would have major implications for the formation history of supermassive black holes. Several candidates have been found but alternative explanations (without resorting to kicks) cannot be ruled out; see [4] for a review.

As part of our DiRAC project, we have found that even larger recoil velocities, up to 4,300 km/s are realized when the black holes inspiral on moderately eccentric (rather than quasi-circular) orbits. In general, moderate eccentricity can amplify the recoil effect by up to 25% with possible effects on the retention rate of black holes in globular clusters, second-generation populations of black-hole merger events, as well as the black-hole occupation fraction of galaxies.

• [1] B. P. Abbott et al, Phys.Rev.Lett. 116 (2016) 061102, arXiv:1602.03837 [gr-qc]
• [2] J. A. Gonzalez et al, Phys.Rev.Lett. 98 (2007) 231101, gr-qc/0702052
• [3] M. Campanelli et al, Astrophys.J. 659 (2007) L5-L8, gr-qc/0701164
• [4] S. Komossa, Adv.Astron. 2012, 364073, arXiv:1202.1977 [astro-ph]
• [5] U. Sperhake et al, Phys.Rev.D 101 (2020) 024044 arXiv:1910.01598 [gr-qc]

---

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 v_F/c ≈ 1/300, where the pitch of the cones v_F 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 N_c: determination of N_c 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 L_s to ∞. Simulations have been performed using the RHMC algorithm on 12³ and 16³ systems with L_s ranging from 8 to 48. It turns out that at strong couplings g² and light fermion masses m recovery of GW symmetry is slow – nonetheless at weak-to-moderate couplings our results for <ψψ> as a function of g², m obtained at L_s = 48 are compatible with an extrapolation to the large L_s 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/g² ≈ 0.28 (see figure). The fit is less convincing at the strongest couplings (lying within the broken phase), which we attribute to the L_s to ∞ 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 N_c>1. This is a significant step towards understanding the formulation of fermion quantum fields in strongly-interacting theories.

---

The universal structure of dark matter halos over a mass range of 10²⁰

Cosmological models in which the dark matter consists of cold elementary particles predict that the population of dark matter haloes in the Universe should extend to masses comparable to the Earth’s, many orders of magnitude below those where galaxies can form. Using a novel multi-zoom technique, we carried out a consistent cosmological simulation of the formation of present-day haloes over the full mass range populated when dark matter is assumed to be a Weakly Interacting Massive Particle (WIMP) of mass ~100 GeV. The simulation has a dynamic range of 30 orders of magnitude in mass and resolves the internal structure of hundreds of Earth-mass haloes in as much detail as that of hundreds of rich galaxy clusters. We find that halo density profiles are universal over the entire mass range and are well described by simple two-parameter fitting formulae, such as the well-known Navarro-Frenk-White profile. Halo mass and concentration are tightly related in a way that depends on cosmology and on the nature of the dark matter. At fixed mass, concentration is independent of local environment for haloes less massive than those of typical galaxies. Haloes over the mass range (10³ – 10¹¹) M_o contribute about equally (per logarithmic interval) to the dark matter annihilation luminosity, which we find to be smaller than all previous estimates by factors ranging up to one thousand. This was a hugely challenging project that took seven years to complete, during which we uncovered and solved a number of technical problems in N-body simulations revealed only in the extreme conditions of our calculation. It was carried out on DiRAC’s Cosmology Machine (COSMA) at Durham.

References:

• Jie Wang, Sownak Bose, Carlos Frenk, Adrian Jenkins, Liang Gao, Volker Springel, and Simon White, to appear in Nature.

---

How primordial magnetic fields shrink galaxies

Little is known about how primordial magnetic fields are generated nor about their initial configurations. Their exact primordial normalization, coherence length and spectral index are undetermined but they are widely believed to be many orders of magnitude lower in strength (3 x 10^-20 G < B₀ < 3 x 10^-10 G) than their microGauss values in present-day galaxies. What is known is that magnetic fields play as large a role in galaxies as turbulence, thermal and cosmic ray energy densities and that therefore they have a crucial role to play in shaping a galaxy’s evolution. To understand how different configurations of primordial magnetic fields might affect the global morphological and dynamical properties of galaxies, we used DiRAC to perform a suite of high-resolution constrained transport magneto-hydrodynamic cosmological zoom simulations where we varied the initial magnetic field strength and configuration along with the prescription for stellar feedback. In Martin-Alvarez et al, 2020, MNRAS, 495, 4475 we report our findings that strong primordial magnetic fields delay the onset of star formation and drain the rotational support of the galaxy, diminishing the radial size of the galactic disk and driving a higher amount of gas towards the centre. A possible mechanism behind such a reduction in angular momentum is magnetic braking. The figure above shows the effect of increasing primordial magnetic field strength on the stars and gas of a galaxy from our simulations. For the strongest primordial magnetic field studied (B₀ = 3 x 10^-10 G ), the gas radial scale length is halved compared to the simulation where the primordial magnetic field is weakest. New Scientist magazine reported on our work: https://www.newscientist.com/article/2244840-the-milky-way-may-have-been-shrunk-down-by-ancient-magnetic-fields/

---

Stellar Hydrodynamics, Evolution and Nucleosynthesis

The stellar hydrodynamic group of Prof Raphael Hirschi at Keele University and their international collaborators computed a series of 3D hydrodynamical simulations of stellar convection during carbon burning (see figure 1). These 3D implicit large eddy simulations were computed with the PROMPI code on the COSMA supercomputer in Durham. The simulations were analysed on-the-fly with PROMPI’s built-in mean-field analysis based on the Reynolds-averaged Navier–Stokes framework (RA-ILES, see Arnett et al 2019 and references therein). Both the vertical convective velocities within the convective region and the bulk Richardson number of the boundaries were found to scale with the driving luminosity as expected from turbulence theory: v ∝ L^1/3 and Ri_B ∝ L^-2/3, respectively. The positions of the convective boundaries were estimated through the composition profiles across them, and the strength of convective boundary mixing was determined by analysing the boundaries within the framework of the entrainment law. The entrainment was found to be approximately inversely proportional to the bulk Richardson number, Ri_B (∝ Ri_B^-a, with a ~0.75). Although the entrainment law does not encompass all the processes occurring at boundaries, the results support the use of the entrainment law to describe convective boundary mixing in 1D models and new 1D models including entrainment show promising results. These large-scale simulations also inform the team’s ongoing theoretical efforts to overcome the shortcomings of the widely used but outdated mixing-length theory (Arnett et al 2019).

---

The first observational evidence that planets
migrate?

The fact that young planets migrate as a result of interaction with their natal disc is well established theoretically and is widely applied when explaining the orbital properties of exoplanets. To date however the very long timescales for planetary migration have meant that there is no direct evidence for this. Moreover theoretical estimates for the rate of planetary migration are rather uncertain.

We have proposed the first potential observational test of planetary migration, pointing out, on the basis of our hydrodynamical simulations of dust and gas in discs, that a migrating planet has a different signature from a stationary planet in terms of the way that it sculpts the disc dust. We have modeled the case of the young disc Elias 24 where the Atacama Large Millimetre Array (ALMA) has demonstrated a pronounced gap (presumed to contain a planet) and a bright emission ring outside. We have shown that the dust emission profile of the disc is better fit by models where the planet migrates . In order to test this hypothesis further we are seeking to acquire data at a range of wavelengths that are sensitive to a range of dust grain sizes. We predict that smaller grains (which are probed by shorter wavelength observations) should be preferentially concentrated in a ring interior to the planet whereas larger grains should instead congregate outside the planet; this therefore changes the emissivity profile at different wavelengths. The predicted profiles depend on how fast dust grains drift inwards relative to how fast the planet is migrating. Thus multi-wavelength observations provide a way to estimate the planet’s migration rate.

---

Computational spectroscopy of exoplanets

ince the first discovery of a planet orbiting a distant star 25 years ago (Nobel Prize in Physics 2019) over 4000 of exoplanets have been discovered. Now we want to know what these exoplanets are made of, what weather is like there, and whether they are habitable. These questions can be answered via characterisation of their atmospheres, which requires detailed models of the planetary atmosphere and large quantities of laboratory data on how possible atmospheric species absorb light.

In 2019, a UCL team comprising Dr Angelos Tsiaras, Dr Ingo Waldmann, Prof Giovanna Tinetti, Prof Jonathan Tennyson, and Prof Sergey Yurchenko made the first detection of water in the atmosphere of an exoplanet, called K2-18b orbiting within the star’s habitable zone. This publication in Nature Astronomy, 3, 1086 (2019) was the highest ranked astronomy paper of 2019 with some 4000 worldwide media reports.

An atmospheric retrieval requires significant input on the spectral properties of the individual atoms and molecules which comprise this atmosphere. The ExoMol group led by Profs. Jonathan Tennyson and Sergey Yurchenko has pioneered quantum mechanical-based theoretical procedures which allow them to compute the appropriate lists of absorption lines. The ExoMol database provides comprehensive and accurate opacities for the key molecules which are likely to exist in the atmospheres of exoplanets. ExoMol line lists are used by exoplanetary models world-wide, to model the behaviour of molecules in the atmospheres of exoplanets. Because the atmospheres of most observable exoplanets are rather hot (often well over 1000 K) the molecular line list become huge. The calculation of such big data can be only accomplished on modern HPC facilities and DiRAC has been a reliable partner for ExoMol in providing these facilities since 2012.

ExoMol database contains extensive line lists for 80 molecules, 190 isotopologues. Most its 700 billion transitions have been produced using DiRAC resources, with 10 line lists resulted from the new thematic project “Spectroscopy of hot exoplanets”. One of these data is the line lists for SiO₂ consisting of over 32 billion transitions (Owens et al., MNRAS, 495, (2020)). SiO₂ is believed to be an important absorber in atmospheres of hot “lava planets” which orbit so close their host star that their rocky surface must be molten. Until our calculations, there existed no laboratory data on the absorption properties of this key molecule.

Owing to the success of ExoMol and the demand for these data ERC have just funded a follow-up project, ExoMolHD, led by Prof Tennyson, on precision spectroscopic data for studies of exoplanets and other hot atmospheres which will run from 2020 to 2025.

---

Exploring Jupiter’s Magnetic Field

The Juno spacecraft reached Jupiter in 2016 and has since been taking measurements of the gas giant’s magnetic and gravitational fields. Many interesting features on Jupiter’s surface, such as the intense isolated magnetic flux patch near the equator dubbed the “Great Blue Spot”, have been observed. Juno’s measurement also extends the Lowes spectrum—magnetic energy per spherical harmonic degree at the planet’s surface— considerably further than old data allowed. Using our numerical model of Jupiter’s dynamo, we investigate the connection between the shape of magnetic energy spectrum and the dynamo radius—the depth below the planetary surface at which dynamo action begins (Yue-Kin Tsang and Chris Jones, Earth Planet. Sci. Lett., vol. 530, 115879, Jan 2020).

Our dynamo model solves the MHD equations describing the behaviour of an electrically conducting fluid in a spherical shell driven by convection. A key point is that the electrical conductivity increases sharply but smoothly with depth so it is not known a priori where the dynamo radius is. The density of Jupiter varies greatly with depth, so an anelastic model is required to represent the convection accurately. Furthermore, the model parameters needs to be in the strong dynamo regime where the magnetic field strongly affects the flow. Both features are essential for realistic Jupiter models, but they increase computational cost significantly, making the use of the DiRAC facilities essential for our work.

The numerical model produces Jupiter-like dipolar dominant magnetic field. The Lowes spectrum at the surface, as well as the magnetic energy spectrum at different depths, are calculated. We find that within the region where magnetic field is being generated, the magnetic energy spectrum is shallow and maintains more or less the same shape. Outside this dynamo region, the spectrum steepens. This transition enables us to identify the dynamo radius in the model. We can compare our numerical spectrum with the spectrum observed by Juno, which shows that the dynamo radius of Jupiter cannot be deeper than 0.83 Jupiter radius.

Although Juno has revealed the intense flux patches predicted by our simulations, such as the great blue spot, the observed field is nevertheless smoother than expected at the high magnetic Reynolds numbers inside Jupiter. One possible explanation is the existence of a stably stratified layer just under the upper non-conducting molecular layer due to “helium rain-out”, high pressure making the helium in the gas giant form droplets which fall under gravity. This stable layer might also help explain the fierce zonal winds seen on Jupiter. We are currently including a stably stratified layer in our models, to see if this can give a dynamically consistent picture fitting the Juno data.

---

‘Fast and furious’ planets around tiny stars

New astronomy research from the University of Central Lancashire suggests giant planets could form around small stars much faster than previously thought. As published in the “Astronomy and Astrophysics Journal” (https://www.aanda.org/10.1051/0004-6361/201936954), Dr Anthony Mercer and Dr Dimitris Stamatellos’ new planet formation research challenges our understanding of planet formation.

Red dwarfs, the most common type of stars in our Galaxy, are small stars, 10% to 50% the size of our Sun. Despite their small mass, they are found to host giant planets up to 10 times bigger than Jupiter, the largest planet in our Solar System. The formation mechanism of these big planets remains an unsolved mystery. Giant planets around stars, like our Sun, are thought to have formed by the gradual build-up of dust particles to progressively bigger bodies. However, red dwarfs are tiny when compared to the Sun, and they do not seem to have enough material around them to form such big planets.

The research team used the UK Distributed Research using Advanced Computing (DiRAC) supercomputing facility to simulate the evolution of protoplanetary discs around red dwarf stars. Protoplanetary discs are rotating structures of dense gas and dust found around all newly-born stars. The researchers found that if these young discs are big enough they can fragment, i.e. break up into pieces, forming gas giant planets. This theory predicts that the formation of giant planets happens within a few thousand years, a timescale which is extremely fast in astrophysical terms. The researchers also found these planets are extremely hot when they form, with temperatures at their cores reaching thousands of degrees. Such hot planets would be relatively easy to observe when they are still young. They do not have an internal energy source, so they become fainter with time, and the window of opportunity to directly observe them is very small. Nevertheless, they can still be indirectly observed by their effect on their host star. Many planets like these have been discovered so far around small stars. Future observations of planets around very young red dwarf stars will test the predictions of this new theory.

A high resolution animation can be downloaded at here, and a related YouTube animation can be found at here.

---

Monte Carlo radiative transfer for exotic thermonuclear supernovae

The last few years have seen many exciting developments in time-domain astronomy. In particular, thanks to modern surveys that monitor large areas of the night sky on a regular basis, there has been a huge increase in the both the quantity and quality of observational data for a wide variety of astrophysical transients, including supernovae. Such data allow us to test our understanding of astrophysical phenomena in ways, and to a degree of accuracy, that was not previously possible. For example, while the basic picture that supernovae result from the explosive deaths of stars is well established, it has become clear that there is a great deal of variation within the supernova population. This applies even to classes of supernovae that have been traditionally considered rather homogenous, such as the thermonuclear “Type Ia” supernovae (Taubenberger 2017).

Our work focuses on using computer simulations to make predictions for the radiation emitted by supernovae and related stellar explosions. By comparing our results to observations, we hope to test the underlying theoretical models and to help interpret new observational data. Thus, driven by observational advances, a major component of our work in recent years has been understanding how the observed diversity can be reconciled with particular explosion scenarios.

In 2020 we published new results based on simulations of a particular class of white-dwarf explosion and showed that these can provide a plausible match to certain unusual supernovae. Specifically, our work explores the “double-detonation” model in which a white dwarf star, which is composed mainly of carbon and oxygen, is able to accumulate a surface layer of helium via interaction with a companion star in a binary system. If nuclear reactions can ignite in the helium layer, it has been suggested that this may lead first to detonation of the helium which, in turn, may trigger detonation of the underlying carbon-oxygen material. This “double-detonation” scenario predicts certain key signatures owing to the structure where nuclear ash from the carbon-oxygen core is overlaid by material that originated in the helium layer. Our simulations show (see Figure; details published in Gronow et al. 2020) that the theoretical spectra predicted from such a model give a good match to particularly unusual transients, such as SN 2016jhr (Jiang et al. 2017). Such comparisons support the conclusion that transients of this sort may indeed be explainable via this class of theoretical model – one piece for the puzzle of understanding the origins of supernovae. We are now developing further such simulations to explore more fully the range of properties that double-detonation models may be able to explain – it remains to be seen whether this class of model can really account for a large fraction of the observed Type Ia supernova population or whether a mix of scenarios must be invoked to account for the variations that are observed.

This work also used the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

References:

• Taubenberger S., 2017, in Handbook of Supernovae, ed. Alsabti & Murdin (New York Springer), 317;
• Jiang J., Doi M., Maeda K., et al. 2017, Nature, 550, 80;
• Gronow S., Collins C., Ohlmann S., Pakmor R., Kromer M., Seitenzahl I., Sim S., Röpke F., A&A, 2020, 635, 169.

---

COSMOS: Studying cosmic strings

COSMOS consortium members (dp002) have made good progress over the past year using DiRAC HPC Facilities, particularly implementing adaptive mesh refinement (AMR) techniques for the study of extreme gravity in a cosmological context. This year we highlight high resolution simulations of radiation from cosmic strings which were generated using GRChombo, that is, the pioneering numerical relativity code developed within the COSMOS consortium and now attracting a growing research community (see the public repository github.com/GRChombo). The improving sensitivity of the LIGO/Virgo gravitational wave detectors ultimately motivated this work because of the discovery potential of these experiments for cosmic strings, but the initial applications are actually on the origin of dark matter axions.

Cosmic strings – heavy line-like defects traversing through space – can arise when the rapid cooling of the hot early Universe triggers a phase transition. For this reason, cosmic strings are strongly tied to fundamental physics and so their detection could offer deep insights into the underlying structure of particle physics and spacetime itself. Accurate signal templates are essential for observational searches, but numerical simulations of cosmic strings are computationally challenging because of the wide range of physical scales involved; the ratio between the string width and the scale of the Universe across which the strings stretch is approximately 10³⁰! This is far beyond the capability of current simulations to resolve, so different approximations are made for simulations to be computationally viable, usually either taking a zero thickness limit or a fixed comoving width; the outcome has been quite different predictions in the literature for both gravitational wave signatures and string axion (dark matter) radiation. A recent consortium paper (arxiv.org/1910.01718.pdf) addresses this controversy by using AMR to simulate global (or axion) string configurations and networks, allowing the numerical grid resolution to adapt to the scale of the problem. Figure 1 shows a simulation of an oscillating sinusoidal string, with massless radiation emanating out in a quadrupole pattern. Performing a careful quantitative analysis of this axion radiation led to the conclusion that the evolution of global strings tends towards that predicted in the thin-string limit, with additional radiation damping or backreaction. This provides a significant step towards resolving this issue and making precise predictions (in this case) of the small axion mass, which is a key parameter defining current dark matter axion searches.

This research is also being taken forward with much larger simulations of full cosmic string networks (see Figure 2), using diagnostic tools developed to measure competing massless (axion or GWs) and massive (Higgs) decay channels. This cosmic string work with the flexible GRChombo code has proved to be an important exemplar for testing new HPC capabilities, particularly for the COSMOS Intel Parallel Computing Centre collaboration on in-situ visualization (with the Intel Viz team, TACC and Kitware); this pioneering work using the OSPRAY ray-tracing visualization package has been demonstrated by consortium members at recent international supercomputer conferences and at SIGGRAPH 2019. We have also initiated a collaboration with the Chombo team at Lawrence Berkeley National Laboratory to speed up core libraries.

---

The evolution of high mass-ratio circumbinary discs

In recent years the detection of a large variety of structures in protoplanetary systems (once believed to be featureless), the detection of a high number of exoplanets with properties much different from those detected in the solar system, and the detection of gravitational waves from binary black holes have all raised a similar fundamental research question: “How does the mutual interaction between a binary and its surrounding gas shape the systems we observe?”

The conservation of angular momentum during the infall of material on to a binary (which in this context can mean two stars, two black holes, or a star and a planet) forces the material to form a disc structure that orbits the binary, and the dynamics of these so-called circumbinary discs have been studied for many years. The binary exerts a so-called “tidal torque” on the accretion disc, carving a gap or even a cavity in the disc. The disc in turn exerts a back-reaction torque on the binary, changing its orbit, causing it to migrate and/or changing its eccentricity. The traditional theoretical approach requires some simplifying assumptions, in order to make the equations analytically tractable, such as assuming that the mass of the secondary, or the mass of the disc, is much smaller than the primary. When these assumptions are not satisfied the theoretical predictions become unreliable, and we must turn to numerical hydrodynamic simulations in order to investigate the problem properly. However, the time-scales involved in this problem are long (at least hundreds of binary orbits), and high numerical resolution is required to resolve the tidal torques correctly. As a result, accurate hydrodynamical simulations of binary-disc interactions are extremely computationally challenging.

By exploiting the power of DiRAC, we were able to perform a suite of 3-D Smoothed Particle Hydrodynamics simulations which follow the evolution of a wide range of different circumbinary discs for long timescales (typically 2000 binary orbits). Our main aim was to understand the physical mechanism responsible for the formation of the so called “horseshoe” structures that have been seen in a large number of recent simulations of circumbinary discs, and which have been proposed as an explanation for the asymmetric structures observed by ALMA in numerous planet-forming systems. Our simulations (shown in Fig.1) showed that horseshoe structures are associated with the rapid growth of the cavity size when the disc becomes eccentric, if the binary companion is sufficiently massive. The long-term evolution of our simulations suggests that these features are actually transient, but they can survive for thousands of binary orbits. In real systems this corresponds to a significant fraction (~10%) of the ~Myr lifetimes of protoplanetary discs, which may explain why these structures are only seen in a sub-set of observations. Our simulations have provided important new insight into an important theoretical problem, as well as raising some interesting new questions. In particular, we unexpectedly discovered a very strong relationship between the growth of the disc cavity and its eccentricity, and our future work with DiRAC will investigate the physics behind this result in much greater detail.

---

Astrophysical Maser Flares from Shock Impact

One mechanism that has been put forward as a generator of maser flares in massive star-forming regions is the impact of a shock of modest speed on a cloud with a typical size of a few astronomical units. The sequence of diagrams on the left show the effect on images of such a cloud if struck by a 7.5km/s shock that moves up the z-axis (bottom to top). There is a frame change at the 50% shocked point from a view where the shock appears to move, at lower precentages, to a view where the remains of the uncompressed cloud flows into a stationary shock at the cloud mid-point (for shocked proportions > 50%). The shock is modelled as isothermal, with a compression factor proportional to the square of the Mach number. In this case, the Mach number is 3, giving a density enhancement of 9 in the shocked part of the cloud, which increases from 10% of the cloud in the top diagram to 80% in the bottom panel. Compression is applied to the node distribution of the model prior to triangulation. The colour palette is the base-10 logarithm of the maser specific intensity relative to the saturation intensity so the brightest parts of the image increase from about 0.001 of the saturation intensity in the undisturbed cloud to of order 104 times this parameter in the 80% panel, so that the shock increases the brightest pixel intensities by of order 10 million times. The effect on the flux density, averaged over the whole cloud is lower, but this impact still represents a very convincing mechanism for generating high-contrast maser flares. If one model unit on the axes of each panel is equal to half an astronomical unit, then the rise time of the flare is approximately the time taken for the shock to cover 1 astronomical unit, or 231 days. The effectiveness of the flare mechanism relies on a switch, due to the shock, from a weakly saturated initial cloud (approximately spherical) to a quite strongly saturated, and highly oblate, final cloud.
Data from DiRAC is available for (original) optical depths ranging from 0.1-30.0, and a number of different shock speeds.

---

Dust destruction in supernova remnants

It is well established that (sub-)micrometre sized dust grains can form in over-dense gas clumps in the expanding ejecta of supernovae remnants. However, highly energetic shock waves occur in the ejecta which can potentially destroy a large fraction of the newly formed dust grains. The gas in the ejecta is heated up to billions of Kelvin and is accelerated to a few hundreds of kilometre per seconds which causes thermal and kinematic sputtering of the dust grains. Moreover, the dust grains can collide among each other with high velocities and get fragmented or even vaporized. Previous predictions for the dust survival rate depend strongly on initial parameters and range from less than 0.1 % to 99 %. The net dust survival rate is crucial for determining whether or not supernovae significantly contribute to the dust budget in the interstellar medium.

In order to model a shock wave interacting with an ejecta clump we performed hydrodynamics simulations using the grid-based code AstroBEAR (Fig. 1). Afterwards, dust motions and dust destruction rates are computed using our newly developed external, post-processing code Paperboats, which includes gas and plasma drag, grain charging, kinematic and thermal sputtering as well as grain-grain collisions. We used DiRAC HPC Facilities to determine the dust survival rates for the oxygen-rich supernova remnant Cassiopeia A for a huge range of parameters, including initial grain sizes, dust materials and clump gas densities.

We find that up to 40 % of the silicate (Fig. 2) and up to 30 % of the carbon dust mass is able to survive the passage of the reverse shock. The survival rates depend strongly on the initial grain size distribution, with ∼10−50 nm and ∼0.5−1.5 μm as the grain radii that show the highest surviving dust masses. The dust processing causes a rearranging of the initial grain size distribution. Our results show that grain-grain collisions and sputtering are synergistic and that grain-grain collisions can play a vital role in determining the surviving dust budget in supernova remnants.

References:

• Kirchschlager, F., Schmidt, F. D., Barlow, M. J., Fogerty, E. L., Bevan, A., & Priestley, F. D. 2019, MNRAS, 489, p.4465-4496

---

Gaseous Dynamical Friction in the Early Universe

Dynamical friction is the process by which a massive perturber, moving through some background medium, gravitationally interacts with that medium, producing a net retarding force to its motion. When the background medium is gaseous, the pressure forces present in the gas must be included in modelling the response of the medium, and so impact the resultant force. To numerically model this process, we used the gravo-hydrodynamic code GIZMO, which uses state of the art Lagrangian hydro solvers, to test highly idealised wind tunnel-like setups against the analytic predictions for this net force. We found that the numerical results significantly under produce the force predicted by analytic treatments of the problem. These solvers did not reproduce the predicted structure of the wake, instead smoothing out predicted density profiles, and producing unexpected bow wave structures. These results were found in regimes where the analytic prediction should hold. As mentioned above, dynamical friction is critical in driving the mergers of dark matter sub-halos into their larger hosts, as it provides a crucial mechanism by which angular momentum can be transferred away from the orbiting substructure, allowing it to spiral into the central object. This is particularly important in the early Universe, when gas fractions were higher. The reduced retarding force (shown below in Figure 1 for different perturber velocity Mach numbers) suggests that modern cosmological simulations, that are run using these solvers, could be underestimating the merger rates of early dark matter structure, and so incorrectly recovering the evolution of early galaxies.

---

The b1 resonance in coupled-channel scattering from Lattice QCD

The vast majority of hadrons are resonances, seen experimentally through decays into a multitude of other hadrons in various angular distributions and corresponding to pole singularities in scattering amplitudes. Resonances which decay into two pseudoscalar mesons (mesons with J^P=0^–, where J is spin and P is parity) have been studied extensively using first-principles lattice QCD calculations over the last decade. However, all but the simplest resonances have decay channels which involve hadrons with non-zero spin and, in particular, those with ‘unnatural’ spin parity, J^P=0^–, 1⁺, 2^–, are forbidden to decay to two psuedoscalars. Calculations of these scattering channels are significantly more complicated.

In a recent publication [arXiv:1904.04136, Phys. Rev. D 100, 054506 (2019)], we studied the axial-vector resonance (J^P=1⁺). Experimentally its decay channels include πω and πφ, where π is a pseudoscalar meson (J^P=0^–) and ω, φ are vector mesons (J^P=1^–). The non-zero spin of the vector mesons means that the decay may feature a superposition of orbital angular momentum L=0 (S-wave) and L=2 (D-wave). Working with heavier-than-physical light quarks such that the π mass is 391 MeV, we performed the first lattice QCD computation of the relevant coupled scattering amplitudes – the resulting squared-amplitudes for πω in S-wave and D-wave, and πφ in D-wave are shown in the figure below (upper panels) labelled by ^2S+1L_J. Notably, we observe a bump-like enhancement in S-wave πω, suggestive of a narrow resonance, and we indeed find a pole in the analytic continuation of the amplitudes to complex energies (lower panel) – the real and imaginary parts are related to the mass (m_R) and width (Γ_R) respectively, and the strengths of the coupling to each channel (|c|) are also shown. Interestingly, the absence of significant enhancement in S-wave πφ suggests no evidence for a Z_S, proposed as an analogue of the enigmatic Z_c(3900) observed in πJ/ψ. This work opens the door to investigating the π₁(1600) exotic light resonance, nucleon resonances such as the Roper, and puzzling structures observed in the charmonium and bottomonium sectors.

---

Electron acceleration and transport in turbulence

Highly energetic electrons are present in many different astrophysical systems, from solar flares and supernovae to planetary magnetospheres and the intra-cluster medium. Plasma turbulence is ubiquitous in these systems, and it is fundamental for the production of energetic particles and their transport properties.

The problem of particle acceleration in astrophysical turbulence was first considered by Fermi, who considered particles interacting with randomly moving scattering centres, thereby gaining energy stochastically [1]. If the scattering centres have some coherent, large-scale motion, the energisation is more rapid [2]. Since Fermi’s original ideas, the topic has been extensively investigated.

Using DiRAC we have studied the acceleration and transport properties of transrelativistic electrons in plasma turbulence [3] using hybrid PIC and test particle simulations. We have discovered that turbulence strength is a key parameter that controls electron acceleration, with rapid acceleration due to particle trapping happening when the turbulence level is high. Figure 1a shows two examples of open and trapped electron trajectories. Figure 1b shows a detail of electron trajectory with the trapped part highlighted in red. It can be seen, in the figure inset, that, when the electron is trapped, fast quasi-monotonic energisation is achieved. The fast energisation stage resembles a first-order Fermi process as observed by other authors in simulations of magnetic reconnection (see Ref [3] for further details).

The importance of turbulence strength to activate fast electron acceleration has important consequences for heliospheric and astrophysical plasmas, with implications for understanding, for example, the physics of solar flares and radio emission from galaxy clusters.

References:

• [1] E. Fermi, On the origin of cosmic radiation, Physical Review 75:1169-1174 (1949)
• [2] E. Fermi, Galactic magnetic fields and the origin of cosmic radiation, The Astrophysical Journal 119:1 (1954)
• [3] D. Trotta, L. Franci, D. Burgess and P. Hellinger, Fast acceleration of transrelativistic electrons in astrophysical turbulence, The Astrophysical Journal, 894:2-136 (2020)

---

Extreme QCD: Quantifying the QCD Phase Diagram

The FASTSUM collaboration uses 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 inside these hadrons is Quantum ChromoDynamics, “QCD”. We are particularly interested in the behavior of QCD as the temperature increases to billions, and even trillions of Kelvin. 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 (near New York).

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 pull quarks apart, whereas at temperatures above the “confining” temperature T_c, 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. Understanding how these masses change with temperature can give deep insight into the thermal nature of QCD and its symmetry structure.

Our most recent results are summarized in the plot where we show the proton mass as a function of temperature up to around T_c in the top left panel. The other panels show results for hadrons containing various numbers of strange quarks. As can be seen, the mass of positive parity states have fairly constant masses, whereas their negative parity partners’ masses decrease substantially as the temperature increases until they become nearly degenerate with their positive parity partners at around T_c, as dictated by the global chiral symmetry of the QCD Lagrangian. Beyond T_c masses are difficult to determine which is consistent with the QCD force becoming so weak that the hadrons break apart. This calculation supplements our theoretical understanding of QCD with quantitative detail..

---

Deep inelastic structure of the proton

Hadron structure functions are ubiquitous in the description of leptonic interactions with hadrons, encoding: elastic form factors, inclusive electro- (and photo-) production of resonances, diffractive processes and Regge phenomena, and partonic structure in deep inelastic scattering. Lattice QCD, however, has really only been able to probe some limited kinematic corners of the all-encompassing structure functions primarily being limited to elastic form factors and low(est) moments of leading-twist parton distributions. The particular interest in partonic structure has motivated a number of strategies to overcome limitations in the lattice formulation, including: the Euclidean hadron tensor, lattice OPE, heavy-quark Compton amplitude, symmetry-improved operator construction, factorisable matrix elements, and quasi-PDFs and related quantities. We have embarked on a complementary program to extract the forward Compton amplitude in the unphysical region [1]. Computationally, we are able to take advantage of the efficiency of the Feynman-Hellmann approach to hadron structure and avoid the need to compute 4-point functions. Here we highlight [2,3] some recent progress towards revealing scaling and higher twist-phenomena in the low-order moments of the Compton amplitude from lattice QCD.

References:

• [1] A. Chambers et al., Phys. Rev. Lett. 118 (2017) 242001.
• [2] A. Hannaford-Gunn et al., PoS (LATTICE2019) 278.
• [3] R. Horsley et al., ibid 137.

---

New LIGO-Virgo gravitational waves

The third LIGO-Virgo observing run (O3) began in April 2019, and has since announced three gravitational-wave observations – the binary-neutron star merger GW190425, the binary-black-hole merger GW190412, and most recently GW190812, the merger of a 23-solar-masss black hole with a 2.6-solar-mass object, that could be either the lowest-mass black hole yet observed through gravitational-waves, or the highest-mass neutron star ever found. Measuring the properties of these sources, in particular the masses of the objects, relied on several theoretical signal models, including the family of phenomenological models that is tuned to fully general relativistic numerical simulations of black-hole binaries that are performed in this project.

Over the last year, the models have been upgraded to include subdominant harmonics (which were essential to the analysis of GW190412), and to treat black-hole—neutron-star systems (of which GW190812 may be one). O3 data continue to be analysed using these models, and a full catalogue of all confirmed observations will be published when it is complete.

---

UNITY: Multi-wavelength simulations of relativistic cosmology

Gravitation in our Universe is best-described by Einstein’s theory of General Relativity. Cosmological simulations typically use a Newtonian approximation for the gravitational field however, which is far simpler to solve, and which proves to be an excellent approximation over a wide range of distance scales. Forthcoming surveys of the large-scale structure of the Universe (e.g. galaxy surveys such as Euclid and LSST, and intensity mapping surveys with MeerKAT and SKA) are now becoming large enough to capture the very largest distance scales in the observable Universe however, which inherently require a relativistic treatment.

To help us precisely understand the effects of the relativistic corrections on observations that will be made by these surveys, we have produced the largest ensemble of fully-relativistic N-body simulations to date, using the cutting-edge gevolution code. This includes a rigorous treatment of relativistic degrees of freedom, such as massive neutrinos, as well as solving for photon paths using the full set of relativistic fields that exist on large scales in the Universe. The result is a particularly accurate representation of a range of different cosmological large-scale structure observables that fully accounts for subtle GR effects.

As part of the UNITY project, we ran an ensemble of 51 gevolution simulations on the DiRAC Data Intensive system in Leicester, with corresponding lightcones and tomographic maps of various observables taken from each and stored for further processing to add survey-specific observational features, such as different populations of galaxies. The simulations were divided into two sets: an ensemble containing different random realisations of the initial conditions, which can be used to study the statistical properties of observables and methods used to recover them, and a set with fixed initial conditions but systematically different values of cosmological parameters, allowing us to perform a precision study of the cosmology-dependence of the various relativistic correction terms.

Two key works are based on these simulations are currently being prepared for publication: a technical description of the simulations and study of the relativistic dipole in the galaxy distribution (L. Coates et al, in prep.), and a detailed study of GR-induced multipoles in the power spectra and two-point functions of multiple galaxy populations (C. Guandalin et al, in prep.).

---

Impact of Reionization on the Intergalactic Medium

The latest measurements of CMB electron scattering optical depth reported by the Planck satellite significantly reduces the allowed range of HI reionization models, pointing towards a later ending and/or less extended phase transition than previously assumed. During HI reionization the intergalactic medium (IGM) is photoheated to ~10⁴ K, and owing to long cooling and dynamical times of the diffuse gas, comparable to the Hubble time, its vestiges persist in the IGM to much later times. Existing and upcoming observations of the Lyman-α (Ly-α) forest at high redshift can be used to detect these extant thermal signatures of HI reionization by comparing them with accurate hydrodynamical models of HI reionization. This can provide new independent constraints not only on the timing of HI reionization but also on the sources that reionized the universe. However to achieve these goals reliable and precise models of the Ly-α forest are needed to be compared with the observations. The challenge for these models is to reach enough dynamical range to attain the necessary precision in the Ly-α forest predictions and at the same time to obtain the number of models required to guarantee a good sampling of the parameter space.

We have run a grid of state-of-the-art cosmological hydrodynamical simulations in order to explore the physical parameter space characterized by the onset of reionization, its duration and the heat injected into the IGM. The new methodology used allowed us to take into account the complex inhomogeneities of the reionization process and to explore the full physical parameter space very efficiently. The creation of an accurate emulator of the Ly-α statistics with these simulations will be crucial to provide accurate fits to Ly-α observations.

---

Tidal disruption by supermassive black holes

Tidal disruption events (TDEs) occur when a star passes too close to a supermassive black hole (SMBH). By “too close” we mean that the tidal field resulting from the black hole’s gravity is so strong that it can overcome the star’s self gravity and rip the star apart. Stars that pass close to the black hole typically do so on near-parabolic orbits, and thus approximately half of the debris from the disrupted star is on bound orbits, and the other half is on unbound orbits with respect to the black hole. The bound portion of the stream is thought to return to the black hole and circularise into an accretion disc which powers a luminous flare in an otherwise dormant galactic nucleus. Thus TDEs can be used to probe the SMBH mass function, the properties of individual stars, and stellar dynamics in galactic nuclei. Upcoming missions will detect thousands of TDEs, and accurate theoretical modelling is required to interpret the data with precision. In Golightly et al. (2019 ApJL 882, 2, L26) we present numerical simulations performed on Data Intensive at Leicester (DI@L) of TDEs with accurate structures for the simulated stars that are taken from the MESA stellar evolution code. This is in contrast to previous numerical simulations which employed simple polytropic models for the stars. We find that the more accurate stellar density profiles yield qualitatively-different fallback rates when compared to polytropic density profiles, with only the fallback rates from low-mass (less than the mass of our Sun), zero-age main sequence stars are well fit by a standard polytopic model. Stellar age has a strong affect on the shape of the fallback curve, and can produce characteristic timescales (e.g., the time to the peak of the fallback rate) that greatly differ from the polytropic values. We further investigate whether polytropic models for the fallback rate, which are routinely used to infer properties of observed systems from observational data, can lead to inaccurate measurements of e.g. the black hole mass. In the Figure we show the results of two simulations with different black hole masses (differing by a factor of 5) in one of which the star was modelled as a polytrope and in the other the star was modelled with an accurate structure from the MESA stellar evolution code. The similarity of the fallback curves for such a different black hole mass suggests that detailed modelling that captures realistic stellar density profiles is required to accurately interpret observed lightcurves of TDEs.

---

DiRAC enables prediction for matter-anti-matter asymmetry in the Standard Model

The Standard Model (SM) of Elementary Particles has been extremely successful in correctly predicting and describing the properties of elementary particles studied at experimental facilities all around the world. The SM itself covers particles like quarks and electrons, for instance, and their anti-particles. While both particles and anti-particles behave largely in the same way, the SM allows for tiny differences in their behaviour.

This is however extremely important since it is an essential ingredient for being able to describe the abundance of matter over anti-matter observed in our universe. One of the big questions in modern particle physics is therefore, whether the SM correctly describes the observed symmetry breaking, or whether new, not yet known and undiscovered physics is driving it.

DiRAC resources have now enabled a detailed study of the symmetry breaking by considering the decay of kaons into two pions, as studied experimentally at CERN and Fermilab two decades ago. Both particles are mesons, i.e. they consist of a quark and an anti-quark, respectively. Researchers from Brookhaven National Laboratory, CERN, Columbia University, the University of Connecticut, the University of Edinburgh, the Massachusetts Institute of Technology, the University of Regensburg and the University of Southampton (RBC/UKQCD Collaboration) set out to compute the amount of symmetry violation present in the SM by means of large-scale computer simulations of Lattice Quantum Chromodynamics (QCD). QCD is the part of the SM that describes the strong interactions of quarks (and anti-quarks) via the exchange of force carriers called gluons. The group of researchers 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.

To date RBC/UKQCD’s prediction is the only one of its kind and therefore highly anticipated. At the level of precision achieved in their calculation the result agrees with the experimental measurement. On the one hand this result is yet another piece of evidence for the unparalleled success of the SM as the theory of elementary particle physics. On the other hand, we expect the SM to be incomplete in describing all aspects of particle physics and CP violation remains an excellent place to look for limitations of the SM. The Collaboration is therefore now working on ways to improve the precision of their result.

Bibliography:

• “Direct CP violation and the ΔI=1/2 rule in K→ππ decay from the Standard Model”, arXiv: 2004:09440
• “Standard Model Prediction for Direct CP Violation in K → ππ Decay” Phys.Rev.Lett. 115 (2015) 21, 212001, arXiv:1505.07863

---

Ion-neutral decoupling in the nonlinear Kelvin—Helmholtz instability: Case of field-aligned flow

In many areas of the solar atmosphere, including prominences, and many other astrophysical systems the majority of the fluid is neutral with only a small fraction being plasma. As with fully ionized plasmas, plasmas that are only partially ionized can be subject to flow instabilities which result in turbulence developing in the plasma. However, the coupling between the neutral component of the fluid and the magnetic field is not perfect, with the coupling coming as a result of collisions between neutral and charged species. The consequence of this is the fluids are strongly coupled to each other (resulting in the neutrals strongly feeling the Lorentz force) for slow dynamics, but if dynamic motions are quick enough the neutral fluid will not feel the magnetic field and behave purely hydrodynamically.

In this study we performed an extremely high-resolution 2D calculation, which allowed the dynamics at both coupled and decoupled scales to be simultaneously resolved. The basic set up had a velocity and density jump and a magnetic field aligned with the flow. The fluid was only 1% ionized.

The most surprising results of this study came from the thermodynamics. We found that heat transfer in the system could occur strongly across the magnetic field because the neutral fluid

How turbulent heating occurred also threw up a surprise. In standard turbulence, it is normally the smallest scales which drive the dissipation, but we found the frictional heating in partially ionized plasma is determined by the larger vortex scale instead of the small scale.

Based on results published in “Ion-neutral decoupling in the nonlinear Kelvin—Helmholtz instability: Case of field-aligned flow” by Andrew Hillier, Physics of Plasmas (2019). (The article can be accessed at https://doi.org/10.1063/1.5103248.)

---

Simulation of Sp(2N) gauge theories for Composite Higgs models

Understanding the nature of the standard model Higgs boson is still an open problem. An appealing possibility is that the Higgs boson be a composite particle resulting from a novel strong interaction. The lack of observation of otherwise unexplained particles interacting through this conjectured novel force and with mass comparable to that of the Higgs boson would then require a mechanism that keeps the Higgs parametrically light. Such a mechanism can be provided by the spontaneous breaking of the global flavour symmetry in the new interaction. In this framework, the Higgs boson is interpreted as a Pseudo-Nambu-Goldstone Boson (PNGB) of the novel interaction. Among candidate realisations of PNGB compositeness are Sp(2N) gauge theories. The two simplest theories in this class are Sp(2), coinciding with SU(2), and Sp(4). Given the strong nature of the novel force, in order for the former to be validated experimentally as viable theories, first-principle calculations would need to be performed that allows us to determine quantitative predictions from these models. These calculations can be executed numerically formulating the models on a spacetime lattice.

Our project has provided the first numerical calculation in the Sp(4) gauge model with two fundamental flavours and the first determination of its lowest-lying meson spectrum and of the corresponding decay constants in the continuum limit. This has enabled us to determine the coefficients of an effective field theory that can be used to compare the predictions of the model with experimental results [1]. Our study of the Sp(2N) Yang-Mills model combined with results in the literature for SU(N) and SO(N) gauge theories has unravelled that in both three and four spacetime dimensions the ratio of the mass of the tensor glueball over the mass of the scalar glueball is independent of the gauge group [2]. This result provides further insights on the phenomenon of confinement, which characterises non-Abelian gauge models and provides striking experimental signatures such as jets in high-energy experiments involving the standard model strong force.

Our project also provided the first calculation in a composite Higgs scenario of the scattering amplitude of two PNGBs in the SU(2) model. The amplitude controls the production cross section of a new resonance expected to be produced at the LHC in a Composite Higgs scenario. The results therefore constrain the viability of the model and shed light on a key observable to test composite Higgs models using lattice simulations.

References:

• [1] E. Bennett et al., Sp(4) gauge theories on the lattice: Nf=2 dynamical fundamental fermions, JHEP 12 (2019) 053 [arXiv:1909.12662]
• [2] E. Bennett et al., On the colour dependence of tensor and scalar glueball masses in Yang-Mills theories, Phys. Rev. D Rapid Communications to appear [arXiv:2004.11064]
• [3] V. Drach et al., Resonance study of SU(2) model with 2 fundamental flavours of fermions, Proceedings of Science [arxiv:1910.13847] – paper in preparation.

Some of the Research Highlights from our users in 2018:

• HPQCD: The weighty issue of quarks.
• QCDSF-UKQCD: Light quark flavour breaking effects in B meson decay constants.
• Gravitational Waves 2018.
• Exomol: Computational spectroscopy of Exoplanets.
• UK MHD Consortium: 1 Solar and Planetary Interiors.
• Leicester Highlight: Planet-driven warps in protoplanetary discs.
• COSMOS Consortium.
• Modelling Planetary Collisions with the SWIFT code.
• The FABLE simulations of galaxies, groups and clusters.
• Fast and energetic AGN-driven outflows in simulated dwarf galaxies.
• Studying the origin and amplification of magnetic fields in galaxies.
• Turbulence, Shocks and Dissipation in Space Plasmas
• Longest ever simulation of a planet in a protoplanetary disc
• Galactic Archaeology in the Gaia Era

2018 list of publications

---

The weighty issue of quarks

Electrons make circular tracks when they move perpendicular to a magnetic field in a particle detector. From the radius of the circle, knowing the electric charge, we can determine the electron mass. Quarks, on the other hand, are never seen as free particles. Instead, all we see in our particle detectors are the bound states of quarks, called hadrons. We can measure hadron masses, but how do we then determine, or even define, the mass of a quark?

The answer is that we must define the mass as the parameter that enters the equations of QCD, the theory of the strong interaction that binds quarks into hadrons. We can solve the equations of QCD using the numerical approach of lattice QCD and hence determine the quark mass parameter that we need to put into the theory to give the correct mass for a calibration hadron containing that quark. Since experimentalists can measure some hadron masses very well, we can do this accurately. For example, to fix the charm quark mass, m_c, we use the mass of the η_c (pictured in Figure 1).

To be of use to other theorists, however, we must convert the lattice quark mass into that appropriate to QCD defined in continuous space-time. We must do this accurately too, to avoid introducing a sizeable uncertainty into the result.

HPQCD developed one of the main techniques used to determine quark masses (the JJ method) and most recently used this, on DiRAC, to determine m_c to 1% in 2014 (arXiv:1408.4169). This year we used a completely different technique (called RI-SMOM) to ‘match’ our lattice quark mass to that of continuum QCD, directly on the lattice (arXiv:1805.06225). We were careful to quantify and remove lattice artefacts, which had been ignored in previous calculations. Our new 1% accurate result agrees very well with our 2014 value and also with a third, MRS, method developed by FNAL/MILC/TUM in 2018. This gives a world-average for m_c of 1.2753(65) GeV (1.36 times the mass of the proton). The quark masses need to be accurately known so that we can test whether the newly-discovered Higgs boson is behaving as we expect when it decays to hadrons via quark-antiquark pairs of different flavours.

---

Light quark flavour breaking effects in B meson decay constants

In the standard model of particle physics, there are six types of quarks whimsically called up u, down d, strange s, charm c, bottom b and top t in order of increasing mass. These bind together to form two classes of hadrons, mesons consisting of a quark – antiquark pair and baryons with three quarks. An example of a baryon is the proton (uud) while examples of mesons are the pion (ud), kaon (us), B (ub) and Bs (sb). These mesons decay to a lepton (usually the electron) and a neutrino. The decays are parametrised by ‘meson decay constants’, from which, together with experimentally determined leptonic decay rates elements of the CKM matrix can be found. These are fundamental parameters of the Standard Model and their determination can help to isolate possible new physics. Presently recent B physics experiments have left us with a number of heavy flavour physics puzzles which will be more intensively investigated by Belle II.

In our Lattice QCD programme we have developed a method by expanding about a point in the light quark mass parameter space where the u, d, s quark masses are initially set equal, and then extrapolating to the physical quark masses while keeping the average quark mass constant. The original expansions were for hadron masses, and lead to a ‘fan’ plot where the different hadron masses radiate from their common point. They were extended to decay constants for the pion and kaon mesons, [1], as shown in the left panel of the figure. We have now further extended this programme to the heavier quark masses and in particular the b, [2]. In the right panel of the figure we give the equivalent plot for the B and Bs mesons, together with the present (FLAG16) values.

[1] V.G. Bornyakov et al., Phys. Lett. B 767 (2017) 366.

[2] S. Hollitt et al., PoS(LATTICE2018)268, arXiv:1811.06677 [hep-lat].

---

Gravitational Waves 2018

The Advanced LIGO and Virgo detectors were undergoing an upgrade during 2018, so the year was devoted to completing analyses of data from the first two observing runs (O1 and O2), and preparations for the third observing run (O3), due to start in April 2019.

The chief results from the LIGO and Virgo collaborations were the final estimates of the physical properties of the binary-neutron-star merger, GW170817, including stronger constraints on the neutron star tidal deformability, and therefore its equation of state and radius; and the release of the first catalogue of observations from O1 and O2, based on a re-analysis of all data, which in turn had undergone improved characterisation and cleansing. This re-analysis uncovered new signals in the data, and doubled the number of observed black hole binary mergers from five to ten. This included re-classifying the marginal “LIGO-Virgo transient” from October 12, 2015 as a bona-fide GW signal, and among the new signals, GW170729, which is likely the most massive binary yet observed, at around 85 solar masses. All of this work used the waveform models that were calibrated to numerical-relativity simulations performed on DiRAC.

Modelling work also continued. The first inspiral-merger-ringdown model to include higher signal harmonics (PhenomHM) was published in Physical Review Letters, and is now being used to update the measurements of the properties of GW170729. The bulk of the most extensive and systematic set of simulations of precessing binaries was completed on Cosma6, and a preliminary model of the precession effects through merger and ringdown (the first of its kind) was completed, with publication expected in 2019, for use in analysis of O3 observations.

---

Exomol: Computational spectroscopy of Exoplanets

Molecular spectra play a critical role for the retrieval of atmospheres of exoplanets and cool starts. However, the lack of laboratory data on many of molecular species severely limits models for objects as diverse as Jupiter, exoplanets and brown dwarfs. The UCL ExoMol team is world leader in providing molecular line lists for exoplanet and other hot atmospheres. The ExoMol procedure uses a mixture of ab initio calculations and available laboratory data. These line lists form the input for opacity models for cool stars and brown dwarfs as well as for radiative transport models involving exoplanets. So far ExoMol has collected molecular opacities for more than 50 molecules (130 isotopologues): 30 line lists have been generated using the ExoMol methodology and other line lists were taken from external sources or from our work predating the ExoMol project.

Ethylene and Acetylene are important absorber in hot Jupiter exoplanets, brown dwarfs and cool carbon stars. As part of our Hydrocarbons DiRAC project, we have new line lists for these molecules, each of which contain over 30 billion transitions. In order to accomplish this data intensive task, some of the UK’s most advanced supercomputers have been used, provided by the DiRAC project. Even larger line lists have been completed. For example, the methyl chloride’s line list contains over 300 billion lines.

Many of those exoplanets discovered thus far are categorized as rocky objects with an atmosphere. Most of these objects are however hot due to their short orbital period. Models suggest that water is the dominant species in their atmospheres. The hot temperatures are expected to turn these atmospheres into a (high pressure) steam bath containing remains of melted rock. The spectroscopy of these hot rocky objects will be very different from that of cooler objects or hot gas giants. Molecules suggested to be important for the spectroscopy of these objects are reviewed together with the current status of the corresponding spectroscopic data. We have started building a comprehensive database of linelist/cross sections applicable for atmospheric models of rocky super-Earths as part of the ExoMol project.

An efficient usage of line lists comprising billions of lines is extremely challenging. To address this challenge of Big Data, we have developed a new Fortran 2003 program ExoCross. ExoCross takes line lists as input and returns pressure- and temperature-dependent cross sections as well a variety of other derived molecular properties which depend on the underlying spectroscopic data. These include state-dependent lifetimes, temperature-dependent cooling functions, and thermodynamic properties such as partition functions and specific heats and is designed to work with the recently proposed method of super-lines. It is capable of simulating non-LTE spectra using a simple two-temperature approach.

Supported by the DiRAC high performance facilities, we have organised and run a workshop “Digital Exoplanets” in Prague, January 2019 (www.digitalexoplanets.org). The workshop gathered over 50 computational exoplanetary scientists (with more than half career researchers) and focused on the data for atmospheric models of exoplanets and cool stars. The main theme of the workshop was the open access to the data and computer codes. This was the first workshops of this kind, focusing not only on observations and theory, but mostly on ‘how to actually run these codes’. The participants had the opportunity to learn how to use (install, run, debug etc) different main stream atmospheric and spectroscopic codes during the hands-on sessions, which were the central part of the workshop. These practical sessions were run on the resources provided by DiRAC team, with the help of the DiRAC team, which is greatly appreciated.

---

UK MHD Consortium: 1 Solar and Planetary Interiors

One of the great challenges in the study of planetary interiors is to understand the generation of magnetic fields in the interiors of the giant gas planets. Wieland Dietrich and Chris Jones from the University of Leeds have been addressing this problem through high resolution magnetohydrodynamic models, with the novel feature of incorporating realistic radially-dependent profiles for the electrical conductivity (Icarus 305 (2018) 15-32). This allows the study of a range of models in which the conductivity ranges from Jupiter-like, with only a thin outer non-conducting shell to Saturn-like, with a thick non-conducting shell. They find Jupiter-like steady dipolar fields, quadrupolar dynamos for profiles between those of Jupiter and Saturn, and dipolar fields again for Saturn-like profiles. The non-axisymmetric components of the Saturn model are, though, more pronounced than in Saturn itself, prompting the suggestion that a stably stratified region – postulated to exist in Saturn, but not yet incorporated in the model – may be of significance.

---

Leicester Highlight: Planet-driven warps in protoplanetary discs

Since their discovery in the 1990s our knowledge of extra-solar planets has increased at a staggering rate, and we now know that most, perhaps all, stars host planetary systems. These planets form in cold discs of dust and gas around young stars, and in recent years our observations of these protoplanetary discs have also advanced dramatically. We typically treat planet-forming discs as flat, smooth structures but recent observations have revealed a wealth of sub-structures in these discs, including spiral arms, gaps, misalignments and warps. While we can explain the origin of some of these features, many are still open questions. One of the best-studied protoplanetary discs is around the nearby star TW Hya, and this almost face-on disc shows evidence of gaps (possibly due to planets) as well as a localised dark region in the outer disc. Observations with the Hubble Space Telescope found that this dark region is moving around the outer disc, but it is seen to move much more quickly than gas orbits at that radius. This strongly suggests that the dark region is actually a shadow cast by something in the inner disc (where the gas orbits faster). This was interpreted as evidence of a disc warp, such that the inner part of the disc is inclined with respect to the outer disc. The relative misalignment between the inner and outer disc then results in a shadow being cast on the outer disc. Such a configuration could be caused if a misaligned planet is present in this system, where the planets orbit is tilted away from the mid-plane of the disc. The motion of the inner, tilted disc and hence the shadow would thus be governed by the gravity of this misaligned planet, and it might be possible to infer properties of the planet from the observed shadow.

Simulating such a configuration is technically challenging, because we need to resolve the out-of-plane interactions between the disc and the planet, as well as including the influence of the outer disc. In practice this requires high resolution 3-D hydrodynamic simulations covering a large dynamic range. Using DiRAC, we were able to run a suite of simulations to study the interaction between misaligned planets and their parent discs in detail, investigating the importance of different planet masses and inclinations on the subsequent disc structure. For planets that were massive enough to carve the disc into an inner disc and an outer disc (as seen in Fig.1), we confirmed it was possible to drive the configuration predicted from the observations of TW Hya. Importantly, the large spatial scales we were able to simulate allowed us to show that this was possible even for small planet inclinations. We subsequently ran detailed radiative transfer simulations to calculate how our simulated discs would appear to our telescopes, and confirmed that such a configuration is broadly consistent with what is observed in TW Hya. New observations are finding that misalignments of this type are common in planet-forming discs, and in the coming years DiRAC will allow us to understand the dynamics and evolution of these nascent planetary systems in detail.

---

COSMOS Consortium

The COSMOS consortium (dp002) continues to exploit DiRAC HPC Facilities and has made progress towards ambitious sciences in extreme gravity, inflation and the early universe, the study of the cosmic microwave sky and large-scale structure (including the Final Data Release for the ESA Planck satellite). Here, we highlight recent breakthroughs in general relativistic simulations of black holes and inflation using GRChombo, a flexible AMR numerical relativity code publicly released by COSMOS consortium members in early 2018 (see grchombo.org and github.com/GRChombo).

The groundbreaking discovery of a gravitational signal from a binary black hole merger in 2016 by LIGO is now followed by an exciting series of detections, including that of the first neutron star binary merger in 2018. With the current LIGO/VIRGO O3 science run with upgraded detectors detecting multiple events weekly, the new era of gravitational wave astronomy is well underway.

The GRChombo AMR code, originally developed and now anchored by COSMOS consortium members, is well-placed to exploit these new observations of gravitational waves (GWs). We have continued our long-standing industrial collaboration through the COSMOS Intel Parallel Computing Centre (IPCC) to modernize our research codes in cosmology and general relativity, focusing on the optimization of GRChombo for multi/many-core systems; for example, GRChombo is one of the few codes that can take full advantage of the many-core Xeon Phi systems. In addition, GRChombo continues to grow in functional capabilities, with the addition of a powerful MultiGrid initial conditions solver and the pioneering many-core OSPRay realtime ray tracer and volume renderer with in-situ capabilities (supported by the Intel Viz team).

The advent of GW astronomy has opened up a window into the one of the key questions in modern cosmology about the nature of dark matter. COSMOS investigators are taking a leading role in this search. One of the prime candidates for dark matter is an ultralight field which are sometimes called “axion-like particles” or ALPs. In recent papers (notably arXiv:1609.04724 and arXiv:1806.09367), we have investigated the formation of highly compact “axion dark stars” from these light fields, and the potential GW signals created by their collisions with neutron stars and black holes (arXiv:1808.04668, shown in the figure below), and with other axion stars (arXiv:1802.06733). In (arXiv:1904.12783), we showed that these fields can also interact directly with black holes to grow “scalar hair”, leading to a possibly detectable dephasing of the GW signals from BH binary mergers.

Beyond GW astronomy, COSMOS consortium investigators are also leading in the numerical study of fundamental questions in physics. For example, while inflation is now the paradigmatic theory of the early universe, providing a dynamical mechanism to explain why the universe is homogenous over cosmological scales and provide a source of the seeds of large scale structure. Nevertheless, can inflation begin in the first place, if the initial conditions are not homogenous? In recent work (e.g. arXiv:1712.07352), we classified the models of inflation which are not robust to inhomogenous initial conditions. Intriguingly, we found that low scale inflation – preferred by many string theory phenomenological models ­– is not robust, failing to inflate in the presence of large scale but subdominant inhomogeneities. 

In addition, we continue to take the lead in investigating the nature of gravity in higher dimensions. Following our work in showing the failure of weak cosmic censorship in 5D due to black ring instabilities, we continue to study the end states of higher dimensional rotating (Myers-Perry) black holes. We find that, similarly to the black ring systems, these black holes exhibit the Gregory-Laflamme instability, leading to a fractal instability with the ultimate end state likely to exhibit naked singularities in an asymptotically flat space time (shown above). Given these results in extra dimensions, it is clear that, if cosmic censorship is indeed true, then it must be a property specific to the three-dimensional space that we live in; without this essential property, Einstein’s theory of general relativity could lose its predictive power.

---

Modelling Planetary Collisions with the SWIFT code

We have produced the highest resolution simulations of the giant impact that is thought to have knocked Uranus’ spin axis into the plane of the solar system. The two orders of magnitude increase in resolution was enabled by our use of the SWIFT hydrodynamics and gravity code on the DiRAC facility in Durham. We incorporated equations of state for various different materials into SWIFT, which has primarily been designed for cosmological applications.  The high-resolution simulation of this planetary giant impact shows features that are entirely absent from previous smaller calculations. For instance after the initial impact, the subsequent sloshing around of the system can be sufficiently vigorous to unbind some of the icy mantle from the proto-Uranus – an effect that was not seen in earlier million-particle simulations. Determining what material was placed into orbit as a result of this impact is crucial for deciding if the currently observed Uranian satellites could have formed from this debris. Again, the higher resolution simulations now allow us to make reliable predictions for this. These represent the first scientific results with the SWIFT code, and demonstrate how it will open up new ways to constrain models of planetary giant impacts relevant to most bodies in the solar system.

This calculation used the Dell 940 system with four Intel skylake with 112 cores.  This system was provided by Intel and Dell as a proof-of-concept system as part of the Dell/Intel and DiRAC@Durham industrial engagement at list price of approximately £100k.

---

The FABLE simulations of galaxies, groups and clusters

Clusters of galaxies collapse from the densest regions of the early Universe to become the largest objects in existence by the present day. Their unparalleled sizes make them powerful tools to probe the growth of large-scale structure and further our understanding of cosmology. At the same time, clusters represent unique astrophysical playgrounds in which galaxies interact with each other and with the intracluster gas, which itself is churned up by the powerful energy output of supermassive black holes. This process, known as AGN feedback, stems from black holes at the centres of galaxies swallowing gas, converting it into outflows of energy and impacting not only their galaxy, but also the surrounding intracluster medium millions of light-years away.

The inclusion of AGN feedback in cosmological hydrodynamical simulations of galaxy formation has proven to be a key ingredient for producing realistic galaxies in a representative region of the universe. Projects such as the Illustris simulation have met with great success in this area and thus have made great strides in our understanding of galaxy formation and evolution, as well as guiding observations. Yet whilst this model matches a wide range of observed galaxy properties, it struggles to reproduce other important properties of galaxy clusters and groups such as the contents of the intracluster medium.

With the FABLE project we have built upon the success of Illustris by improving the models for feedback from AGN and stars so that we can simultaneously produce realistic galaxies, groups and clusters, thereby completing our picture of the Universe across this vast range of scales. Consisting of high-resolution simulations of more than thirty groups and clusters performed with the state-of-the-art moving-mesh code AREPO, FABLE reproduces a number of important observables such as the build-up of stars in galaxies and the thermodynamic properties of the intracluster medium.

FABLE further reproduces a range of integrated cluster properties and their interdependencies as represented via the oft-studied cluster scaling relations. These relations are vital for the cosmological interpretation of observed cluster data, the availability of which is set to be greatly expanded by numerous ongoing and future surveys, for example with observations of the Sunyaev-Zel’dovich effect via SPT-3G and ACTpol or in X-rays with the eROSITA and Athena observatories. Yet as larger cluster surveys beat down the statistical uncertainties, we are increasingly limited by our incomplete understanding of cluster physics and its impact on the cluster scaling relations, which are required to link cluster observables to their fundamental masses.

Cosmological hydrodynamical simulations are in a unique position to quantify these effects and thereby facilitate the use of clusters as precise cosmological probes. Indeed, careful comparison of our results with observed scaling relations hints towards a non-negligible bias in observational studies that attempt to measure the masses of clusters from their X-ray emission under simplifying assumptions about the clusters’ dynamical state.

Furthermore, the increased size and depth of future surveys means that many new clusters will lie at vast distances from which their detectable emission may have taken many billions of years to reach us. This means it is vital to understand how the cluster scaling relations change over time so that they can be applied reliably to the most distant objects in upcoming datasets. While current observations are severely limited in their ability to constrain such evolution, simulations allow us to study the scaling relations for a well-defined sample of objects over the whole history of the Universe. Upcoming results from FABLE demonstrate a significant departure of the evolution of the relations from the simple model typically assumed in observational studies, which has important implications for interpreting the incoming wealth of new cluster data.

---

Fast and energetic AGN-driven outflows in simulated dwarf galaxies

Dwarf galaxies, the least massive galaxies in the Universe, are the smallest building blocks of the cosmic structure formation process. They are interesting cosmological probes, as there are a number of apparent discrepancies when comparing observations of dwarf galaxies with the predictions of structure formation models. For example, the large disparity in the number of observed dwarfs with respect to the predicted number of dark matter haloes that may be hosting these systems (the missing satellite problem) or the lack of observed massive dwarf galaxies (the too-big-to-fail problem). Some have suggested that these issues could be resolved by changing our models of dark matter, which provides the skeleton for cosmic structure formation. Others argue that the focus should be on improving our models of the luminous, ‘baryonic’ matter that we can observe directly.

One of the most important baryonic processes, reionization, stems from the radiation produced by the first galaxies in the Universe. This radiation shuts down star formation in the low-mass dwarf galaxies rendering them undetectable. For the more massive dwarf galaxies it has been suggested that violent stellar explosions – supernovae (SNe) – could shut down star formation, however this remains controversial. Most simulations indicate that another ingredient would be needed for the theoretical models to explain the sparsity of massive dwarf galaxies in our Universe.

Recently, it has been proposed that active galactic nuclei (AGN) – actively-accreting massive black holes – might be this missing ingredient. Previously, theorists did not include AGN in their simulations of dwarf galaxies as these were only observed in more massive galaxies. However, the systematic analysis of optical surveys has revealed a population of dwarf galaxies hosting AGN, which have been confirmed by X-ray follow-up observations. All of these AGN are extremely luminous relative to their black hole mass and are likely only the tip of the iceberg. There has also been some first evidence of direct AGN feedback from the MaNGA survey, which identified six dwarf galaxies that appear to have an AGN that is preventing on-going star formation. These observational results are very exciting as they suggest that AGN may be able to solve the too-big-to-fail problem.

It is therefore timely to study the physical properties of dwarf galaxies, in particular whether the presence of an AGN can affect their evolution. Using the moving mesh code arepo, we have investigated different models of AGN activity, ranging from simple energy-driven spherical winds to collimated, mass-loaded, bipolar outflows in high resolution simulations of isolated dwarf galaxies hosting an active black hole. Our simulations also include a novel implementation of star formation and mechanical SN feedback.

Figure 1 shows the visuals of our simulation suite. Here, edge-on projections of surface density, temperature, vertical velocity and metallicity at t = 300 Myr are plotted. The control runs with neither SN nor AGN feedback (’NoFeedback’ runs) do not produce any outflows. With added SN feedback, we obtain relatively cold and slow-moving outflows. When we also include AGN activity in our simulations, the outflows are much faster as well as much more energetic, reaching temperatures of up to 10 9 K.

These results provide tantalizing evidence that the observed correlation between large kinematic offsets and AGN activity may be due to high-velocity AGN-driven outflows (see Figure 2 for an example of the different line-of-sight velocity maps observed with MaNGA). In our simulations, the same correlation arises because only the AGN-driven gas outflows are fast enough to alter the velocity patterns set by the rotational motion of the galaxy. The relatively slow-moving supernova-driven outflows on the other hand do not offset the rotational motion of the gas so that it stays aligned with the stars (see Figure 3, left panel). This then results in a significant difference between the kinematic position angles (PAs) of stars and gas for dwarf galaxies hosting AGN (see Figure 3, right panel).

We also investigated the effects of these AGN-driven outflows on star formation and found that AGN activity has a small but systematic effect on the central star formation rates (SFRs) for all set-ups explored, while substantial effects on the global SFR are only obtained with strong SNe and a sustained high-luminosity AGN with an isotropic wind. Our findings from the isolated galaxy set-up therefore indicate that AGN are unlikely to directly affect global dwarf SFRs. However, in realistic cosmological environments inflows are known to be important especially for
high-redshift dwarfs. It is hence possible that AGN-boosted outflows may prevent some of this cosmic ‘pristine’ gas reaching the dwarfs in the first place, providing a mechanism for indirect star formation regulation. This possibility will be addressed in a follow-up study using cosmological zoom-in simulations.

---

Horizon: Studying the origin and amplification of magnetic fields in galaxies

Magnetic fields permeate the Universe from cosmological to molecular cloud scales, but neither their origin nor the mechanism of how they grow in amplitude is pinned down. We have used DiRAC to make headway on this on two fronts. Firstly in Martin-Alvarez et al. 2018, MNRAS, 479,3348 we sought to understand whether hierarchical galaxy growth could amplify extremely weak (~ 10^-19 Gauss) primordial magnetic fields to micro-Gauss levels in present day galaxies. For this purpose, we ran a suite of constrained transport magnetohydrodynamical adaptive mesh refinement cosmological zoom simulations with different stellar feedback prescriptions. We found that cosmological accretion at high redshift excited a galactic turbulent dynamo, which amplified the magnetic field even in the absence of stellar feedback. Our convergence tests are the first in the literature to demonstrate the capture of the elusive turbulent dynamo in cosmological simulations of galaxies.

Secondly, in Katz et al. 2019, MNRAS, 484, 2620, we presented an algorithm to trace magnetic fields produced by different sources, ranging from primordial fields generated for example by inflation to magnetic fields injected by supernovae and black holes. The figure above displays the magnetic energy density, color-coded by origin (green corresponds to primordial magnetic fields, red denotes fields injected by supernovae, and blue indicates the interaction of the two) in one of our simulation outputs of a 150 kpc on a side cosmological volume at redshift 4. With this algorithm, we will be able to methodically explore the contributions of different sources to the magnetization of the Universe.

---

Turbulence, Shocks and Dissipation in Space Plasmas

In situ observations from solar and heliospheric spacecraft missions measure power spectra of the solar wind plasma and electromagnetic field fluctuations which show a power-law behaviour over several decades in frequency. This is interpreted as the manifestation of a turbulent cascade of energy from large fluid scales to small kinetic scales, characteristic of the motion of the ions and the electrons. Direct numerical simulations of turbulent plasmas are not only important for interpreting the nature of the solar wind fluctuations at large scales, but also a fundamental tool to understand how energy is channelled to protons and electrons. In particular, the hybrid particle-in-cell model, which treats the ions as particles and the electrons as a charge-neutralizing massless fluid, has proved to be the most appropriate and efficient method to follow the turbulent cascade down to sub-ion scales, retaining ion kinetic effects [1].

We have used DiRAC resources to perform a collection of state-of-the-art two-dimensional (2D, Fig. 1 left panel) and three-dimensional (3D, Fig. 1 middle) hybrid simulations of plasma turbulence, investigating the processes responsible for coronal heating and solar wind heating and acceleration. We have analysed, among other things, the spectral properties (Fig. 1, right), the spectral anisotropy [2], the intermittency properties, and the nature of fluctuations at kinetic scales. We also investigated the interplay between plasma turbulence and the proton firehose instabilities in a slowly expanding plasma, performing the first 3D hybrid simulation employing the expanding-box model, to take into account the solar wind expansion [3].

Our numerical results are being used to interpret in situ observations in the Earth’s magnetosheath and in the solar wind from the Magnetospheric Multiscale (MMS) mission [4], as well as to predict upcoming results from Parker Solar Probe (PSP) and Solar Orbiter.

More recently, we have produced high-resolution high-cadence 2D and 3D data exploring the conditions encountered by PSP on its first perihelion. Analysis of the results is ongoing and it will provide a large amount of information on many different aspects of plasma turbulence, to be directly compared with present and upcoming spacecraft observations.

[1] Franci et al., Solar Wind Turbulent Cascade from MHD to Sub-ion Scales: Large-size 3D Hybrid Particle-in-cell Simulations, ApJ 853(1):26, 2018
[2] Landi et al., Spectral anisotropies and intermittency of plasma turbulence at ion kinetic scales, submitted to ApJL, 2019
[3] Hellinger et al., Turbulence vs. fire hose instabilities: 3-D hybrid expanding box simulations, submitted to ApJ, 2019
[4] Franci et al., Predicting and interpreting data of solar wind turbulence with numerical modelling: hybrid simulations vs. MMS observations, in preparation, 2019

---

Longest ever simulation of a planet in a protoplanetary disc

Ragusa et al (MNRAS 474,460: 2018) analysed the longest ever simulation of a planet in a protoplanetary disc, using the FARGO-3D code (Ben´ıtez-Llambay & Masset 2016) to follow a planet’s evolution over a record 300000 orbits. During this time the planet’s eccentricity (blue in figure below) underwent a complex evolutionary variation as it exchanged angular momentum with the disc. The simulations demonstrate that a relatively small (factor three) change in the disc mass has a major effect on the final eccentricity of the planet; the complex evolutionary pattern can be understood in terms of a superposition of precessing eccentric modes of the planet + disc system. These findings are relevant to explaining the significant eccentricity of the first hot Jupiter observed in a protoplanetary disc.

Rosotti & Juhasz (MNRAS 474,L32; 2018) proposed a novel way to probe the vertical temperature profile in protoplanetary discs from measuring the pitch angle of spiral structure in discs in both the infrared and submm. Their hydrodynamic simulations, using the PLUTO code (Mignone et al. 2007) of a planet induced spiral in a disc were post-processed with the radiative transfer code RADMC-3D and were able to demonstrate that spiral structure is significantly more tightly wound in the submm than the infrared. This finding is directly linked to the fact that submm emission traces cold mid-plane regions of the disc. This work means that multi-wavelength observations of spiral structure in real discs can be used as a `dynamical thermometer’ of disc mid-plane regions.

Winter et al (MNRAS 475,552; 2018) presented long timescale Smoothed Particle Hydrodynamical simulations using the GANDALF code (co-developped in house: Hubber et al 2018; MNRAS 473,1603) of a remarkable young stellar system whose widely spaced components are associated with extremely extended far infrared emission. These simulations demonstrated that the morphology of the emission and the kinematics of the stellar components can be explained through a dynamical history involving a previous star-disc encounter. This successful reproduction of an observed system provides some of the first evidence that young stars indeed from dynamically complex few-body systems, as suggested by star formation simulations.

---

Galactic Archaeology in the Gaia Era

On 25^th April 2018, the ESA’s Gaia mission has made its second data release (Gaia DR2). Inspired by our N-body simulations of the Milky Way-like disk in this and our previous projects, Kawata et al. (2014, MNRAS, 443, 2757) and Kawata et al. (2018b, MNRAS, 473, 867), we studied the rotational and vertical velocity structure of the Galactic disk as a function of the Galactocentric radius, across a radial range of 5<Rgal<12 kpc, using the Gaia DR2 data. We find many diagonal ridge features (Fig. 1). We have found also radial wave-like oscillations of the peak of the vertical distribution in the Gaia DR2 data.

We noticed that the diagonal ridge features found in Gaia DR2 looks qualitatively similar to what we see in our N-body/SPH simulations (an example shown in Fig.2). The simulation shows transient spiral arms and we considered that the transient and winding spiral arms act as perturbers in the disk to make these ridge structures. We followed up this idea in Hunt, Hong, Bovy Kawata & Grand (2019, MNRAS, 481, 3794), and using a toy models of transient and winding spiral arms, we showed that the rotation velocity ridge features can be explained with the transient and winding spiral arms.

• New models give insight into the heart of the Rosette Nebula;
• Gravitational Waves 2017;
• Shining Light on the Structure of Hadrons;
• Stellar Structure and Nucleosyntheis;
• COSMOS Consortium: Flashes of light on dark matter;
• Virgo Consortium; New results from the Eagle Project
• Breaking resonances between migrating planets;
• Extreme QCD: Quantifying the QCD Phase Diagram;
• Hadron structure;
• UKQCD: Lattice QCD tests Standard Model at high precision;
• Modelling black hole spins and jets in the new era of gravitational wave astronomy;
• The formation of massive stars;
• Radiative Transfer Simulations for Type I Supernovae;
• Understanding our Milky Way;
• Protoplanetary disks;
• Excited Charmonia and Charmonium Resonances from Lattice QCD;
• Exomol;
• Galactic Archaeology in the Gaia Era;
• Ab initio calculation of the potential bubble nucleus 34Si;
• Understanding the origin of strong galactic outflows.

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 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.”

“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.

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.

---

 

---

 

Stellar Structure and Nucleosyntheis

Movie showing the evolution of the velocity magnitude of a 1024³ 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

GRChombo goes public!

---

 

Consortium

 

2014 was another productive year for our PPAN theory community researchers and some of the Science Challenges they tackled are highlighted below.

• UKMHD Consortium: The origins of magnetism in the quiet sun
• Kinetic Plasma Turbulence
• Leicester: Merging Supermassive Black Holes
• UKQCD Beyond Standard Model
• HPQCD: Searching for new physics in the magnetic moment of the muon
• ExoMol
• CLOVER: Baryon Mass Splittings
• Black Holes
• VIRGO
• Extreme QCD: Towards Quantitative Understanding of the QCD Phase Diagram
• COSMOS
• Domain Wall
• ECOGAL

2014 list of publications

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

• HPQCD: Stress Testing the Standard Model
• Kinetic Plasma Turbulence: Determining the shape of Type Ia supernova explosions
• HORIZON: Galaxy alignments from cosmological hydrodynamics simulations
• Leicester: First pre-infall mass estimate of a Milky Way satellite galaxy
• Black Holes
• VIRGO: Modern cosmological test of gravity
• Exeter

2015 list of publications

---

 

Stress Testing the Standard Model

• Gravitational Waves: DiRAC simulations play a key role in gravitational-wave discovery
• Lattice QCD and the Search for New Physics
• Post-periapsis Pancakes: Sustenance for Self-gravity in Tidal Disruption Events
• Molecular Data to Model (Hot) super-Venuses
• COSMOS Consortium
• Understanding the Properties and Stability of Matter
• Black Holes as the Nemesis of Galaxy Formation
• Black Holes: Illuminating the Early Universe with Radiative Transfer Simulations
• Hadron Spectroscopy: Baryon Mass Splittings
• HORIZON: Galaxy Formation Simulations
• Heavy Elements: Linking the Nuclear Interaction to the Structure of Heavy Elements

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---

 

---