## Abstract

Observations of the redshift *z* = 7.085 quasar J1120+0641 are used to search for variations of the fine structure constant, α, over the redshift range 5:5 to 7:1. Observations at *z* = 7:1 probe the physics of the universe at only 0.8 billion years old. These are the most distant direct measurements of α to date and the first measurements using a near-IR spectrograph. A new AI analysis method is employed. Four measurements from the x-shooter spectrograph on the Very Large Telescope (VLT) constrain changes in a relative to the terrestrial value (α_{0}). The weighted mean electromagnetic force in this location in the universe deviates from the terrestrial value by Δα/α = (α* _{z}* − α

_{0})/α

_{0}= (−2:18 ± 7:27) × 10

^{−5}, consistent with no temporal change. Combining these measurements with existing data, we find a spatial variation is preferred over a no-variation model at the 3:9σ level.

## INTRODUCTION

### Theoretical background

What fundamental aspects of the universe give rise to the laws of Nature? Are the laws finely tuned from the outset, immutable in time and space, or do they vary in space or time such that our local patch of the universe is particularly suited to our own existence? We characterize the laws of Nature using the numerical values of the fundamental constants, for which increasingly precise and ever-distant measurements are accessible using quasar absorption spectra.

The quantity that we focus on here is the fine-structure constant, which can be expressed as α = *e*^{2}/4πε_{0}ħ*c*, where *e* is the electron charge, ε_{0} is the permittivity of free space, ħ is the reduced Planck constant, and *c* is the speed of light. The dimensionless quantity described by α is the ratio of the speed of an electron in the lowest-energy orbit of the Bohr-Sommerfeld atom to the speed of light. α may be considered to relate quantum mechanics (through ħ) to electromagnetism (through the remaining constants in the ratio).

The quest to determine whether the bare fine-structure constant, α, is a constant in space and time has received impetus from the recognition that there might be additional dimensions of space or that our constants are partly or wholly determined by symmetry breaking at ultrahigh energies in the very early universe. The first proposals for time variation in α by Stanyukovich (*1*), Teller (*2*), and Gamow (*3*) were motivated by the large-number coincidences noted by Dirac (*4*, *5*) but were quickly ruled out by observations (*6*). This has led to an extensive literature on varying constants that is reviewed in (*7*–*11*).

There are also interesting new problems that have been about extreme fine-tuning of quantum corrections in theories with variation of α by Donoghue (*12*) and Marsh (*13*). Accordingly, self-consistent theories of gravity and electromagnetism, which incorporate the fine-structure “constant” as a self-gravitating scalar field with self-consistent dynamics that couple to the geometry of spacetime, have been formulated in (*14*–*20*) and extended to the Weinberg-Salam theory in (*21*, *22*). They generalize Maxwell’s equations and general relativity in the way that Jordan-Brans-Dicke gravity theory (*23*, *24*) extends general relativity to include space or time variations of the Newtonian gravitational constant, *G*, by upgrading it to become a scalar field. This enables different constraints on a changing α(*z*) at different redshifts, *z*, to be coordinated; it supersedes the traditional approach (*25*) to constraining varying α by simply allowing α to become a variable in the physical laws for constant α. Further discussions relating spatial variations of α to inhomogeneous cosmological models can be found in (*26*, *27*).

Direct measurements of α are also important for testing dynamical dark energy models, since they help to constrain the dynamics of the underlying scalar field (*11*), and thus, dynamics can be constrained (through α) even at epochs where dark energy is still not dominating the universe. The possibility of doing these measurements deep into the matter era is particularly useful, since most other cosmological datasets (coming from type Ia supernovas, galaxy clustering, etc.) are limited to lower redshifts.

The inputs needed for these theories come from a variety of different types of astronomical observations: high-precision observations of the instantaneous value of α(*z _{i}*) characterizing quasar spectra at various redshifts

*z*to test possible time variations, both local and nonlocal measurements of

_{i}*28*–

*30*) to search for spatial variation, the cosmic microwave background (

*31*,

*32*), the Oklo natural reactor (

*33*–

*36*), atomic clocks (

*7*,

*37*,

*38*), compact objects in which the local gravitational potential may be different (

*39*,

*40*), or atomic line separations in white dwarf atmospheres (

*41*).

These detailed studies, using large numbers of quasar spectra, hint at a spatial variation of α at a level of ∼4σ (*42*–*44*). The apparent spatial signal was initially claimed to be caused by long-range wavelength distortions (*45*), but a more detailed analysis showed this to be incorrect (*44*). This apparently persistent spatial signal motivates further direct measurements, especially by extending the measurement redshift range.

### Observations and artificial intelligence algorithm

The relative wavelengths of absorption lines imprinted on the spectra of background quasars are sensitive to the fine-structure constant. Comparing quasar measurements with high-precision terrestrial experiments provides stringent constraints on any possible spacetime variations of the fine-structure constant, as predicted by some theoretical models (*7*, *11*, *16*, *46*, *47*). The quasar J1120+0641 (*48*) is of particular interest in this context because of its very high redshift. Its emission redshift is *z* = 7.085, corresponding to a look-back time of 12.96 billion years in standard ΛCDM (lambda cold dark matter) cosmology. J1120+0641 is one of the most luminous quasars known (*49*), enabling high spectral resolution at high signal-to-noise ratio. We make use of spectra obtained using the x-shooter spectrograph (*50*) on the European Southern Observatory’s Very Large Telescope (VLT), with nominal spectral resolution *51*). The total integration time is 30 hours.

The x-shooter instrument provides a broad spectral wavelength coverage. This maximizes the discovery probability of absorption systems along the sightline, enabling the identification of potential coincidences (i.e., blends) between absorption species at different redshifts, an essential step in making a reliable measurement of α. In all, 11 absorption systems are detected (Table 1). Desirable characteristics of an absorption system are a selection of transitions with different sensitivities to a change in α and a velocity structure in the absorbing medium that is as simple as possible.

Of the 11 absorption systems identified along the J1120+0641 sightline, 4 are found to be suitable for a measurement of α, at redshifts *z*_{abs} = 7.059,6.171,5.951, and 5.507. The atomic transitions used to measure α in these four systems are highlighted in Table 1. The highest-redshift system has, of the four, the least sensitivity to varying α. No other direct quasar absorption α measurements have previously been made at such high redshift. Before the measurements described in this paper, the highest-redshift quasar absorption direct measurement of α was at *z* = 4.1798 (*52*). Voigt profile models for each of the four absorption systems were automatically constructed using a genetic algorithm, gvpfit, which requires no human decision-making beyond initial setup parameters (*53*). The genetic part of the procedure controls the evolution of the model development. vpfit (*54*) is called multiple times within each generation to refine the model, which then becomes the parent for subsequent generations. Absorption model complexity increases with each generation. A description of gvpfit can be found in (*53*), where it was used for the analysis of an absorption system at *z*_{abs} = 1.839 toward the quasar J110325−264515. That particular system had previously been analyzed by several groups and so provided important comparative information between gvpfit and previous methods. Further assessments of gvpfit’s performance are given in (*55*). The procedure outperforms human interactive methods in that it gives objective, reproducible, robust results and introduces no additional systematic uncertainties. The method is computationally demanding, requiring supercomputers. New procedures have been introduced for the analysis in this paper, beyond those described in (*53*), and so are described here.

The analysis of each of the four absorption systems took place in four stages. Throughout, Δα/α is kept as a free-fitting parameter, making use of the Many-Multiplet Method (*28*, *38*). In stage 1, we imposed the requirement that all velocity components are present in all species being fitted, irrespective of line strength. Without this requirement, an absorbing component in one species might fall below the detection threshold determined by the spectral data quality, but not in another. This requirement was only applied in the first stage because it was found in practice to help model stability by discouraging the fitting procedure from finding a model with physically implausible cloud parameters in one or more components. By “physically implausible,” we mean either large linewidth, i.e., a Voigt profile *b* parameter of tens of kilometers per second, or an improbably high column density. The requirement is dropped subsequently. gvpfit was allowed to evolve (that is, the complexity of the model was allowed to increase) for the number of generations required to pass through a minimum value of the corrected Akaike Information Criterion statistic (AICc) (*56*, *57*). The model resulting from this first stage of the analysis is the model at which AICc is at a minimum and is already quite good but is not final.

In stage 2, we use the model from stage 1 as the parent model input to gvpfit but now drop the requirement that all velocity components are present. The other requirements from stage 1 were carried over to stage 2. At this stage, one further increase in model complexity is introduced. Although the spectral continuum model was derived before the line fitting process, we allow for residual uncertainties in continuum estimation where needed by including additional free parameters allowing the local continuum for each region to vary using a simple linear correction as described in the vpfit manual (http://www.ast.cam.ac.uk/~rfc/vpfit11.1.pdf). The minimum AICc model from this stage is again taken as the parent model for the next stage.

In stage 3, we check to see whether any interloping absorption lines from other redshift systems may be present within any of the spectral regions used to measure α. When interloper parameters are introduced, degeneracy can occur with other parameters associated with the metal lines used to measure α. To avoid this problem, all previous parameters are temporarily fixed, and gvpfit is used in a first pass to identify places in the data where the current model is inadequate. Interlopers, modeled as unidentified atomic species, are added automatically by gvpfit to improve the current fit.

In stage 4, the model resulting from this third stage is used as the input model for the fourth and final part of the process, which entails running gvpfit again but this time with all parameters free to vary (subject to the physical constraint that all *b* parameters are tied, and all redshifts of corresponding absorbing components are tied, as was the case throughout all stages).

## RESULTS

In previous non–artificial intelligence (AI) analyses, the general approach was to construct absorption system models based on turbulent broadening (*30*) and then to build a thermal model from the turbulent parameters. One important advantage of the AI approach is that it is straightforward to construct turbulent and thermal models independently, and this has been done for all four absorption systems reported here. While this was possible before gvpfit automation, it was very time consuming to do manually and therefore was generally not done. The Doppler or *b* parameters of different ionic species are related by *i*th ionic species has mass *m*, *k* is Boltzmann’s constant, and *T* is the temperature of the absorption cloud. The first term describes the thermal contribution to the broadening of the *b* parameters, and the second term describes the contribution of bulk, turbulent motions. If the linewidths for a particular absorption cloud are dominated by thermal broadening, the second term of the equation is zero and vice versa if the broadening mechanism is predominantly turbulent. These two cases are the limits of possibilities for the values of the *b* parameters.

We have modeled each absorption system using the two limiting cases: first assuming that the lines are thermally broadened and then assuming turbulent broadening. Modeling in this way results in two measurements of the fine-structure constant for each absorption system. Table 2 gives the results, which show that both measurements, for all four absorption systems, are consistent with each other. Rather than discarding the highest χ^{2} model, since both models are statistically acceptable, we give a single value of α from our results using the method of moments estimator to determine the most likely value. The method of moments estimator compares the weighted relative goodness-of-fit differences between the thermal and turbulent models. This method is conservative, in that it only ever increases the uncertainty estimate on α from the smallest one and accounts for cases where the fits are consistent (like our results), and where inconsistent, the value chosen is more heavily weighted to the model with a lower χ^{2} [see (*30*)].

Figure 1 illustrates one model for the lowest-redshift system analyzed (see caption for details). All final model parameters associated with the four high-redshift absorption systems modeled here are provided in the Supplementary Materials, which also provides information relating to upper limits on potential systematic effects due to wavelength distortions. This is the first time multiple absorption systems along a given sightline have been simultaneously modeled to constrain the presence and impact of long-range wavelength distortions across a large wavelength range. This is important because in this way, the distortion model parameters are more tightly constrained, and hence, the possible additional systematic error on Δα/α is minimized. We find that in this case, the additional systematic is smaller than the statistical uncertainty on Δα/α.

The x-shooter spectral resolution may not resolve individual absorbing components. However, we simultaneously fit multiple transitions at the same redshift, with tied parameter constraints, such that the Voigt profile parameters and α measurements are reasonably well constrained. Nevertheless, the lower spectral resolution of x-shooter [compared to echelle spectrographs such as UVES (Ultraviolet and Visual Echelle Spectrograph) and HIRES (High-Resolution Echelle Spectrometer)] would lead us to expect that some absorption components are missed. In a small number of cases, elevated *b* parameters in the final models reinforce that expectation (full model parameter details and estimated uncertainties for all four absorption systems are provided via the Supplementary Materials). Even so, Fig. 2 indicates that α is likely to be insensitive to missing components because α stabilizes in relatively early model generations and subsequently varies only slightly as model complexity increases. The same insensitivity of α to missing components was borne out in the numerically simulated spectral simulations described in (*55*). The impact of subsequent higher spectral resolution would evidently reduce the Δα/α error bar, but the x-shooter results presented here should not be systematically biased by the lower resolution.

The fine-structure constant has been measured in four high-redshift absorption systems using an x-shooter spectrum of the *z*_{em} = 7.084 quasar J1120+0641. These four measurements are the highest-redshift direct measurements of α to date. The final results are summarized in Table 2, giving both statistical and systematic parameter uncertainties. The weighted mean value of α is consistent with the terrestrial value and is Δα/α = − 2.18 ± 7.27 × 10^{−5}.

## DISCUSSION

Our analysis compares relative shifts between atomic species and therefore directly tests for any possible change in α. No molecular species have been used. We thus have no sensitivity to changes in the proton-to-electron mass ratio, μ = *m*_{p}/*m*_{e} or combinations of constants such as α^{2}μ and α^{2}*g*_{p}/μ; hence, no model-dependent coupling constant assumptions have been necessary.

To update the parameters associated with the spatial dipole discussed in (*42*,*30*), we form a new combined sample of α measurements as follows:

1) the four new x-shooter measurements from this paper;

2) the large sample of 293 measurements from (*30*);

3) 20 measurements from (*58*), 14 of which were remeasurements of points already in (*30*), the points from (*58*) taking priority; and

4) 21 recent measurements as compiled in (*43*).

Together, our final sample comprises a total of 323 measurements spanning the redshift range 0.2 < *z*_{abs} < 7.1. These data points are illustrated in Fig. 3. The new statistical sample defined here enables an updated estimate of the spatial dipole model reported in (*30*): the updated dipole amplitude, *A* = 0.72 ± 0.16 × 10^{−5}, and the dipole sky location is right ascension 16.76 ± 1.17 hours and declination −63.79 ° ± 10.30°.

Using the bootstrap method described in (*30*) to estimate statistical significance, this deviates from a null result at a level of ∼3.9σ. We can also directly compare the dipole model prediction (using the new parameters above) with the actual weighted mean from the four new x-shooter measurements: The dipole prediction for the weighted mean is Δα/α = − 0.19 × 10^{−5}, in agreement with the actual measurement of Δα/α = − 2.18 ± 7.27 × 10^{−5}.

The x-shooter data presented here highlight an important benefit that is generally not available with higher-resolution echelle spectra of quasars: The extended wavelength coverage increases our ability to detect absorption systems along the line of sight simply because more transitions at the same redshift appear. Systems that might otherwise remain undiscovered or uncertain become clear. This is important because potential blends in transitions of interest are revealed (Table 1), and hence, systematic effects on the measurement of α that those blends may cause can be accounted for. A second advantage of the extended x-shooter wavelength coverage is that since there are more transitions falling within the observed spectral range, a more stringent constraint on Δα/α is achieved. Ultimately, the precision of the very high redshift measurements reported in this paper will be improved by obtaining higher spectral resolution using new instrumentation such as HIRES on the ELT (Extremely Large Telescope).

## MATERIALS AND METHODS

### Observations and data reduction

Observations of the quasar J1120+0641 were obtained using the x-shooter spectrograph on the European Observatory’s VLT. The total integration time was 30 hours spanning a period from March 2011 to April 2014 (European Southern Observatory programs 286.A-5025, 089.A-0814, and 093.A-0707). All exposures were taken with slit widths of 0.9 arc sec for the visual (vis) and near-infrared (nir) arms of the x-shooter spectrograph, giving spectral resolutions of *R* = λ/Δλ = 7450 and 5300, respectively. However, inspection of telluric absorption lines indicate a higher spectral resolution, suggesting that the atmospheric seeing was better than the slit width used. Atmospheric absorption lines were measured as having an *R* ≃10,000 for the vis arm and *R*≃7000 for the nir arm, consistent with a seeing full width at half maximum ≃0.7 arc sec. These values are consistent with the more detailed discussion about x-shooter resolution given in (*59*). The total wavelength coverage is approximately 5505 to 22,740 Å. The spectral signal-to-noise varies across the spectrum and is approximately 21 per 10 km s^{−1} pixel at 11,191 Å.

Data reduction was performed using custom Interactive Data Language (IDL) routines. The procedures include flat fielding of the exposures and sky subtraction using the optimal extraction method as described in (*60*). The extracted one-dimensional spectra, rebinned to 10 km s^{−1} pixels, were flux-calibrated using response curves derived from standard stars. Absolute flux calibration was performed by scaling the corrected spectrum to match the VLT/ Focal Reducer and Low Dispersion Spectrograph 2 (FORS2) and Gemini Near-Infrared Spectrograph (GNIRS) on the Gemini North Telescope spectra of J1120+0641 obtained by Mortlock *et al.* (*48*). Atmospheric line removal was performed using SkyCalc atmospheric transmission models (www.eso.org/sci/software/pipelines/skytools/). A comprehensive description of the observations and data reduction is given in (*51*).

### Continuum fitting

Before profile fitting the absorption systems of interest, we need a reliable estimate of the unabsorbed quasar continuum. This was obtained using the iraf (http://ast.noao.edu/data/software) task continuum. Small spectral regions flanking each absorption line were selected, and the pixels containing absorption were masked. The spectral regions used to estimate the underlying continuum (by fitting cubic splines, typically of order 3) each contained ∼100 to 300 pixels.

### Identification of absorption systems

Identification of absorption systems and atomic species present was carried out using qscan (*61*), an interactive Python program to display the spectrum on a velocity scale such that transitions occurring at the same redshift align visually in velocity space. The spectral ranges chosen for fitting were selected as described in (*58*). Eleven absorption systems were detected (Table 1). Of these 11, 4 were selected for their sensitivity to fine-structure constant variation: *z*_{abs} = 7.05852, *z*_{abs} = 6.17097, *z*_{abs} = 5.95074, and *z*_{abs} = 5.50793. We excluded Si IV and C IV in the determination of α in the three latter systems as the ionization potentials of these transitions are significantly higher than the other (more sensitive) transitions available. Mg II 2796 in absorption system at *z*_{abs} = 5.50794 falls in a region of the spectrum containing an incompletely removed telluric line and was not included in the modeling. For the highest-redshift system, *z*_{abs} = 7.059, only high-ionization species were available. The lower sensitivity of these results in a substantially larger error on α but including the system has the advantage of producing a tighter constraint on any possible long-range distortion in the spectrum, hence improving the overall result.

### Atomic data and sensitivity coefficients

We used multiplets from different atomic species simultaneously to constrain any possible variation of the fine-structure constant α. The method used, the Many-Multiplet Method, was introduced in (*28*, *38*). Sensitivity coefficients (*q* coefficients, parameterizing the sensitivity of an observed wavelength to variation of the fine-structure constant) are compiled in (*62*) along with laboratory wavelengths, oscillator strengths, and hyperfine structure and spontaneous emission rates (Γ values). The large nonordered range in *q* coefficients and their different signs create a unique varying α signature and assist in overcoming simple systematic effects. Figure 4 shows how the transition wavelengths of Si II, Al II, Fe II, and Mg II (the transitions used in this analysis) depend on α. The range in α is exaggerated for illustration.

### Further details and final model adjustments

The gvpft modeling produces near-final fits, but additional physical considerations, not coded into the AI methodology, are helpful in deriving the final absorption line models. These relatively minor tweaks to the nonlinear least-squares input-guess parameters were done (i) to remove or minimize the presence of parameters that are physically implausible and (ii) to further improve the model by reducing the overall χ^{2}. The notes here record those considerations and justify final model parameters.

*z _{abs} = 5.507261*. No changes to the gvpfit models were required. The following discussion applies to both the thermal and the turbulent fits.

The gvpfit model for this system included an interloper at approximately +100 km s^{−1} in Fe II 2600 Å (see Fig. 1) for which the column density and *b* parameter were poorly constrained. Removing this interloper resulted in a negligible change to either the overall χ^{2} or any of the remaining model parameters, so the interloper has been excluded.

Visual inspection shows a very broad shallow continuum depression over the Fe II 2600-Å absorption line at approximately −80 km s^{−1}, as Fig. 1 illustrates. gvpfit modeled this using a high-*b* component (*b* = 81 km s^{−1}) straddling the whole absorption complex. There is no species identification for this interloper. It may be due to real unidentified absorption or it may be some observational artifact. The associated degeneracy in the final model due to these additional parameters is minimal and affects negligibly on the other model parameters.

*z _{abs} = 5.950744*. Thermal fit: The Mg II lines reveal three components. The gvpfit model found an interloper heavily blended with the left-hand Mg II component at approximately −20 km s

^{−1}. Cosmic ray events on the detector spoil the Fe II lines in this region, so the leftmost component provides only a very weak constraint on Δα/α, which is thus constrained almost entirely by the right-hand component at +75 km s

^{−1}. Once the initial model fits were available, it became apparent that several pixels in the Mg II 2796- and 2803-Å lines were significantly deviant. We assumed that this was due to weak cosmic rays contaminating the spectrum. We thus manually clipped three pixels in the 2796-Å line and two pixels in the 2803-Å line, as can be seen in fig. S2.

Turbulent fit: The independently derived turbulent model found by gvpfit differs to the thermal model in that no interloper was assigned to the leftmost Mg II component. While this makes the model rather different to the thermal one, the impact on Δα/α and its uncertainty is minimal, and both thermal and turbulent models yield a consistent result for Δα/α.

*z _{abs} = 6.170969*. No changes to the gvpfit model were made for this absorption system. This absorption system was modeled by gvpfit as single component. The

*b*parameters for the transitions seen in this system are comparable with those found when modeling higher-resolution data. No interlopers were identified.

*z _{abs} = 7.058521*. The CI V and N V lines are strong, visually comprising two components, but which are found by gvpfit to break into further components. Si IV is weak so contributes little to the α constraint despite being more sensitive to a change in α. No model changes were made to the automated fit.

### Checking for long-range wavelength distortions

Long-range wavelength distortions have been found and measured in the Keck/HIRES, VLT/UVES, and SubaruHDS spectrographs (*45*, *63*–*65*). While no analogous distortions have been identified in x-shooter, in this analysis, we apply caution and assume that they could be present. It has recently been shown that such distortions can be modeled independently of any additional calibration exposures (*44*). We use the same method to model a putative distortion, taking advantage of the presence of multiple absorption systems along the same line of sight. The presence of many transitions spread over a wide range in observed wavelength allows us to place tight constraints on any possible distortion.

Previous studies have found that the functional form of the long-range distortions (found in solar twin and asteroid measurements) are approximately linear, with no shifts found at the central wavelength of the science exposure (*45*, *64*). That is, they can be parameterized as a linear fit, with slope γ (m s^{−1} Å^{−1}), passing through the central wavelength of the exposure. We have adopted these assumptions for determining a best-fit distortion model for J1120+0641.

We solve externally for the slope of a simple distortion model. The slope is varied in steps of δγ = 0.05 m/s/Å about the best-fit value, and vpfit is used to solve for the absorption model parameters at each step. The overall χ^{2} is then minimized as a function of distortion slope, as described in (*44*). Using all four absorption systems simultaneously to solve for the distortion slope, we find γ = 0.15 ± 0.19 m s^{−1} Å^{−1} (fig. S8). We thus find no significant evidence for long-range distortion. We nevertheless include an additional term in the final error Δα/α budget corresponding to this γ value.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/17/eaay9672/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**Results are based on observations collected at the European Southern Observatory, Chile, programs 286.A-5025(A), 089.A-0814(A), and 093.A-0707(A). We are grateful for the award of computing time for this research on the gStar and OzStar supercomputing facilities. M.R.W. acknowledges support from an Australian Postgraduate Award. J.K.W. thanks the John Templeton Foundation, the Department of Applied Mathematics and Theoretical Physics and the Institute of Astronomy at the Cambridge University for hospitality and support, and Clare Hall for a Visiting Fellowship. We thank J. King for useful discussions.

**Funding:**The work of A.C.L. and C.J.A.P.M. was financed by the FEDER—Fundo Europeu de Desenvolvimento Regional funds through the COMPETE 2020—Operational Programme for Competitiveness and Internationalisation (POCI) and by Portuguese funds through FCT—Fundação para a Ciência e a Tecnologia in the framework of the project POCI-01-0145-FEDER-028987. A.C.L. was supported by an FCT fellowship (SFRH/BD/113746/2015) under the FCT Doctoral Program PhD::SPACE (PD/00040/2012). J.D.B. thanks the STFC for support. S.E.I.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 669253).

**Author contributions:**J.K.W. conceived and managed the project. J.K.W. and J.D.B. wrote the paper. M.R.W., V.D., D.M., and J.K.W. prepared the figures and tables. S.E.I.B. carried out the data reduction. M.R.W., C.-C.L., V.D., and J.K.W. carried out the spectral analysis. J.K.W. and C.-C.L. carried out the administration, submission, and referee communications/responses. All other coauthors contributed through discussions/ideas at various stages throughout.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).