ABSTRACT
We outline our techniques to characterize photospheric granulation as an astrophysical noise source. A four-component parameterization of granulation is developed that can be used to reconstruct stellar line asymmetries and radial velocity shifts due to photospheric convective motions. The four components are made up of absorption line profiles calculated for granules, magnetic intergranular lanes, non-magnetic intergranular lanes, and magnetic bright points at disk center. These components are constructed by averaging Fe i 6302 Å magnetically sensitive absorption line profiles output from detailed radiative transport calculations of the solar photosphere. Each of the four categories adopted is based on magnetic field and continuum intensity limits determined from examining three-dimensional magnetohydrodynamic simulations with an average magnetic flux of 200 G. Using these four-component line profiles we accurately reconstruct granulation profiles, produced from modeling 12 × 12 Mm2 areas on the solar surface, to within ∼ ±20 cm s−1 on a ∼100 m s−1 granulation signal. We have also successfully reconstructed granulation profiles from a 50 G simulation using the parameterized line profiles from the 200 G average magnetic field simulation. This test demonstrates applicability of the characterization to a range of magnetic stellar activity levels.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
One of the consequences of plasma motions in the outer layers of low-mass stars with convective envelopes is radial velocity (RV) shifts due to variable stellar line profile asymmetries, known as astrophysical noise (or stellar "jitter"). This can pose a major problem for planet hunters as RV follow-up is mandatory for most planet confirmation and characterization. As these net RV shifts, produced by photospheric convective motions, are on the m s−1 level (Schrijver & Zwaan 2000), stellar granulation will severely impact the confirmation of low-mass planets around solar-like stars. Additionally, the promise of cm s−1 RV precision from spectrographs like Espresso and Codex makes it ever more pertinent that we understand granulation both in terms of the physical processes involved as well as an important source of astrophysical noise.
Meunier et al. (2010) conclude that the limiting factor for detecting Earth-like planets in the habitable zone is, in fact, the convective contribution to the RVs and not the larger amplitude contribution of spots and plages. Meunier et al. (2010) go on to affirm that the amplitude of convective noise in solar-like stars must decrease by one order of magnitude or more to make possible the detection of low-mass planets.
To date, there have been only a handful of attempts to decrease the effects of granulation for the benefit of exoplanet detection and confirmation. Dumusque et al. (2011) propose adjustments to the number of measurements and exposure time in order to average out some of the effects of granulation. Although this method has had success (see Pepe et al. 2011 and Dumusque et al. 2012), it is both time and cost intensive, and does not provide information on the nature of granulation. Aigrain et al. (2012) take a different approach and use photometry to estimate RV variations due to stellar activity, including a parameter to estimate the perturbation due to convection. They successfully demonstrate that they can reproduce variations well below the m s−1 level. However, as low-mass planets around solar-like stars require cm s−1 precision, there is still an urgent need to further characterize and remove the RV variations of surface convection.
Although there have been few attempts to remove or reduce the effects of granulation, there have been numerous observations and theories devoted to understanding the nature of granulation (for a comprehensive review see Nordlund et al. 2009, Stein 2012, and references therein). We were particularly inspired by the four-component model of Dravins (1990), used to reconstruct the observed asymmetries in stellar line profiles due to granulation. Similar to Dravins (1990), our aim is to measure and understand the RV signal produced by convective motions integrated across a stellar disk. However, our main goal is to simulate numerous granulation observations and ultimately develop a statistically significant noise removal method.
In this paper, the first in a series, we outline a technique to characterize granulation as an astrophysical noise source. The backbone of this characterization is a state-of-the-art three-dimensional (3D) magnetohydrodynamic (MHD) solar simulation, coupled with detailed wavelength-dependent radiative transfer. Due to the time-intensive nature of detailed radiative diagnostics, producing enough realistic granulation patterns to cover an entire stellar disk with this method is not feasible. As such, we use the output from a short time series of simulations to parameterize the granulation signal.
In Section 2, we outline our procedure to separate absorption line profiles from solar simulations into four categories based on continuum intensity and magnetic field strength to categorize the different physical components present in photospheric magneto-convection. In Section 3, we use these components to reconstruct the absorption line profiles from the simulated models. In Section 4 we show that this parameterization can be successfully applied to other magnetic field strengths. We discuss our procedure to generate new realistic granulation line profiles in Section 5 and present the success of our parameterization as well as its current limitations in Section 6.
2. CHARACTERIZATION OF CONVECTION COMPONENTS
The radiative MHD simulations of solar magneto-convection reproduce the properties of granulation observed at the solar surface in detail and for a wide range of magnetic activity levels (for a comparison of simulation and observation, see Figure 1, and for further comparisons with observables see Shelyag et al. 2007). However, there is no established to method to remove the granulation signal from observational data. We use such simulations, examining the Sun as a star in this paper, to develop a parameterized model of granulation that will eventually be used to remove convective noise from observations and can later be applied to other spectral types. Applicability to other spectral types is essential to developing a robust noise removal technique for exoplanet surveys due to the variety of planet-hosting stars. The simulations also provide many advantages toward developing such a technique including controllability of magnetic flux, seeing/resolution, and rotation.
We simulate solar surface granulation using the MURaM code (Vögler et al. 2005). The code solves large-eddy radiative 3D MHD equations on a Cartesian grid, and employs a fourth-order central difference scheme to calculate spatial derivatives. The numerical domain has a physical size of 12 × 12 Mm2 in the horizontal direction and 1.4 Mm in the vertical direction and is resolved by 480 × 480 × 100 grid cells (see Shelyag et al. 2011 for details).
We first introduce a uniform vertical 200 G magnetic field and produce a sequence of 190 snapshots of solar magnetic granulation, with a ∼30 s cadence. The sequence covers approximately 80 minutes of physical time, corresponding to ∼10–20 granular lifetimes. This does not include the first 20 minutes of the simulation, which is discarded, in order to allow the large amplitude oscillation, induced from the introduction of magnetic field, to dampen and the model to stabilize (Vögler et al. 2005; Vögler & Schüssler 2007). During this first part of the simulation, the magnetic field is advected into the intergranular lanes by the convective plasma motions, and the maximum field strength rises from its initial value of 200 G to a few kilogauss in the intergranular lanes. The 200 G field is sufficiently strong to ensure each physical aspect of convection found in the quiet Sun is present in the simulations.
We also produce simulations in a similar manner for a uniform vertical 50 G field. The cadence again is ∼30 s, and this simulation covers approximately 45 minutes of physical time. A similar oscillation in the 50 G model is present throughout the simulation since the lower magnetic field does not provide as much damping compared to the 200 G model. This weak magnetic field simulation provides a valuable test case for our parameterized model. An accurate reconstruction of the 50 G model from the 200 G components demonstrates the robustness of our parameterization across magnetic field strengths.
The calculated models contain basic solar photospheric plasma parameters such as density, pressure, temperature, magnetic field strengths, and velocity. We perform the radiative transport diagnostics with these parameters, using a modified version of the STOPRO code (as described in Shelyag et al. 2007), for the simulated photosphere with the Fe i 6302 Å magnetically sensitive absorption line in local thermodynamic equilibrium approximation. This line is widely used in both solar observations and simulations as a diagnostic tool for magnetic field and temperature (Shelyag et al. 2007). The full Stokes vector is calculated for each pixel in the snapshots throughout the simulation; for this study we use only the I component of Stokes vector. The wavelength range used for calculation is ±0.3 Å and is resolved by 400 wavelength points.
2.1. Four-component Model
The premise for our parameterization of granulation line profiles is that they are a summation of absorption line profiles corresponding to each physical component of solar surface granulation. As such, we began our parameterization by examining the different physical components of photospheric magneto-convection in the 200 G simulations. Granules are non-magnetic, hot, bright bubbles of rising gas. The plasma in granules cools, darkens, and falls back down into the intergranular lanes. The magnetic bright points (MBPs) present within the dark intergranular lanes are due to the enhanced continuum intensity and decreased radiation absorption in evacuated, highly magnetic flux tubes. Thus, a four-component parameterization, which divides the granulation pattern into granules, magnetic intergranular lanes, non-magnetic intergranular lanes, and MBPs, is a natural step to move forward.
The simulation provides information on magnetic field and continuum intensity (an example of which is shown in the left panel of Figure 1), alongside the Stokes profiles. Since the four components differ inherently in brightness and magnetic field, it is logical to impose cutoff limits based on continuum intensity and magnetic field (measured at the continuum formation height) to separate each component within the simulation. As a result, for each snapshot we average together the Stokes I profiles that correspond to each category to create our component profiles. Thus, for each of the 190 snapshots in the time sequence there are four line profiles, each representing a component within the granulation pattern.
Figure 1. Example of our simulation of solar granulation (left) in the white-light continuum (417 nm) compared to observations (right). The observations were taken with the Rapid Oscillation in the Solar Atmosphere (ROSA) imaging system on the Dunn Solar Telescope. The contrast in the observations is reduced due to seeing effects.
Download figure:
Standard image High-resolution imageTo determine the optimal magnetic field and continuum intensity limits for the four components we explored the parameter space from 100 to 1200 G and from 0.8 to 1.2 times the average continuum intensity. The final selection of limits, 0.9 average continuum intensity and 1000 G, was based on how well the four profile components could reconstruct the original simulated line profiles; the procedure of reconstruction is discussed in Section 3. The four profiles used to characterize the different components of photospheric magneto-convection are shown in Figure 2. The breadth about each component's respective means is representative of the precision obtained. The characteristic bump seen in the MBP and magnetic intergranular lane profiles is due to the Zeeman splitting from the high magnetic fields present in these components.
Figure 2. Line profiles from the entire time series for the different physical components of granulation: granules (green), MBPs (red), magnetic (blue), and non-magnetic (purple) intergranular lanes. The time-average line profiles are shown in black. All components have been shifted by their bisector mean wavelength to the bisector mean wavelength of their time-average profile to remove the oscillation signal. The breadth about each component's respective means is representative of the precision obtained.
Download figure:
Standard image High-resolution image2.2. Removing Oscillations
Throughout the time sequence, all profiles experience an RV shift. An examination of Fourier power in the pressure domain of the simulations reveals two peaks at around 2 and 7 mHz. We attribute these variations to the solar p-modes that contribute to the well-known "5 minute" oscillations. In order to disentangle the granulation signal, this oscillation needs to be removed from the profiles. Removing the oscillation allows the convective motions to be viewed in detail and prevents time-average component profiles from being skewed by oscillation-induced broadening, thus allowing a more accurate granulation parameterization.
In order to remove the effects of the oscillation, we shifted each component profile to fixed reference points. The reference points were determined by first creating initial time-average component profiles. The mean values of the line bisector wavelengths from these initial time-average profiles served as the reference points because they had no knowledge of the oscillation. We then shifted the components by their bisector mean wavelengths to preserve shape variation due to the granulation, while simultaneously removing the oscillation effects. Once the profiles were shifted to remove the oscillation, we then used a final time average for each component to be representative of each granulation feature. As seen in Figure 2, there is little variation within each component category throughout the time sequence, providing confidence that the oscillation signal does not significantly affect component profile shape and that it has been successfully removed.
3. PARAMETERIZED RECONSTRUCTION
We examined the accuracy of the four time-average line profiles, representing the different components of photospheric magneto-convection, by reconstructing the line profiles produced by the simulations. For each granulation snapshot, we recorded the proportion of the area occupied by each component to the total area (filling factor) based on the number of Stokes I profiles in a given category to total number in the snapshot; mean values for these within 1σ are shown in Table 1. To reconstruct the line profile for an entire snapshot, we multiplied each component profile by its corresponding filling factor and then added all four-component profiles together. Examples of the reconstructions are shown in Figure 3.
Figure 3. Reconstructed line profiles (red) are compared to the original line profiles (black) from the 200 G simulations. Left: worst reconstruction achieved; right: best reconstruction achieved.
Download figure:
Standard image High-resolution imageTable 1. Mean Filling Factor for Each Component with ±1σ
| Component | Mean Filling Factor ±1σ |
|---|---|
| Granule | 0.559 ± 0.016 |
| Non-magnetic intergranular lane | 0.348 ± 0.019 |
| MBP | 0.054 ± 0.004 |
| Magnetic intergranular lane | 0.039 ± 0.006 |
Download table as: ASCIITypeset image
We calculated the average relative error, 〈ΔI/〈I〉〉 (where I is the continuum intensity), between the reconstructed profiles and the original profiles produced by the simulations. Even at the very worst, the reconstructed and original profiles still agree to at least ∼2%, as seen in Figure 3. On average the reconstructed and original profiles are in agreement to ∼0.6% or better.
3.1. Impact of Additional Components
We have also investigated the impact of assigning additional categories to our parameterization scheme. Upon examining the probability density distribution of the modulus of the vertical component of magnetic field (measured roughly at the continuum formation layer in the simulation box) versus the continuum intensity, shown in Figure 4, we realized a potential for an additional two components. There are two higher density clusters, one at high and one at low magnetic fields with a range of continuum intensities. These clusters indicate the presence of magnetic and non-magnetic features that are both bright and dark. In our four-component model we separate these regions based on a magnetic field limit of 1000 G and a continuum intensity limit of 0.9. However, it is clear from Figure 4 that many profiles with intermediate magnetic fields were being included in the non-magnetic components (granules and non-magnetic intergranular lanes). We chose to test a model where there was an additional magnetic field limit at 200 G to separate the non-magnetic features from those with intermediate magnetic field. This additional limit then introduced two new components that we term semi-magnetic bright and semi-magnetic dark regions, shown in Figure 5, with magnetic field between 200 G and 1000 G. We found that the six component model did not produce significantly smaller average relative errors when compared to the original line profiles from the simulations. In fact, at times the four-component model even faired better, as seen in the top panel of Figure 6. Even at the most extreme, the average relative errors between four and six component model differed by only ∼0.0005 (<1%). For simplicity we proceed with the four-component model.
Figure 4. Probability density distribution of the modulus of the vertical component of the magnetic field approximately at the continuum formation level vs. continuum intensity for an arbitrarily chosen snapshot from the 200 G simulation. The average dependence of the continuum intensity on the modulus of the vertical magnetic field strength is shown by the black dot-dashed curve. The cutoff limits for each of the physical granulation components are shown in black (∼0.9 continuum intensity, 1000 G for both four and six component models and, additionally, 200 G for the six component model).
Download figure:
Standard image High-resolution image
Figure 5. Line profiles from the entire time series for the different physical components of granulation, as in Figure 2, with two additional components: granules (in green), MBPs (in red), magnetic (in blue) and non-magnetic (in purple) intergranular lanes, and semi-magnetic bright (in turquoise), and semi-magnetic dark (in pink) regions. As before, oscillation effects have been removed and the breadth about each component's respective means (shown in black) is representative of the precision obtained.
Download figure:
Standard image High-resolution image
Figure 6. Top panel: the average relative error between the original line profiles from the simulations and the reconstructed line profiles over time. Middle panel: the RVs of the reconstructed line profiles over time. Four-component model is shown in black and six component model is shown in red. Bottom panel: filling factor vs. time for the two largest components, granules (in green) and non-magnetic intergranular lanes (in purple) in four-component model. Any correlation between RV and average relative error is likely due to changes in filling factors that simultaneously impact the RV and line shape of a profile.
Download figure:
Standard image High-resolution image3.2. Impact of Errors
Ultimately we want to measure the RV signal produced by the convective motions integrated across a stellar disk. In order to do so, we must be confident in our ability to accurately reproduce the RVs from these convective motions using our four-component parameterization. We need to ensure that any resultant RV variations are not due to errors introduced in the reconstruction process. To do so we began by measuring the RVs of the reconstructed profiles.
We measured the RVs from reconstructed profiles by cross-correlating a single, randomly selected profile with the entire time series. At times, there is a possible correlation between average relative error (top panel, Figure 6) and RV (middle panel, Figure 6). We believe this is due to the relationship between filling factor and RV. For example, if the granule profile has a larger error, compared to the profiles from the other components, then a reconstruction with a higher granule filling factor would have a more negative RV and also a larger error. We see this in Figure 6, where changes in component filling factor (bottom panel, shown in green for granules and in purple for non-magnetic intergranular lanes) correlate with RV (middle panel) and error (top panel). However, such a correlation could also indicate that the measured RVs may originate from errors within the reconstruction process rather than the granulation signal.
To confirm that the measured RVs for the reconstructed profiles do originate from the granulation signal, and not from reconstruction errors, we compared the reconstructed RVs to those of the original line profiles. To measure the RVs of the original profiles we cross-correlated them with a randomly selected reconstructed profile. We used a reconstructed profile, arbitrarily chosen, as a reference because it has no knowledge of the oscillation (first panel, Figure 7). As mentioned in Section 2.2, we removed the oscillation for our parameterization and as a result we must also remove it from the original profiles to directly compare with the reconstructed profiles. To do so, we began by estimating the oscillation signal. If the reconstruction profiles match perfectly the simulated granulation then a cross-correlation of the original profiles with their respective reconstructed profiles would provide the oscillation signal alone. While the reconstruction process is not perfect, we assume it is good enough to give a good approximation of the oscillation signal (second panel, Figure 7). This approximate oscillation signal is then subtracted from the original profile RVs, and therefore any remaining RV signal should be due to the convective motions and any reconstruction errors (third panel, Figure 7). (We note here that if the oscillation signal were to significantly alter the component profile shape, then remnant signatures would be present after this stage and in the reconstructions, which we do not find.) We then compared this remaining RV signal to the reconstructed RVs (fourth panel, Figure 7), and calculated the difference between them to produce the residuals (fifth panel, Figure 7). The residuals are approximately ±20 cm s−1 for more than 95% of the time sequence. As such, we conclude that our categorization is good to ±20 cm s−1 on a ∼100 m s−1 granulation RV variation calculated over a 12 × 12 Mm2 area. We therefore conclude that any errors in the parameterization do not significantly impact our ability to generate line profiles that represent realistic granulation patterns.
Figure 7. First panel: RVs of the original profiles over the time series; the oscillation is clearly evident in the first half of the series. Second panel: RVs of the approximate oscillation signal. Black and red curves represent use of the four and six component reconstructions, respectively. Third panel: original profile RVs after the approximate oscillation signal has been subtracted. Fourth panel: RVs of the reconstructed profiles. Fifth panel: difference between original and reconstructed RVs, after removal of oscillation. As shown in the fifth panel, our reconstruction is good to ∼ ±20 cm s−1 on a ∼100 m s−1 granulation signal (green dashed lines are plotted only to guide the eye).
Download figure:
Standard image High-resolution image4. 50 G RECONSTRUCTION
Our four-component model parameterizes granulation line profiles from the 200 G simulations, where the full spectrum of photospheric magneto-convection features is present. In accordance with our premise (Section 2.1), such a parameterization should be capable of accurately reproducing granulation line profiles throughout the quiet Sun, where average surface magnetic fields vary from near 0 up to ∼200 G. We confirm this by using our four-component 200 G parameterization to reconstruct granulation profiles produced from the simulation with a 50 G average field.
The 50 G simulation was parameterized in the same manner as the 200 G (see Section 2.1). From the simulation and parameterization, we found the filling factors of each of the four granulation components in a 50 G model. As such, we use these filling factors in conjunction with the four-component profiles from the 200 G parameterization to reconstruct the profiles from the 50 G simulations. As shown in Figure 8, we have achieved excellent accuracy in our reconstruction, ranging from at worst ∼1% to at best ∼0.3%. From these results, we are confident that our parameterization of Stokes I profiles is robust for the various magnetic field strengths found in the quiet Sun.
Figure 8. Four-component line profiles from 200 G magnetic field simulation are used in conjunction with filling factors from 50 G simulation to reconstruct granulation line profiles in a 50 G field. The reconstructed line profiles (in red) are compared to the original line profiles (in black). Left: worst reconstruction achieved; right: best reconstruction achieved.
Download figure:
Standard image High-resolution image5. GENERATING NEW GRANULATION LINE PROFILES
To create new stellar line profiles that represent realistic granulation patterns over a 12 × 12 Mm2 patch, we generate a probability distribution for the granules based on filling factors throughout the time sequence. This allows us to randomly and realistically select the granule filling factor for a given profile. As seen in Figure 9, there is a strong linear anti-correlation (with slope −1.19) between granule and non-magnetic intergranular lane filling factor (with a Pearson's correlation coefficient of −0.98). As a result, the granule filling factor is used in conjunction with this anti-correlation to determine the filling factor for the non-magnetic intergranular lanes. Together, granules and non-magnetic intergranular lanes make up approximately 90% of the total granulation profile (Table 1). Since there are no other strong correlations amongst the filling factors, a probability distribution is again generated, this time for the magnetic intergranular lane filling factor. As the four-component filling factors must total to one, the MBP filling factor is then determined by the difference between the other three filling factors and one. Once all four filling factors are determined we use these in conjunction with our four-component profiles determined above to create new realistic granulation line profiles.
Figure 9. There is a strong anti-correlation between granule filling factor and the non-magnetic intergranular lane filling factor. The solid red curve is a linear fit. This fit is used in conjunction with a randomly selected granule filling factor to determine the non-magnetic intergranular lane filling factor used to generate new stellar line profiles representative of realistic granulation patterns.
Download figure:
Standard image High-resolution image
6. CONCLUDING REMARKS
Simulating realistic solar granulation can take up to months for a sequence of snapshots that last little more than an hour. Computing the Stokes profiles from these simulations for just one snapshot can take hours on its own right. Currently, it is not feasible to tile a stellar disk with each 12 × 12 Mm2 tile having an independent granulation pattern produced by simulations. However, using our parameterized model described here we can produce accurate granulation line profiles for magnetic field scenarios between 50–200 G. In addition, this method reproduces RVs which are accurate to ∼ ±20 cm s−1 on a ∼100 m s−1 granulation signal.
Throughout this paper we outlined the methodology for our parameterization of granulation at disk center. Naturally, the granulation component characterization may change as we apply our parameterization to various disk positions; however, our methodology will remain the same. Applying our parameterized method to off-center regions of the stellar disk is more difficult. Photospheric magneto-convection results in a corrugated surface. An apt analogy is to think of the granules as "hills" and the intergranular lanes as "valleys" (Dravins 2008). Beyond zero degrees, we begin to see different aspects of these "hills" and "valleys." As we incline the simulations the granular walls become visible, while the intergranular lanes become ever more hidden, and MBPs disappear from view toward the limb. Consequently, there are changes in profile shape for granulation patterns at various inclinations. Furthermore, higher inclinations lead to a physically longer line of sight and require finer resolution (in terms of grid cell number) along the inclined line of sight due to the appearance of strong under-resolved gradients in plasma parameters. Thus, in order to maintain the same precision, an even greater amount of computational time is required. In a future paper, we will extend our parameterization across the stellar disk using the methodology outlined here.
Once we have a model that applies to various stellar disk positions, we will tile a stellar disk with line profiles that represent random, independent granulation patterns and integrate the wavelength-dependent intensity over the disk. These "simulated star" observations will allow us to search for a way to characterize the granulation signal as a whole by examining relationships among line asymmetries, bisector changes, and RV. As this method stems from the radiative 3D MHD simulations, we plan to expand these simulations to other spectral types thereby allowing our parameterized method to be tested beyond the solar spectrum. In future work, we will also explore such a parameterization and search for noise removal with other stellar absorption lines—we begin here with one line to validate our method. There are a number of potential routes for achieving RV detection of Earth-like planets. Currently, RV measurements of exoplanets are done via cross-correlation with a spectral type template mask, requiring several thousand lines. However, it may be that for Earth-like confirmation we may have to target specific, well-behaved unblended lines, in which case we may only need a handle of lines modeled. Additionally, if we do find correlations between granular noise and observables then we may be able to search for said correlations empirically for a variety of lines, which may limit the need for such detailed modeling of other stellar lines.
With our current parameterized four-component model, we are able to accurately match the synthesized Stokes I profiles obtained from radiative diagnostics of 3D MHD radiative simulations of the solar photosphere and mimic observations of solar granulation, while maintaining the opportunity to apply our model to other spectral types. As we move ever closer toward the confirmation of Earth-like worlds, granulation and other low-level noise sources such as meridional flows or variable gravitational redshift (Beckers 2007; Makarov 2010; Cegla et al. 2012) become an imminent limitation. While our primary goal is to remove/reduce granulation effects in the exoplanet searches, a characterization of granulation is multidisciplinary and may be beneficial for solar physics and stellar astrophysics. Furthermore, as we expand our approach to other spectral types we will not only be able to assist exoplanet surveys around a larger variety of host stars, but also understand the convection process more thoroughly across the H-R diagram.
This research has made use of NASA's Astrophysics Data System Bibliographic Services. H.M.C. acknowledges support from a Queen's University Belfast university scholarship. C.A.W. acknowledges support by STFC grant ST/I001123/1. Financial support from the Leverhulme Trust (grant RPG-249) is gratefully acknowledged. We thank the anonymous referee for their comments which have significantly helped improve the clarity of this paper.