Uncategorized

Understanding and obtaining predictions from Quantum Chromodynamics (QCD), the theory of strongly-interacting quarks and gluons, is one of the foremost problems in particle physics.

The vast majority of ordinary strongly-interacting matter is composed of protons and neutrons, themselves composed of light up and down quarks. However, the spectrum of QCD also contains excitations involving other types of quarks such as charm and strange. The lowest mass QCD states appear to be made of either three quarks, like protons and neutrons, or a quark-antiquark combination, like pions and kaons. Other combinations could exist, such as four-quark states or states with intrinsic gluons. These are thought to be unstable resonances that decay to lighter strongly-interacting particles and this obscures their nature.

Short-lived hadron resonances containing a charm quark are readily observed at hadron physics experiments such as LHCb. They are also accessible theoretically, such as in lattice QCD, a first-principles approach to QCD where calculations are performed numerically, and so are a particularly useful comparison point.

In this work [Phys. Rev. Lett. 129, 252001 (2022); arXiv:2205.05026], we used lattice QCD to study the scattering of a pion (composed of a light quark and a light antiquark) with a spin-1 (vector) D* meson (composed of a charm and an antilight quark), employing the DiRAC Data Analytic facility in Cambridge for some of the calculations. The intrinsic spin of the D* meson can combine with orbital angular momentum in different ways to produce states with the same total angular momentum. This produces a rich, dense spectrum of states with differing properties.

Working with an unphysically-large light-quark mass to suppress multihadron channels, we computed various scattering amplitudes and from these obtained the masses, widths, couplings and thus decay modes of two axial-vector D_1 mesons and a spin-2 tensor D_2 resonance. The upper panel of the figure shows the scattering amplitudes relevant for the axial-vector mesons. These two mesons have very different character: one couples dominantly to a decay mode with no orbital angular momentum (S-wave) producing a broad feature (red band), while the other couples to a decay mode with orbital angular momentum 2 (D-wave) and produces a narrow peak (green band). There is very little coupling between the two channels (orange band).

Excited hadrons can be described by pole singularities at complex energies in the scattering amplitude. The positions of the axial-vector meson poles in the complex energy plane are plotted in the middle panel of the figure. The couplings obtained from these poles are illustrated by the histograms at the bottom of the figure and demonstrate again that each couples dominantly to one decay mode. The tensor resonance was found to decay to D pi and D* pi in approximately equal amounts.

A useful theoretical comparison point is the limit where the charm quark is made infinitely heavy. In this limit four mesons arise in two doublets: the two in one doublet decay via S-wave and the other two decay via D-wave. Even though the charm quark is only 10 or so times larger than the scale of QCD interactions, we found that our results closely follow the pattern predicted from this limit. This study is the first time that it has been possible to study this in a first-principles approach to QCD.

Andrew Hillier, Ben Snow, Giulia Murtas

The solar atmosphere is a complex and chaotic plasma environment, rife with explosive phenomena

driven by the complex magnetic structure undergoing topological changes and releasing this

stored energy through a process known as magnetic reconnection. This process readily occurs in

the twisted (kinked) magnetic flux ropes that exist throughout the solar atmosphere, however the behaviour

of this process has not been well studied in the lower solar atmosphere where the medium is

partially ionised, i.e., consists of both ionised and neutral species. Partial ionisation can accelerate

the rate at which energy is released through magnetic reconnection and lead to additional heating.

In this project, we performed on the DiRAC COSMA7 system 3D multi-fluid, partially ionised simulations

of the kink-instability to understand how magnetic reconnection occurs in the lower solar

atmosphere. Numerical simulations are performed using the (PIP) code that evolves two-fluid (neutral,

ion+electron) equations using the non-dimensional form. The species are coupled through both

thermal collisions and ionisation/recombination. The energy required to ionise the system exists as

an energy sink in the plasma species and is balanced initially by a heating term

Our initial conditions are a 3D force-free magnetic flux rope that is unstable to the kink instability.

As the system evolves, the magnetic field twists, creating current sheets that undergo magnetic reconnection,

converting magnetic energy to kinetic and thermal energy. Since magnetic reconnection

is a plasma process, it does not directly affect the neutral species. The plasma couples to the neutral

species through thermal collisions and ionisation/recombination. We find that in a partially ionised

plasma, the instability is accelerated increasing the explosive nature of the dynamics.

Three-dimensional partially ionised simulations become numerically taxing due to the small timestep

(determined by the interactions between ionised and neutral fluids) and the complex physics

involved (thermal collisions, ionisation, recombination, radiative losses). As such, 3D simulations

of two-fluid plasma become computationally expensive. This project was only possible due to the

generous computational time provided on the COSMA7 system (DiRAC@Durham) as part of the

STFC DiRAC HPC Facility (ST/P002293/1, ST/R002371/1 and ST/S002502/1, ST/R000832/1).

Project dp048 – PI I. Baraffe – Data Intensive Cluster System in Leicester & Cambridge

Despite being the most observed star in the Universe, many open questions remain regarding the properties of convection and waves in the deep interior of the Sun. Sophisticated numerical simulations can unveil some important properties and provide a better understanding of the Sun’s convection and internal waves (Vlaykov et al. 2022; Le Saux et al 2022). Performed with the DIRAC facilities, our numerical simulations of the Sun’s convective envelope and stable, radiative central region (see image) revealed in particular a local heating process resulting from downward plume penetration below the Sun’s convective envelope. This process can have important impact on the Sun’s structure and help solving a well known problem, namely the sound-speed discrepancy between solar models and the observed structure of the Sun inferred from helioseismology (Baraffe et al. 2022).

The conventional picture of giant planet migration in protoplanetary discs involves an inexorable inspiral towards the star through the process of Type II migration, thus providing a popular scenario for the creation of hot Jupiters. In Scardoni et al 2022[1] we explored very long term (0.6M planetary orbits)[2] simulations of Type II migration in the limit that the planet is relatively massive compared with the local disc and uncovered the fact that there are situations where planets migrate outwards or else stall. The reason is that in this (`light disc’) limit, the planet migration time is relatively long compared with the viscous flow time in the disc and there is therefore time for the disc to achieve a condition close to a steady state interior to the planet. In this steady state the disc edge interior to the planet self-adjusts to deliver a (positive) torque on the planet which depends only on the accretion rate through the disc, irrespective of the thermal properties of the disc. On the other hand, the disc exterior to the planet is not in a steady state and therefore the (negative) torque that it imposes on the planet depends on the structure of the inner edge of the exterior disc, being smaller in the case of the wide and deep gaps associated with colder (geometrically thinner) discs.

In practice this means that when planets enter the inner disc (which is geometrically thinner) the balance of torques shifts towards there being a net positive torque on the planet and, in consequence, outward migration. The change in sign of planet migration as a function of the gap width parameter, K (itself a function of the geometric aspect ratio) is shown in the left hand panel below, while the right hand panel depicts how this sign change would affect planet migration over timescales of Myr.  It is notable that this effect is expected to drive planetary migration towards the location in the disc where the torques from the inner and outer disc are in balance. The location of this attractor is expected to be at a few A.U. from the star, depending on the disc temperature.

While future calculations will be required to discover whether this effect persists in 3D simulations, these results have profound consequences for the orbital architecture of exoplanetary systems, explaining the prevalence of gas giants at a few A.U. from their host stars. A corollary of this result is that hot Jupiters would not be expected to form via disc mediated migration but would instead likely result from planet scattering at a later evolutionary stage.

---

[1] MNRAS 514,5478

[2] These calculations made use of the CSD3 and DIaL services

The image shows a vertical cross-section taken from a 3D stellar simulation run on the DiRAC Memory Intensive Service. It is possible to see the turbulent motions of the plasma inside the star and how the chemical composition is transported and mixed, represented by the abundance of neon in colour scale (red more abundant, blue less abundant). This snapshot shows how the turbulent motions brings neon-rich material (dark red plumes) from the top stable region into the turbulent region via a process called entrainment.

Running and analysing stellar multi-dimensional simulations help us understand hidden aspects of important astrophysical processes, such as the evolution and fate of stars, their explosions and collapse to neutron starts and black holes, but also the stellar production of chemical elements (called nucleosynthesis) and their later expulsion to form new stars, planets and other celestial objects.

Fionntán Callan, Stuart Sim, Christine Collins, Luke Shingles, Joshua Pollin

Astrophysics Research Centre, Queen’s University Belfast, Northern Ireland

Type Ia supernovae (SNe Ia) are extremely luminous astrophysical events that arise from the thermonuclear explosions of white dwarf stars. They play a number of key roles in astrophysics including acting as cosmological distance indicators, contributing substantially to cosmic nucleosynthesis and injecting energy in galaxy evolution. Despite the significant interest in SNe Ia, the mechanisms which produce them are still poorly understood with many different explosion scenarios proposed to explain them. Modern transient surveys have also made it clear that SNe Ia are a diverse population with multiple explosion scenarios likely required to explain the various sub-classes of SNe Ia. In our DiRAC project we carry out radiative transfer computer simulations to produce predictions of light curves and spectra for various theoretical models of SNe Ia. This allows us to make direct comparisons between theoretical models and observed SNe Ia in order to better understand both the physical conditions in multiple sub-classes of SNe Ia and ultimately their origins.

Type Iax supernovae (SNe Iax) are the largest sub-class of SNe Ia estimated by Foley et al. (2013) to make up ~ 30% of the total SNe Ia rate. SNe Iax show spectroscopic differences to “normal” SNe Ia, have lower peak magnitudes and show a much larger spread in their brightness, differing by more than 4 magnitudes at peak. Previous studies have shown pure deflagration models of Chandrasekhar mass carbon-oxygen white dwarfs reproduce many of the observed characteristics of SNe Iax, making them the most promising explosion scenario to explain SNe Iax. We recently carried out a parameter study of such pure deflagration models (Lach et al., 2022) and found our model sequence was able to produce good agreement with observed SNe Iax over a wide range of luminosities, with poor agreement only found with the faintest SNe Iax. However, as was the case in previous studies (e.g. Fink et al., 2014) there are still some systematic differences between pure deflagration models and observed SNe Iax. In particular, model light curves and spectra produce good agreement with observed SNe Iax at early times but at later times the light curves decline too quickly in comparison to observed SNe Iax. This systematic difference becomes increasingly prominent when comparing to fainter SNe Iax and in redder filter bands.   

Late time observations of SNe Iax have suggested there is a luminous remnant left behind after the explosion (e.g. Foley et al., 2014). Such a luminous remnant is also predicted by pure deflagration explosion models including those from our parameter study (Lach et al., 2022). However, previous studies have not included radiation from this luminous remnant in full radiative transfer simulations. In our project we have recently been investigating whether including the radiation emitted by such a luminous remnant can help explain the systematic differences between our pure deflagration models and observed SNe Iax. From the figure it is clear that including the radiation predicted for the luminous remnant in our radiative transfer simulations significantly improves agreement between our models and observed SNe Iax, for both bright and intermediate luminosity SNe Iax. This strengthens the case that SNe Iax originate from deflagration explosions which do not fully unbind the white dwarf and leave behind a luminous remnant. This work will be submitted for publication in the coming weeks (Callan et al., 2023, in prep). We also intend to carry out future work focussing on how our radiative transfer simulations are impacted by a more detailed treatment of the remnant in the explosion simulations.

This work made use of 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:

Foley et al., 2013, ApJ, 767, 57

Lach et al., 2022, A&A, 658, A179

Fink et al., 2014, MNRAS, 438, 1762

Foley et al., 2014, ApJ, 792, 29

Phillips et al., 2007, PASP, 119, 360

Barna et al. 2021, MNRAS, 501, 1078

 A key aim of the worldwide particle physics programme is to find evidence of new physics beyond our current Standard Model (SM) that would allow us to develop a more complete theory of fundamental physics. B meson decays are a good place to look because some of these are very rare in the SM but the presence of new particles could boost their rates. The process in which a B meson (containing a b quark) decays to a K meson (containing an s quark) and a lepton/anti-lepton (electron, muon or tau) pair is a good example. It must proceed in the SM via a loop made of W bosons and top quarks, as in Fig. 1, and is highly suppressed. New particles could shortcut this loop and give a very different rate (smaller or larger depending on the combination with the SM process). Theorists in the HPQCD collaboration have been spearheading the international effort to calculate B meson SM decay rates from lattice QCD. Our efficient method for handling quarks on a spacetime lattice makes the DiRAC CSD3 supercomputer at Cambridge ideal for our calculations and enables us to achieve world-leading accuracy. 

Fig. 2 below shows our results (arXiv:2207.13371, 2207.12468) for the fraction of 𝐵+ that decay to 𝐾+𝑙+𝑙− compared to the experimental data. The LHCb experiment at CERN has the most accurate data, and we see that there is a significant difference between their values and the lattice QCD calculation, especially at low values of 𝑞2, where 𝑞2 is the invariant mass-squared of the 𝑙+𝑙− pair. Between 𝑞2 = 1 and 6 GeV2 the difference exceeds 4 times its uncertainty, which could be a hint of new physics in this decay process. The low 𝑞2 region is challenging for lattice QCD because the K meson has large momentum there. HPQCD’s calculation represents the first time that this has been successfully tackled. As we can see in Fig. 2, the low 𝑞2 region is important because experimental results are often better there. 

HPQCD was also able to predict the SM rate for 𝐵 → 𝐾𝜈𝜈̅ (𝜈 is a neutrino). This process has not been seen by experiment yet but should be visible at the BelleII experiment at superKEKB in Japan in future. It provides further exciting opportunities for new physics searches. 

Matthew Bate (University of Exeter)

The distribution of stellar masses, known as the initial mass function (IMF), is of central importance in astrophysics, due to the fact that the radiative, chemical and mechanical feedback from a star depends strongly on its mass. Despite many observational studies, there is little evidence for variation of the IMF in our Galaxy (Bastian et al. 2010 ARA&A 48 339). During the past decade, radiation hydrodynamical simulations of star cluster formation have been able to produce stellar populations that match the typical stellar properties of Galactic stellar populations quite well (Bate 2012 MNRAS 419 3115). In agreement with observations, they have found that the distributions of stellar masses produced by present-day Galactic star formation are surprisingly invariable. For example, star cluster formation calculations have shown that varying the metallicity of star-forming gas between 1/100 and 3 times solar metallicity has little effect on stellar properties for present-day star formation (Myers et al. 2011 ApJ 735 49; Bate 2014 MNRAS 442 285; Bate 2019 MNRAS 484 2341).  The only effect of metallicity for present-day star formation that has yet been identified is on the frequency of close binary systems – a triumph of the simulations of (Bate 2014, 2019), performed using DiRAC, is that they do produce the recently observed anti-correlation of close binary frequency with metallicity (Badenes et al. 2018 ApJ 854 147; El-Badry & Rix 2019 482 L139; Moe, Kratter & Badenes 2019 ApJ 875 61). 

The strongest evidence for variation of the IMF is from stellar populations that formed much earlier in the Universe (Smith 2020 ARA&A 58 577). Recently, DiRAC has been used to perform the first radiation hydrodynamical simulations of star cluster formation at high redshift: z = 5 (Bate 2023 MNRAS 519 688). The main result from this study is that although the IMF does not depend significantly on metallicity for present-day star formation (z = 0), it does depend on metallicity at redshift z = 5: high metallicity gives a `bottom-light’ IMF in which low-mass stars are much rarer than in present-day star formation (Fig. 1). The difference at high-redshift is that the cosmic microwave background radiation is much warmer. High-metallicity gas is unable to cool to as lower temperatures at z = 5 as at z = 0, resulting in less fragmentation and larger typical stellar masses.

Figure 1: The cumulative stellar initial mass functions (IMFs) obtained from radiation hydrodynamical simulations of star cluster formation performed at redshifts z = 0 (left; Bate 2019) and z = 5 (right; Bate 2023). The IMFs at z = 0 are quite insensitive to the metallicity of the molecular gas (ranging from Z = 0.01 − 3 Z_⊙), but at redshift z = 5 the IMFs become metallicity dependent, with high metallicity gas producing a bottom-light IMF.

High-redshift surveys have so far discovered over a hundred quasars above redshift of six. This number will likely increase significantly in the coming years, due to ongoing and planned deep, wide-field surveys, such as eROSITA in X-rays, the Vera Rubin Observatory at optical wavelengths as well as Euclid in infrared. Such observations will also push these discoveries to lower luminosities too, giving a more complete picture of the build up of the black hole population in the early Universe. In fact it could even be possible to detect 10⁵ solar mass black holes at z = 10 with mega second exposures from NIRCam on JWST or from the proposed Lynx X-ray telescope.

Focusing on the most extreme examples, supermassive black holes in excess of billion solar masses have been detected above z = 7, challenging theoretical models of the growth of such objects. The current record holder in terms of luminosity, and hence black hole mass is SDSS J010013.02+ 280225.8, with an inferred mass exceeding 10¹⁰ solar masses at z = 6.3. While there have been a number of theoretical studies of these objects, they still struggle to form some of the most massive observed black holes at z > 6.

As a part of our DiRAC project dp012 (Bennett et al., MNRAS, to be submitted), using the computational facilities in Cambridge and Durham, we have investigated the most promising scenarios for building up the most massive known supermassive black holes in the early Universe. Allowing for mildly super-Eddington accretion and earlier seeding redshift which is still entirely compatible with the direct collapse seeding model, we have found that it is possible to assemble a 10¹⁰ solar masses black hole by z = 6.

As shown in Fig. 1 the mass growth of our simulated black hole is consistent with observations, hinting that episodic super-Eddington accretion may be required in the early Universe to grow the most extreme mass black holes within a Gyr of cosmic time. Importantly we found that the feedback impact of this black hole is very significant.

This is shown in Fig. 2 where we plot the combined thermal and kinetic SZ decrement (top) and gas radial velocity (bottom) for the original FABLE model (left) and our new simulations (right). Hot, fast outflows are pushing gas way beyond the virial radius, diminishing the SZ signal on small scales and enchanting it on large scales. Moreover, by studying the redshift evolution of the simulated hot halo we make detailed predictions for mock maps, from radio to X-rays, that should help us disentangle how mass is accumulated onto these gargantuan black holes over cosmic time e.g. an early (sustained) growth and feedback scenario vs. a late and rapid growth and feedback scenario.

New discoveries in particle physics can be made in different ways. One can let particles collide at higher and higher energies and look for direct evidence of new particles and interactions, or search instead for tiny but measurable discrepancies between experimental measurements and theoretical predictions, at the so-called precision frontier. This project is pushing the latter by studying with high precision the decays of mesons like pions and kaons, which are the lightest particles made of two quarks and gluons, into leptons. This is a process where the fundamental weak interaction allows the annihilation of the two constituent quarks and the consequent emission of an electron-neutrino or muon-neutrino pair (the muon being a heavier copy of the electron).

The decay of these particles is dominated by the effect of the strong interactions, also called QCD, inside the meson, but the newly aimed-for precision makes it necessary to take into account also other subleading effects. This is the case of those from QED interactions, which result in the exchange of photons between any of the charged elementary particles involved in the process – see e.g. Figure 1, where the decay of a charged pion into muon and neutrino is depicted.

One way to compute QCD-governed processes at the conditions of our current universe is to simulate them on a finitely sized four-dimensional grid or lattice, corresponding to a discretised and finite version of space and time. In our calculation we have been able to simulate particles with masses that match the values as found in nature. This is generally very challenging, and it has been possible thanks to recent algorithmic developments and the DiRAC high-performance computing resources. Photons, however, are massless particles and make the calculation more difficult: since they can travel large distances before interacting with other particles, simulating QED interactions in a finite box will result in sizeable finite-volume effects. The DiRAC Extreme Scaling facilities allow us to simulate QCD+QED on very large grids, but even with this cutting-edge technology the finite-volume effects are still large and have to be studied carefully. In our finite-volume simulation we are able to compute with very high precision the rate of the decay of the light mesons into muon-neutrino pairs including QED corrections, but when comparing with experimental results – which are performed in our vastly larger universe – the final precision of our results is dominated by these finite size effects. This is illustrated in Figure 2, where we show our result for the QED correction to the ratio of the kaon and the pion decay rates, called δR_Kπ. The point on the right is the one we obtain at the volume of our simulation through an extensive data analysis, while the point on the left shows how much larger our uncertainty gets when we try to extrapolate our result to the “infinite” volume limit. This is due to our limited knowledge of how this quantity varies when changing the size of the box. In red, green, and blue we show some scaling trajectories on which our determination of δR_Kπ could lie on, although with our current study we cannot discriminate between them. This means dedicated simulations on multiple and larger boxes are required to reduce the uncertainty on our infinite volume estimate.

Once we will reduce such a large uncertainty with further simulations on DiRAC resources, the precision of our simulations will become comparable to that of experimental measurements, and we will be able to test the validity of our current theory of particles and interactions with unprecedented accuracy.