2020 Highlights

Published by Alastair Williams on

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

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.  

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

References

Haworth et al. 2020 MNRAS 494 4130

HPQCD: LUV is in the air!

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

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

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

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

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

Spontaneous Symmetry Breaking in the 3d Thirring Model

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

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

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

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

Amplification of superkicks in black-hole binaries through orbital eccentricity

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

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

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

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

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

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

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

Exploring fundamental fields with strong gravity

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

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

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

Discovery+ Universe Unravelled and advanced in-situ visualization

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

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

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

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

The Universe Unravelled series is available on the new Discovery+ service here:
www.discoveryplus.co.uk/show/universe-unravelled-with-the-stephen-hawking-centre

Convection influenced by rotation: simulations and theory

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

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

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

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

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

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

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

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

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

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

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

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.

page1image23355312

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.

page1image23357392

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

Small-scale magnetic reconnection from an anisotropic turbulent cascade

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

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

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

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

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

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.

The Q2-dependence of the lowest moment of the nucleon’s isovector forward Compton amplitude.

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.

The ud interference (`cats ears’) contribution to the nucleon’s Compton amplitude.

In a recent paper [2], we extended this investigation and demonstrated for the first time that the Q2-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- Q2) can be significant and will have important consequences for attempts to extract the leading-twist parton distributions from experimental data at low- Q2. 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 Q2, however their size in the small-Q2 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  Q2 which clearly shows an enhancement at small Q2.

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

Figure 1: Net yield, Ynet, that gives the change of number of atoms in the grain per incident oxygen ion as a function of energy E of the incident oxygen ion. A positive yield corresponds to dust destruction; a negative yield corresponds to dust growth. When the energy of the incident oxygen ion is higher than the sputtering threshold energy, dust atoms can be sputtered. Left: trapping conditions are not fulfilled (oxygen particle escaping after the sputtering event) so the resulting net yield is equal to the regular sputtering yield, Ynet = Ysp (red line). Right: the oxygen ion is trapped so the net yield amounts to Ynet = Ysp – 1 (black line). For both cases, the incident oxygen atom is accreted at lower energies (green shaded region).
Figure 2: Surviving dust mass, M, as a function of time, taking into account grain-grain collisions (GG) and sputtering (SP) with or without trapping or accretion of oxygen ions (different colors and line types). KeV oxygen ion trapping and gas accretion reduce the dust destruction significantly.

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

Fig.1: Synthetic ALMA channel maps – “images” of the emission at a particular velocity – for the CO J=3-2 line from our simulated discs. The lower panels show the results for a disc without self-gravity, which shows the “butterfly” pattern resulting from Keplerian rotation. The upper panels show the corresponding channel maps for a gravitationally unstable disc: the spiral density waves manifest as very clear “GI-wiggles” in the channel maps. If detected in real discs, these wiggles would represent unambiguous evidence for gravitational instability.

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.

The circumbinary disc structure for a thin disc (left panel) and a thick disc (right panel). In the thin case, where the disc viscosity is small, the disc orbits are strongly affected by resonances which transfer energy and angular momentum from the binary to the disc. In the thick case, where the disc viscosity is large, the disc can reach inwards to radii where the orbits become unstable and accrete directly on to the binary with significant excess angular momentum.

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.

Figure 1 The NA62 experiment at CERN (image: CERN)

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

The figure shows the first calculation of the scattering amplitude in the singlet channel which would contribute to the properties of the Higgs boson: the lattice result and its error are shown by an orange dashed line, they suggest that scalar resonance is stable is our simulations [2]

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

The figure shows the spectrum of the Sp(2N) theory in the continuum limit for N = 1, 2, 3, 4, and N = ∞, in units of the string tension [3,4]. Continuum quantum numbers are reported at the top. For comparison, we have reported also the masses of the A++, A++* and E++ glueballs for SU(∞), which coincide with the smooth extrapolation of corresponding to Sp(∞) of the Sp(2N) states with the same quantum numbers.

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 Nf=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

 Face-on (top) and edge-on (bottom) stellar light projection of a high-resolution simulated Milky Way-mass spiral galaxy at the present day from the Auriga Project (Grand et al. (2017 MNRAS.467,179). Bluer colours indicate younger stars and redder colours older stars.

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.

Top: MW Galactic rotation curve (symbols with error bars) as a function of radial distance. The solid red line is the best fitting MW mass model assuming a contracted dark matter NFW halo. The dashed blue line is assuming no contraction. Bottom: the difference between the data and the best fitting contracted halo model. The dashed  blue line is the difference between the NFW and a contracted halo model.

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, MDM = 0.97 ± 0.20 × 1012 Mo, 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+17% 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.  

Comparison of light curve properties for deflagration models and observations of Type Ia supernovae: peak magnitude versus time to reach peak for r band. The red and black symbols show observed type Ia supernovae, which are labelled by supernova name. Black indicates spectroscopically normal type Ia supernovae while red marks members of the type Iax sub-class. Results from our sequence of new simulations are shown in blue (previous theoretical calculations from Fink et al. (2014) are shown in green, for comparison). Our simulations show that variations in ignition configuration can produce a wide range of rise times and peak brightnesses, which are comparable to the diversity exhibited by observed SNe Iax. We also note that the models reach almost a magnitude fainter at peak compared to the pure deflagration models of Fink et al. (2014), placing them in the luminosity regime comparable to the faintest members of the Iax class (e.g. 2008ha and 2019gsc). However, the light curve evolution of the fainter models remains too fast when compared to observations – we are now investigating the origin of this discrepancy as a means to place further constraints on the explosion models.

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

Representation of the AREPO Voronoi shown for the ‘base refinement’ (left) and ‘shock refinement’ (right) runs.

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.

From top to bottom: Maps of average  Mach number, temperature, metallicity and turbulent velocity in the halo. The left column shows the ‘base refinement’ run and the right column shows the ‘shock refinement’ run. The black dashed line indicates the virial radius.

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.

3D surface showing accretion shocks from two different orientations, coloured by Mach number. Cosmic filaments are shown in grey.

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.

Figure: DK scattering amplitudes (related to the scattering cross section) computed using lattice QCD shown in red (S-wave) and green (P-wave) at two light-quark masses. The rapid rise in the S-wave amplitudes is found to be due to a bound state that is closer to threshold at the lighter light-quark mass corresponding to a pion mass mp  = 239 MeV (left) compared with the heavier light-quark mass corresponding to mp = 391 MeV (right).

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.

Fig 1: The lower panel shows a temperature slice for a jet with axis perpendicular to the plane of the circumnuclear disc and launched into a high pressure CGM. The top left panel shows a profile of the specific energy of the (left) lobe of the jet, split into kinetic and thermal components. The top right panel shows a mass profile of the (right) lobe of the jet split into contributions from disc, jet and CGM material.

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.

Fig 2: Various slices for a jet that is initially directed into the circumnuclear disc. The three smaller panels on the left show slices of the temperature field, evolving in time with 3 Myr intervals from top to bottom. The main panel shows a slice of the density field after a further 3 Myrs have elapsed. The two inset plots show slices of the temperature field in the centre and the vorticity field in the outflow. In the bottom left of each panel the arrow indicates the direction of the jet and is labelled with the corresponding inclination angle to the vertical.

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.

Figure 1: A Molleweide projection of neutral hydrogen column density (top) and dust (bottom) for a galaxy that is leaking LyC photons and that has an excess of [OIII] to [OII].  Due to the inhomogeneous distribution of hydrogen and dust, an observer would only measure a photon escape fraction greater than 10% for 38% of the viewing angles even though this galaxy has a high excess of [OIII] to [OII] which is not much impacted by viewing angle.

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 (MDM = 1011.8 Msun) at redshift 6.  The simulations include a seven species non-equilibrium chemistry (including H2) coupled to radiative transfer. We reach resolutions of 1000 Msun for the stars, 4 X 104 Msun/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.

Figure 2: From left to right: Gas column density map, stellar surface-mass density map, [OIII] intensity map, [OII] intensity map, and spatial distribution of O32 (the excess of [OIII] to [OII]). The top row shows an example of a galaxy with high O32 and high photon escape fraction (fesc) while the bottom row shows an example of a similar mass galaxy with high O32 but low fesc. Despite both having high O32, structural, metallicity, stellar population, and feedback differences between the two systems mean that they have very different LyC 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.

Categories: Highlights

Related Posts

2019 Highlights

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. Read more…

2018 Highlights

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. Read more…

2017 Highlights

Some of the Research highlights from our users in 2017: 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; Read more…