Some of the Research Highlights from our users in 2018:

2018 list of publications

The weighty issue of quarks

Figure 1. An artist’s impression of an ηc meson containing a charm quark and antiquark (large blue circles) in a strongly interacting background soup of gluons (wavy lines) and quark-antiquark pairs (small circles).

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, mc, 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.

Figure 2. A comparison of lattice QCD results for mc, given in the conventional continuum QCD scheme. The top result is our new one, and the third our 2014 result, both obtained on DiRAC. The dark grey band is the lattice QCD world average, with much better uncertainty than the previous average light grey band (marked PDG).

HPQCD developed one of the main techniques used to determine quark masses (the JJ method) and most recently used this, on DiRAC, to determine mc 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 mc 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.

Figure 1: Left panel: A ‘fan’ plot for the light quark ud and us decay constants (normalised with the average ud, us decay constant) for one lattice spacing against the pion mass; Right panel: The heavy b quark ub, sb decay constants, again normalised, for several lattice spacings also plotted against the pion mass.

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.

Figure 2: The eleven gravitational-wave signals observed in the first and second LIGO-Virgo observing runs (O1 and O2).

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

Figure 3: The figure shows (from left to right) time-averaged meridional sections of the azimuthal flow, meridional stream function, kinetic helicity, radial and azimuthal magnetic fields; the right hand plots are typical snapshots of the radial field and azimuthal flow at the surface. The top row is for the case of uniform electrical conductivity; the bottom row is for a Jupiter-like profile, in which the electrical conductivity drops off sharply at 90% of the planetary radius.

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.

Figure 4: 3-D numerical simulation of en inclined planet in a protoplanetary disc, performed on the DiRAC Data Intensitve @ Leicester system. Our simulations showed that an inclined giant planet can create an observable warp in the disk. (Figure adapted from Nealon et al. 2018).

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 and

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.

A volume rendering of the bosonic energy density (blue) and the BH conformal factor (grey).

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.

A mid-collision snapshot of a grazing impact with 10$^8$ SPH particles coloured by their material and internal energy, showing some of the detailed evolution and mixing that can now be resolved. In the left panel, light and dark grey show the target’s ice and rock material, respectively, and purple and brown show the same for the impactor. Light blue is the target’s H-He atmosphere.

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.

Figure 1: The main panel shows the gas distribution of one of the largest fable clusters, revealing an array of interconnecting filaments and infalling galaxies and groups. Panel 1 shows an observed spiral galaxy with an enormous tail of gas shining in X-rays and panel 2 shows a similar view for one of our simulated galaxies. This gas has been stripped from the galaxy as it falls into the dense cluster environment. Panel 3 highlights a smaller group of galaxies of various shapes and sizes, while panel 4 shows the structure of AREPO’s moving hydrodynamical mesh through a slice of the central galaxy.

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.

Figure 2: The evolution of the fable cluster shown in the first figure, starting from just 2.2 billion years after the Big Bang (left) through to the onset of a major merger between two smaller clusters 5.7 billion years later (right) that will eventually form one of the largest possible gravitationally-bound objects in the simulated (and real) universe.

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.

Figure 1: Edge-on projections of all the simulation runs at t = 300 Myr. The first row shows surface densities, the second row shows temperatures, the third row shows vertical velocities, and the fourth row shows metallicities. The runs with low star formation efficiency are depicted on the left hand side, and the high star formation efficiency runs on the right hand side. With only SN feedback, the outflows are slow, warm, and high-density. With additional AGN feedback, the outflows have an additional hot, low-density component moving at high velocities.

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 2: Comparison between the line-of-sight velocity maps for stars and ionized gas for a galaxy without observed AGN activity (left panel) and a galaxy with AGN activity (right panel). With AGN activity, the ionized gas is visibly offset from the stellar component.

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

Figure 3: Left panel: Comparison between the 2D velocity maps of the stars and gas in the simulations at t = 150 Myr. The colour coding of the pixels indicates the surface density. The spatial coverage of the primary MaNGA sample of 1.5R eff is indicated by a white circle. With added AGN feedback, the central circular motion is significantly altered by the outflows so that the gas would be seen as kinematically offset. With SN feedback alone, the circular motion remains dominant and we do not observe this feature. Right panel: Comparison between the line-of-sight velocity maps for stars and ionized gas at t = 150 Myr. Only gas and stars within 1.5R eff are included (in analogy to MaNGA), and the spatial resolution is matched to the MaNGA resolution at the mean redshift of the primary sample. For the AGN runs, the ionized gas is visibly offset from the stellar component, just like in the observed MaNGA galaxies.

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.

Figure: The magnetic energy density at redshift 4 in a (150 kpc)3 piece of a cosmological volume. Using a new tracer algorithm (Katz et al. 2019), we have color-coded the magnetic energy density by origin where green corresponds to primordial magnetic fields, red denotes fields injected by supernovae and blue indicates the interaction of the two. Primordial magnetic fields dominate the volume of the simulation whereas supernovae-injected magnetic fields are associated with halos (Katz et al. 2019) where they may be amplified by a turbulent dynamo (Martin-Alvarez et al. 2018).

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

Figure 1. Left panel: Pseudocolor plot of the power of the perpendicular (with respect to the mean field) magnetic fluctuations when turbulence is fully developed in a 2D 4096×4096 simulation. Middle panel: The same as in the left panel, but for a 3D 512x512x512 simulation. Right panel: Omnidirectional 1D power spectra of magnetic field B, electric field E, ion bulk velocity u, parallel magnetic field Bz, and ion density n, from the 3D simulation. Power laws with different slopes represent theoretical predictions and/or observational results, as reference. In the insets, a direct comparison from 2D and 3D simulations shows a remarkable agreement.

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.

Figure 1: Eccentricity e as a function of time for light (left-hand panel) and massive (right-hand panel) case. The blue curve shows the planet eccentricity, the green curve the disc eccentricity at R = 4.7 in the light case and at R = 5 (azimuthal averages) in the massive one, while the red curve is a global measurement of the disc eccentricity starting from the AMD.

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 25th 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.

Figure 1: Normalized distribution of the rotation velocity for the selected Gaia DR2 stars with the RVS radial velocity data (upper) and our all selected Gaia DR2 stars (lower) as a function of the Galactocentric radius. The vertical dashed lines show the position of the Scutum, Sagittarius, Local, and Perseus spiral arms from left to right, which are calculated from Reid et al. (2014, ApJ, 783, 130). The vertical solid line is the assumed solar radius. The white lines highlight the identified ridge features.

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.

Figure 2: The rotation velocity vs. the galactocentric radius for the star particles in N-body/SPH simulation from Kawata et al. (2014). The particles are selected from a similar region of the galactic disk to where the Gaia DR2 data are analysed in the left panel.