Artificial intelligence and quasar absorption system modelling; application to fundamental constants at high redshift
Chung-Chi Lee, John K. Webb, R. F. Carswell, Dinko Milaković
MMNRAS , 1–14 (2020) Preprint 7 August 2020 Compiled using MNRAS L A TEX style file v3.0
Artificial intelligence and quasar absorption systemmodelling; application to fundamental constants at highredshift
Chung-Chi Lee (cid:63) , John K. Webb † , R. F. Carswell , and Dinko Milakovi´c . DAMTP, Centre for Mathematical Sciences, University of Cambridge, Cambridge CB3 0WA, UK School of Physics, University of New South Wales Sydney, NSW 2052, Australia Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei M¨unchen, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Exploring the possibility that fundamental constants of Nature might vary temporallyor spatially constitutes one of the key science drivers for the European Southern Ob-servatory’s ESPRESSO spectrograph on the VLT and for the HIRES spectrographon the ELT. High-resolution spectra of quasar absorption systems permit accuratemeasurements of fundamental constants out to high redshifts. The quality of newdata demands completely objective and reproducible methods. We have developeda new fully automated Artificial Intelligence-based method capable of deriving opti-mal models of even the most complex absorption systems known. The AI structure isbuilt around VPFIT, a well-developed and extensively-tested non-linear least-squarescode. The new method forms a sophisticated parallelised system, eliminating humandecision-making and hence bias. Here we describe the workings of such a system andapply it to synthetic spectra, in doing so establishing methods of importance for fu-ture analyses of VLT and ELT data. The results show that modelling line broadeningfor high-redshift absorption components should include both thermal and turbulentcomponents. Failing to do so means it is easy to derive the wrong model and henceincorrect parameter estimates. We also argue that model non-uniqueness can be sig-nificant, such that it is not feasible to expect to derive an unambiguous estimate ofthe fine structure constant α from one or a small number of measurements. No matterhow optimal the modelling method, it is a fundamental requirement to use a largesample of measurements to meaningfully constrain temporal or spatial α variation. Key words:
Cosmology: cosmological parameters; Methods: data analysis, numer-ical, statistical; Techniques: spectroscopic; Quasars: absorption lines; Line: profiles;Abundances
The development of optimal, objective, and reproduciblemethods for the analysis of high-resolution absorption sys-tems is essential for analysing the high quality data nowbeing obtained on new astronomical facilities. Forthcomingfacilities like the ELT clearly demand sophisticated analytictools. Searches for new physics, revealed by temporal or spa-tial variations in fundamental constants, constitute one ofthe three main science drivers for the ESPRESSO spectro-graph on the European Southern Observatory’s VLT (Pepe (cid:63)
E-mail: [email protected] † E-mail: [email protected] et al. 2014), and one of the key goals for the forthcomingELT (Hook 2009; ESO ELT team 2010, 2011; Liske et al.2014; Marconi et al. 2016).Bainbridge & Webb (2017a,b) first developed AI meth-ods that were capable of modelling arbitrarily complex ab-sorption systems without any human interactions. At theheart of the process is VPFIT (Carswell & Webb 2014,2020), a non-linear least squares system with tied parameterconstraints. The additional genetic code acts as a wrapper,guiding the descent direction to an optimal solution on thebasis of a test statistic (the corrected Akaike InformationCriterion). Hybrid algorithms of this sort have been coined“memetic” (Moscato 1989).In the present paper, we develop these ideas further, © a r X i v : . [ a s t r o - ph . C O ] A ug C. C. Lee et al the key aims being to find faster algorithms and improvethe overall performance, particularly for complex absorptionsystems/datasets that may be challenging to model, com-prising a wide range of line strengths, from heavily saturatedto barely detected, allowing for any blended absorption linesfrom other systems.Performance checking focuses on estimating the finestructure constant in quasar absorption systems, becausethis is a particularly challenging problem, more so than sim-ple 3-parameter absorption line fitting. Whilst the newlycreated methods introduced in this paper have been devel-oped in the context of quasar spectroscopy, they are likelyto be useful in any absorption line spectroscopy application,including stellar spectroscopy and terrestrial experiments.In Section 2 we explain the various stages of the AI algo-rithm, which we call “AI-VPFIT”. We partition the processinto 6 stages, describing the design of each in detail. Section3 evaluates the performance of AI-VPFIT using syntheticspectra. Section 4 shows that it is important to properlymodel the absorption line width and that not doing has asignificant impact on parameter estimation. In Section 5, wediscuss how an objective AI method removes biases inherentto previous interactive analyses and show that model non-uniqueness limits the accuracy achievable from a single mea-surement of variation in the fine structure constant ∆ α / α .Extensive use is made of the corrected Akaike InformationCriterion for comparing the relative improvement in modeldevelopment at each stage. The pros and cons of this arediscussed in Section 6. Finally, in Section 7 we summariseAI-VPFIT’s key characteristics and performance, the mainfindings from extensive testing using synthetic data, and theimplications and requirements of future measurements toconstrain ∆ α / α using quasar spectroscopy. The new algorithms described in this paper follow the broadprinciples outlined in Bainbridge & Webb (2017a,b) al-though there are notable differences and many new algo-rithms. These include: (i) as described in the following text,trial line positions are random; (ii) the ordering of the vari-ous tasks carried out here is different here compared to Bain-bridge & Webb (2017a,b); (iii) Bayesian model averaging isnot used here; we next describe each stage of the modellingprocess; (iv) GVPFIT was set up to use either fully turbulentor fully thermal broadening. This turns out to be a signifi-cant issue (discussed in detail later in this paper). The newmethod presented here does not make these assumptions; (v)two new stages are included to refine the model and addressthe overfitting problem.
In general, the combined dataset will be comprised of multi-ple spectral segments, to be fitted simultaneously. In Stage 1,partially for computing speed considerations, we first build apreliminary model prior to modelling the entire dataset. Thechoice of the primary reference transition (or transitions ) ora primary reference species is important. A simple example might be fitting an absorption sys-tem comprised of multiple FeII transitions and the MgII2796/2803 doublet. Let us assume, in this example, thatall transitions fall longwards of the Lyman- α forest. If thisis a damped Lyman- α system in which both MgII lines aremostly saturated, the sensible primary species choice is FeII.On the other hand, for a non-DLA system, with gen-erally lower column densities, the MgII lines may be strongbut unsaturated and the FeII features may be sufficientlyweak that some velocity components seen in MgII are belowthe detection threshold in FeII. In this case, the MgII ab-sorption complex provides the better constraint on velocitystructure and may be more suited as the primary species.To give further examples, the requirement may be touse a single species, e.g. FeII to solve for the fine structureconstant α , or many H lines to solve for the electron toproton mass ratio µ . In these cases, the primary transition (or transitions) are selected to be the data segment (or seg-ments) that are likely to provide the most reliable initialmodel. To clarify, consider the following example. Supposewe have two FeII transitions, 1608 and 2383˚A, the formerfalling in the Lyman- α forest (and hence being blended withHI lines at various redshifts), the latter falling longwards ofthe Lyman- α emission line. In this situation, the FeII 2383˚Aline should be treated as the primary transition and the FeII1608˚A line as the secondary transition.The fine structure constant α (and/or if fitting molec-ular lines, the electron to proton mass ratio µ ) may be in-cluded as a free parameter in this Stage. However, this canonly be done if the primary species or primary transitionsare comprised of 2 or more transitions with sufficiently dif-ferent sensitivities to α variation so as to avoid degeneracywith redshift z . It is important to include these parametersfrom the outset, rather than first deriving a best-fit modelfor the entire dataset and only then adding those free pa-rameters, since this would inevitably bias α (or µ towardsthe terrestrial value).Having decided on the species to be used as the pri-mary species , Monte Carlo methods are used to constructthe model, simultaneously for all primary species’ transi-tions. Initially, a single model absorption line is placed atrandom within the complex. Its absorption line parameters( N and b ) are default initial values, always the same, user-defined and provided by an input parameter file.The algorithm above thus provides an initial first-guessmodel for the primary species. Non-linear least-squares, VP-FIT, (Carswell & Webb 2014, 2020), is used to refine the ini-tial parameters. At this point (initial iteration), we are likelyto have a very poor fit to the data, using only a single profileto model an arbitrarily complex dataset. The goodness-of-fitis quantified using the corrected Akaike Information Crite-rion (AICc), AICc = χ p + nkn − k − , (1)where n and k are the numbers of data points and modelparameters, respectively, and χ p = n (cid:213) i = ( F i , data − F i , model ) σ i (2)where F i , data , F i , model , and σ i are the observed spectralflux, the model fit, and the estimated error on F i , model for MNRAS000
In general, the combined dataset will be comprised of multi-ple spectral segments, to be fitted simultaneously. In Stage 1,partially for computing speed considerations, we first build apreliminary model prior to modelling the entire dataset. Thechoice of the primary reference transition (or transitions ) ora primary reference species is important. A simple example might be fitting an absorption sys-tem comprised of multiple FeII transitions and the MgII2796/2803 doublet. Let us assume, in this example, thatall transitions fall longwards of the Lyman- α forest. If thisis a damped Lyman- α system in which both MgII lines aremostly saturated, the sensible primary species choice is FeII.On the other hand, for a non-DLA system, with gen-erally lower column densities, the MgII lines may be strongbut unsaturated and the FeII features may be sufficientlyweak that some velocity components seen in MgII are belowthe detection threshold in FeII. In this case, the MgII ab-sorption complex provides the better constraint on velocitystructure and may be more suited as the primary species.To give further examples, the requirement may be touse a single species, e.g. FeII to solve for the fine structureconstant α , or many H lines to solve for the electron toproton mass ratio µ . In these cases, the primary transition (or transitions) are selected to be the data segment (or seg-ments) that are likely to provide the most reliable initialmodel. To clarify, consider the following example. Supposewe have two FeII transitions, 1608 and 2383˚A, the formerfalling in the Lyman- α forest (and hence being blended withHI lines at various redshifts), the latter falling longwards ofthe Lyman- α emission line. In this situation, the FeII 2383˚Aline should be treated as the primary transition and the FeII1608˚A line as the secondary transition.The fine structure constant α (and/or if fitting molec-ular lines, the electron to proton mass ratio µ ) may be in-cluded as a free parameter in this Stage. However, this canonly be done if the primary species or primary transitionsare comprised of 2 or more transitions with sufficiently dif-ferent sensitivities to α variation so as to avoid degeneracywith redshift z . It is important to include these parametersfrom the outset, rather than first deriving a best-fit modelfor the entire dataset and only then adding those free pa-rameters, since this would inevitably bias α (or µ towardsthe terrestrial value).Having decided on the species to be used as the pri-mary species , Monte Carlo methods are used to constructthe model, simultaneously for all primary species’ transi-tions. Initially, a single model absorption line is placed atrandom within the complex. Its absorption line parameters( N and b ) are default initial values, always the same, user-defined and provided by an input parameter file.The algorithm above thus provides an initial first-guessmodel for the primary species. Non-linear least-squares, VP-FIT, (Carswell & Webb 2014, 2020), is used to refine the ini-tial parameters. At this point (initial iteration), we are likelyto have a very poor fit to the data, using only a single profileto model an arbitrarily complex dataset. The goodness-of-fitis quantified using the corrected Akaike Information Crite-rion (AICc), AICc = χ p + nkn − k − , (1)where n and k are the numbers of data points and modelparameters, respectively, and χ p = n (cid:213) i = ( F i , data − F i , model ) σ i (2)where F i , data , F i , model , and σ i are the observed spectralflux, the model fit, and the estimated error on F i , model for MNRAS000 , 1–14 (2020)
I-VPFIT the i th pixel and the subscript p indicates a summation forthe primary data only.As is well known, the second term in Eq. (1) penalises χ , to prevent an arbitrarily large number of free parame-ters being introduced. Other statistics such as the BayesianInformation Criterion (BIC) perform a similar function butimpose a stronger penalty. We defer a comparison betweenthe various options to a forthcoming paper.At the end of this process (Generation 1) We now havethe best-fit single-profile model and can thus begin to refinethe model i.e. increase its complexity until we have a sta-tistically acceptable representation of the data. This best-fit3-parameter model thus now becomes the parent for Gen-eration 2, in which a second trial line is placed randomlyin redshift z within the fitting range , assigning the trialline the same initial default parameters as before. Non-linearleast-squares refinement is again carried out and AICc com-puted. If AICc has decreased compared to Generation 1,the current model is accepted as the Generation 2 best-fit,becoming the parent for the next generation. if AICc has increased compared to Generation 1, the model is rejected,the process repeated (i.e. a second trial line is used, againrandomly placed in redshift), until AICc decreases. We calleach iteration of this procedure one generation.Iterating the above gradually increases the model com-plexity and improves the fit. The loop is terminated onlyafter no reduction in AICc can be found after a user-definednumber of trials (parameter N line , set in the same initialisa-tion file). Stage 1 results in a good fit to the primary species alone(MgII in our example). Stage 2 comprises two steps. Firstly,the overall dataset being fitted is increased by adding inall secondary species (FeII in our example). The absorptionsystem structure derived in Stage 1, i.e. the best-fit set ofredshifts, is replicated for each secondary species to be fitted.The initial model for the secondary species is relatively crudein that each redshift component for each of the secondaryspecies is assigned default trial values of the column density N . The velocity dispersion parameters b can be related eitherby the limiting cases of entirely “thermal” broadening ( b s = b p (cid:112) m p / m s ), where the subscripts p and s refer to primaryand secondary species), entirely “turbulent” broadening ( b s = b p ), or in-between (i.e. the gas temperature T is included asa free parameter). This is discussed in detail in Section 2.7).Then, non-linear least-squares refines the parameters for allsecondary species, after which AICc is again computed andstored.If solving for either α or µ , and if this free parameter wasnot introduced in Stage 1 (see earlier discussion), it can beincluded at this Stage. Including α (or µ ) prior to developinga complete model is important to avoid bias.The second step within Stage 2 is to increase the modelcomplexity, i.e. further redshift components are added, oneat a time, until no descent in AICc can be found for N line trials. Once this happens, the Stage 2 model is complete. In practice, the range is defined by the lowest and highest red-shift available over the multiple data segments used.
So far, we have assumed that no unidentified absorption linesare present anywhere in the data. This is generally not thecase. Each quasar sightline typically intersects multiple dis-tinct absorption systems, such that blending between linesarising in different redshift systems is common.There are two types of interlopers: (i) an absorptionsystem at some other known redshift. In this case the in-terloper species and rest-wavelength may be identified andthere may be other transitions in the same ion or from otherions at the same redshift that can be modelled simultane-ously in order to best constrain the interloper parameters,or (ii) the interloper species/origin is unknown. If the inter-loper is identified beforehand as being of type (i), then theoverall model to be fitted can include the appropriate freeparameters from the outset, that is, the set-up prior to Stage1 indicates that 2 redshift systems are to be simultaneouslyfitted.If we have an interloper of type (ii), all previous (non-interloper) line parameters from Stages 1 and 2 are initially(here, in Stage 3) fixed . The reason for doing this is as fol-lows. χ minimisation will act so as to ‘prefer’ the interloper(since the interloper has no tied parameters, unlike the heavyelement line) such that the interloper may tend to replacethe heavy element line. Physically, this is obviously incor-rect.To identify possible interlopers, there are 2 potentialprocedures: (a) Each data segment is treated separately. Aninterloper is added in a random position within each datasegment, one at a time, and VPFIT allowed to iterate toderive the best-fit interloper parameters (or to reject theinterloper entirely), or (b) All data segments are dealt withat the same time i.e. one interloper is placed in each datasegment simultaneously, with wavelength λ iinterloper = ( − R ) × λ imin + R × λ imax (3)where R is the random seed from − and λ imin ( max ) is theminimum (maximum) wavelength of the i th segment. Onecould generate a different random seed for each wavelengthsegment but this offers no advantage.Superficially, (a) is the favoured option because (b) hasa potential systematic trend towards over-fitting the data(i.e. fitting the data with too many absorption components).The reason for this is because introducing one single in-terloper at a time, and using AICc to check whether thatinterloper is justified, results in a final model that only con-tains parameters that are statistically required. On the otherhand, introducing one interloper into each segment simulta-neously does not permit individual AICc testing on eachinterloper (because of tied parameters, it is not possible todefine AICc on a region-by-region basis). This argument sug-gests that procedure (a) is preferable in order to avoid over-fitting. However, procedure (a) is computationally very timeconsuming – e.g. if we have 10 data segments to be fitted,the time required is 10 times as long. In terms of computingtime, it is more efficient to adopt procedure (b) and sub-sequently test for over-fitting and repair the model wherenecessary (Stage 5). MNRAS , 1–14 (2020)
C. C. Lee et al
This Stage comprises two steps.
Step 1:
Firstly, we release all the previously-fixed primaryand secondary parameters such that they can be freely var-ied. Then, additional primary and secondary components areadded, one per generation. No new interlopers are included.As before, this process terminates only when no AICc de-scent occurs within N line trials. Step 2:
At this point we have only considered parametersassociated with absorption lines themselves. We have notyet considered additional important parameters, such as thecontinuum level nor a possible zero level correction. Thesecond step in Stage 4 is therefore to add in these additionalparameters (if they are required), and refine the model againusing AICc and N line . This Stage allows the model to evolve further to allow forthe introduction of new continuum and zero-level parame-ters. By the end of Stage 4, a good model has already beenachieved. However, since that model has been derived priorto including continuum and zero-level parameters, severalprocedures need repeating. Therefore:(i) All continuum and zero-level parameters are temporar-ily fixed;(ii) The algorithm is returned to Stage 2, but the Stage 2starting (i.e. first guess) model parameters are those fromthe end of Stage 4 (step 2), but all interlopers previouslyfound are removed;(iii) Stages 2 to 4 are then repeated. The previously-fixedcontinuum and zero-level parameters are returned to freeparameters are the (repeated) Stage 4/step 2. The loopis closed by a user-set (or default) AI-VPFIT setting,which defines the number of repeats.We now explain the reasons behind the 3 processes im-mediately. If, prior to this stage, the original continuumplacement was slightly wrong, then because we may be fit-ting multiple transitions from a single species (e.g. severalFeII transitions), a continuum error may result in a poor fitonce secondary transitions are included.This point can be clarified by example. Suppose ourprimary species is FeII2383 and one secondary transitionis the weaker (lower oscillator strength) FeII2344. Supposefurther that the original continuum placement for FeII2344was placed slightly too low. In this example, it could easilybe the case that an interloper (coincident with the FeII2344position) provides a greater reduction in AICc comparedto the (real) FeII2344 line. Such an effect clearly producesproblems as follows. The lowest AICc may then be found byreducing the FeII column density so as to achieve a good fitat both lines, but an interloper is placed in the FeII2383 lineto compensate for the continuum placement error.Slight errors in the original continuum placement canalso lead to the following problem. If, again, the originalcontinuum had been slightly too low in places, interloperswhich in reality are present in the spectrum may have beenmissed. Repeating Stages 2 and 3 guards against this possi-bility. Reproduce velocity structure from primary species intoall secondary species. Reason: some of the secondary speciescomponents could have been dropped. Therefore put allsecondary components back but remove all dropping cri-teria. Setting dropping criteria requires striking a balancebetween avoiding insignificant parameters in the final modeland hence having ill-conditioned matrices yet at the sametime avoiding losing relevant parameters. One reason for do-ing this is that in a system with multiple velocity compo-nents, during the iterative sequence, a parameter which isultimately determined to be important, temporarily falls be-low a dropping criterion, so could be incorrectly removed. Byremoving the a dropping criterion, the search direction (andstep) can become unstable through ill-conditioning. Thiswould cause difficulties in the absence of a temporary con-trol. Therefore, during this Stage, we place an upper boundon parameter update step sizes. The second action taken at this stage is to fix any contin-uum and zero level parameters to those best-fit values ob-tained at the end of Stage 4. This turns out to be helpful forthe following reason. Since here, in Stage 5, we re-initialisethe secondary species absorption line parameters (either to aset of default values or otherwise), the model is temporarilyplaced further away from the best-fit solution. The potentialconsequence of this can be that any sudden increase in (forexample) a continuum parameter level can be compensatedby a corresponding increase in absorption line strength (i.e.an increase in either N or b or both). In other words, a sat-isfactory model may sometimes be found corresponding toa local minimum in χ space. This problem is avoided bypreventing (temporarily) further adjustment in either thecontinuum or zero level parameters. The third action taken at this stage is to temporarilyfix all redshifts (for both heavy element lines and interlop-ers). Again, because this Stage temporarily moves the cur-rent model away from the final best-fit model, as in point2 above, the slight decrease in overall stability increases thechance of a column density being momentarily sent to anartificially low value. When this happens, because that linemay become very weak, its Hessian components correspond-ing to its redshift also become small, and its parameter up-date can become large. This means it can (mistakenly) bemoved far from its current (appropriate) position to some ar-bitrary (inappropriate) position. Since this is hard to recoverfrom, an easy solution is to temporarily fix all absorptionline redshifts, whilst the secondary species column densities‘recover’ from default-initialisation values to near-optimalvalues. After minimising χ to refine all free parameters above,the fourth action taken at this Stage is to remove the tempo-rary parameter holds described above and again to minimise χ to refine all interesting parameters. Interlopers are introduced to the model in Stage 3. Ifthe spectral dataset comprises multiple segments (i.e. sev-eral transitions from several different atomic species), one
MNRAS000
MNRAS000 , 1–14 (2020)
I-VPFIT trial interloper is added to every segment simultaneously,the parameters are refined, and interlopers are kept if the overall AICc decreases; in fact, an interloper might increasethe local
AICc (i.e. for the spectral segment it lies in). Forcomputational efficiency, to ensure that spurious interlopersdo not remain in the model, this is dealt with in Stage 5 byremoving each interloper one at a time and only retainingeach if the overall AICc does not increase. It is possible that at this Stage of the whole process thatthe model contains regions where too many heavy elementparameters have been introduced, i.e. over-fitting has oc-curred. There are several ways in which over-fitting couldoccur for the heavy element lines: (i)
In Stage 1, the primarymodel is developed. In Stage 2, the velocity structure fromthe primary species is replicated to ensure that any velocitycomponents present in the primary species are also presentin the secondary species. However, local least-squares fittingin VPFIT could result in one of the secondary species’ com-ponents falling below a threshold and being rejected from themodel. This problem is remedied in the ‘Preparation’ stepabove, in order to keep the model ‘physical’. Thus, the re-placement of secondary species components is based purelyon physical grounds and not on statistical grounds, i.e. infact the replacement could produce a higher value of AICc,but this is not considered. For this reason, the replacementjust described (i.e. requiring all velocity components to bepresent in all species), creates a marginal tendency to overfitthe data. (ii)
Another potential way in which slight overfit-ting can be created is as follows. In Stage 1, when the pri-mary species model is developed, velocity structure is builtup by adding one line at a time. This is done using a MonteCarlo process (i.e. trial lines are placed randomly within thedata fitting region, repeating the procedure several timesand selecting the best (lowest AICc) solution from the set oftrials. Because the number of trials is finite, it is occasionallypossible to find a local rather than an absolute AICc min-imum when accepting a velocity component. In this sense,the order in which velocity components are added into themodel can influence the final best-fit primary species model.This provides another mechanism by which spurious com-ponents may remain in the best-fit primary species model,which are carried forward to the next Stages.
If the incorrect broadening model is used to generate theabsorption profiles (when fitting multiple atomic species si-multaneously), this can impact adversely on the derived ve-locity structure. Consider a simple absorption system com-prising only MgII2796 and FeII2383 profiles. Suppose theintrinsic gas broadening is thermal but that we use turbu-lent broadening to model the observed profiles. In this case,the model MgII profile will be slightly too narrow and themodel FeII profile will be slightly too broad. Two mecha-nisms can easily compensate for this: (a) a model with 2velocity components can match the observed profiles by ad-justing the relative column densities accordingly, or (b) aninterloper can be blended with the MgII profile and b turb adjusted accordingly to match the data. When α is an addi-tional fitting parameter, the effects just described can add asignificant systematic uncertainty.The discussion above illustrates the importance of using the correct broadening method. The procedure to achievethis is:(i) Where applicable (i.e. where different atomic specieswith sufficiently different atomic masses are present) the b -parameters are changed such that they are related(tied) by b = b turb + b th , where b th = kT / m , m is atomicmass, and T is the gas temperature. Clearly this requiresa new free parameter, T , for each velocity component inthe absorption complex.(ii) The temperature parameter T can be included, if ap-propriate, at the beginning of Stage 2. A fully automated modelling process such as the one de-scribed in this paper opens new opportunities to scrutinisethe limitations present in previous analyses. In particular, wedo not know a priori whether the broadening of any partic-ular absorption component is purely thermal, purely turbu-lent, or a hybrid. Many earlier analyses of quasar absorptionsystems assumed the broadening mechanism from the outset(usually assumed to be turbulent) and constructed modelsbased on that assumption. The justification for doing so wasthat the assumption should not bias the measured ∆ α / α oneway or the other. Whilst this is true, an important point hasbeen overlooked; fitting the wrong model inevitably not onlyleads to incorrect velocity structure, but importantly, pro-duces systematically over-fitted models. This has importantand very undesirable consequences, as we now demonstrateusing synthetic spectra. We base synthetic spectra on the z abs = . absorptioncomplex in the spectrum of the bright quasar HE 0515-0414.The choice is unimportant and some other system couldequally have been used. This particular absorption complexspans an unusually large redshift range, . < z abs < . , corresponding to a velocity range ∼ km/s. Itis not necessary to use such an extensive range for the pur-poses required here. We therefore use a subset of the systemspanning the redshift range . < z abs < . , corre-sponding to a velocity range ∼ km/s. For simplicity wesimulate only three transitions, Mg II 2796 and 2803 ˚A andFe II 2383 ˚A.The simulated spectra are Voigt profiles convolved witha Gaussian instrumental profile with σ res = . km/s. Thepixel size is approximately . km/s and the signal to noiseratio per pixel is 100 (the noise is taken to be Gaussian). Theresolution and pixel size correspond to those of the HARPSinstrument on the ESO 3.6m telescope. We derive the in-put absorption line parameters as follows. In the redshiftrange . < z abs < . , using AI-VPFIT to analyseexisting HARPS spectra, we find a total of eight heavy ele-ment absorption components plus one interloper. We modelthe (real) data including a temperature parameter for eachabsorption component in the complex i.e. each componenthas both a turbulent and a thermal component to its ob-served b -parameter. This AI-VPFIT model becomes one of MNRAS , 1–14 (2020)
C. C. Lee et al the synthetic models. We call this the “temperature” model.The detailed parameters used to create the temperature syn-thetic spectrum are given in Table 1.
Fig.1 illustrates the evolution of the model at various pointsduring Stage 1. Each panel shows a different generation. Theinitial (Scratch) model is a single randomly-placed line. Atthe end of Stage 1, the line has broadened and the columndensity increased to minimise χ . The subsequent genera-tions show how the model gradually evolves as more compo-nents are added. Only the primary species is shown (MgIIin this case), as secondary species are not yet included. Bygeneration 10, the model is already good, although no inter-lopers have been considered so the feature at ∼ − km/s is(incorrectly, but knowingly) modelled as MgII. This prob-lem is caught at Stage 3. In the example shown, AI-VPFITcontinued iterating up to generation 29 but (on the basisof AICc), decided no additional velocity components wererequired, i.e. the model did not change for the last 19 gen-erations.Fig.2 shows the development of the model through eachStage. For Stages 1 through 5, the model illustrated is fromthe end of those stages. In the example used, generation10 is the end-point of Stage 1 so is replicated as the top-leftpanel in Fig.2. The model at the end of Stage 1 comprises 10components. In Stage 2, the secondary species (in this caseonly FeII) is included. The panel illustrating Stage 6 showsthe initial condition of that Stage, where all primary species’velocity components present by the end of Stage 5 are in-cluded but all corresponding secondary species componentsare replaced . In the example used, 3 weak FeII components(at about − , + , and + km/s) are replaced as the ini-tial condition for Stage 6. The motivation for doing so is tocheck that a genuine heavy element velocity component inthe secondary species has not been mistakenly replaced byan interloper in Stage 5 (this is possible since the decisionmaking is based purely on AICc).The final best-fit velocity structure is virtually indistin-guishable from true model, labelled (“Fiducial”) and high-lighted with a red box to indicate this is the input and not afitted model: the number of derived absorption componentsis 8 with one interloper in the MgII2796 line. The decision as to which type of absorption line broad-ening to use when fitting real data is important. As alludedto earlier, previous analyses have generally assumed turbu-lent broadening on the basis that the assumption cannotbias estimating ∆ α / α . This assumption seems reasonable in-tuitively and is likely to be correct in a statistical sense.However, in fact it turns out that the choice of line broad-ening in the modelling procedure is very important.In Section 3, AI-VPFIT was applied to a syntheticspectrum generated with temperature broadening. This is of course not the only option. Most modelling of high res-olution quasar spectra has previously been carried out as-suming either turbulent or thermal broadening or both. Herewe show the pitfalls of these assumptions and argue that, forreliable results, it is important to include an additional freeparameter, T, for each velocity component in the model.We begin by generating 2 additional synthetic spectra;to create the turbulent synthetic model, all parameters arekept the same except for the b -parameter of each absorptioncomponent, which is now set to be equal to the Mg II tur-bulent component. To create the thermal synthetic model,we again use the observed MgII b -parameters, although nowthe FeII b -parameters are adjusted according to their atomicmass. The result of the above is the creation of 3 syntheticspectra, one where the line broadening takes into accountboth turbulent and thermal components of the b -parameter,one where the broadening is assumed to be purely turbulent,and one where the broadening is assumed to be purely ther-mal. The detailed parameters used to create the turbulent and thermal synthetic spectra are given in Table 1. Thosemodel spectra (and the corresponding temperature model),are as shown in Fig.3.Fig. 4 shows a set of nine results. Each of the threesynthetic spectra (turbulent, thermal, and temperature linebroadening) was modelled assuming each line broadeningmechanism. The spectrum and model in each panel is asdescribed in the figure caption. The continuous orange ver-tical lines and blue dashed lines indicate the line centres ofthe heavy element and interlopers absorption components.Figs. 5, 6, and 7 present the central values and errors for the log N , z and b -parameters from each best-fit model in Fig.4. These figures illustrate that modelling a turbulent spec-trum with a thermal model leads to errors, and vice versa .AI-VPFIT sometimes needs two components to achieve thebest-fit model for what is in reality a single component. Theparameter error estimates are typically poorly determinedfor these “double components”.Whilst ensuring that the model solves simultaneouslyfor both thermal and turbulent components of the observed b -parameter, this does come with a cost, and may not al-ways be possible. If the species fitted simultaneously (i.e. theprimary and secondary species) have fairly similar atomicmasses, b turb and b th approach degeneracy, leading to ahuge uncertainty estimates on the gas temperature. In thesynthetic models explored here, the only two elements “ob-served” are Mg (A=24) and Fe (A=56). Even then, slightdegeneracy translates into quite large uncertainties, as is il-lustrated in the bottom right panel in Figs. 5, 6, and 7.Table 2 summarises the results discussed above. Whenthe synthetic spectrum broadening and fitted model broad-ening match, we can see that the number of absorption com-ponents needed to satisfactorily fit the data is minimised(top black line in each section in the table). Conversely,when the “wrong” model is used, additional components areneeded. It is now well-established that the line broadeningmechanism seen in quasar absorption systems is generally(perhaps never) entirely turbulent or entirely thermal i.e. itis important to use a temperature model. Importantly, for allthree synthetic spectra, a temperature model performs verywell in terms of finding the correct number of components(last column), as would be expected.Fig.8 illustrates the results obtained from modelling the MNRAS000
Fig.1 illustrates the evolution of the model at various pointsduring Stage 1. Each panel shows a different generation. Theinitial (Scratch) model is a single randomly-placed line. Atthe end of Stage 1, the line has broadened and the columndensity increased to minimise χ . The subsequent genera-tions show how the model gradually evolves as more compo-nents are added. Only the primary species is shown (MgIIin this case), as secondary species are not yet included. Bygeneration 10, the model is already good, although no inter-lopers have been considered so the feature at ∼ − km/s is(incorrectly, but knowingly) modelled as MgII. This prob-lem is caught at Stage 3. In the example shown, AI-VPFITcontinued iterating up to generation 29 but (on the basisof AICc), decided no additional velocity components wererequired, i.e. the model did not change for the last 19 gen-erations.Fig.2 shows the development of the model through eachStage. For Stages 1 through 5, the model illustrated is fromthe end of those stages. In the example used, generation10 is the end-point of Stage 1 so is replicated as the top-leftpanel in Fig.2. The model at the end of Stage 1 comprises 10components. In Stage 2, the secondary species (in this caseonly FeII) is included. The panel illustrating Stage 6 showsthe initial condition of that Stage, where all primary species’velocity components present by the end of Stage 5 are in-cluded but all corresponding secondary species componentsare replaced . In the example used, 3 weak FeII components(at about − , + , and + km/s) are replaced as the ini-tial condition for Stage 6. The motivation for doing so is tocheck that a genuine heavy element velocity component inthe secondary species has not been mistakenly replaced byan interloper in Stage 5 (this is possible since the decisionmaking is based purely on AICc).The final best-fit velocity structure is virtually indistin-guishable from true model, labelled (“Fiducial”) and high-lighted with a red box to indicate this is the input and not afitted model: the number of derived absorption componentsis 8 with one interloper in the MgII2796 line. The decision as to which type of absorption line broad-ening to use when fitting real data is important. As alludedto earlier, previous analyses have generally assumed turbu-lent broadening on the basis that the assumption cannotbias estimating ∆ α / α . This assumption seems reasonable in-tuitively and is likely to be correct in a statistical sense.However, in fact it turns out that the choice of line broad-ening in the modelling procedure is very important.In Section 3, AI-VPFIT was applied to a syntheticspectrum generated with temperature broadening. This is of course not the only option. Most modelling of high res-olution quasar spectra has previously been carried out as-suming either turbulent or thermal broadening or both. Herewe show the pitfalls of these assumptions and argue that, forreliable results, it is important to include an additional freeparameter, T, for each velocity component in the model.We begin by generating 2 additional synthetic spectra;to create the turbulent synthetic model, all parameters arekept the same except for the b -parameter of each absorptioncomponent, which is now set to be equal to the Mg II tur-bulent component. To create the thermal synthetic model,we again use the observed MgII b -parameters, although nowthe FeII b -parameters are adjusted according to their atomicmass. The result of the above is the creation of 3 syntheticspectra, one where the line broadening takes into accountboth turbulent and thermal components of the b -parameter,one where the broadening is assumed to be purely turbulent,and one where the broadening is assumed to be purely ther-mal. The detailed parameters used to create the turbulent and thermal synthetic spectra are given in Table 1. Thosemodel spectra (and the corresponding temperature model),are as shown in Fig.3.Fig. 4 shows a set of nine results. Each of the threesynthetic spectra (turbulent, thermal, and temperature linebroadening) was modelled assuming each line broadeningmechanism. The spectrum and model in each panel is asdescribed in the figure caption. The continuous orange ver-tical lines and blue dashed lines indicate the line centres ofthe heavy element and interlopers absorption components.Figs. 5, 6, and 7 present the central values and errors for the log N , z and b -parameters from each best-fit model in Fig.4. These figures illustrate that modelling a turbulent spec-trum with a thermal model leads to errors, and vice versa .AI-VPFIT sometimes needs two components to achieve thebest-fit model for what is in reality a single component. Theparameter error estimates are typically poorly determinedfor these “double components”.Whilst ensuring that the model solves simultaneouslyfor both thermal and turbulent components of the observed b -parameter, this does come with a cost, and may not al-ways be possible. If the species fitted simultaneously (i.e. theprimary and secondary species) have fairly similar atomicmasses, b turb and b th approach degeneracy, leading to ahuge uncertainty estimates on the gas temperature. In thesynthetic models explored here, the only two elements “ob-served” are Mg (A=24) and Fe (A=56). Even then, slightdegeneracy translates into quite large uncertainties, as is il-lustrated in the bottom right panel in Figs. 5, 6, and 7.Table 2 summarises the results discussed above. Whenthe synthetic spectrum broadening and fitted model broad-ening match, we can see that the number of absorption com-ponents needed to satisfactorily fit the data is minimised(top black line in each section in the table). Conversely,when the “wrong” model is used, additional components areneeded. It is now well-established that the line broadeningmechanism seen in quasar absorption systems is generally(perhaps never) entirely turbulent or entirely thermal i.e. itis important to use a temperature model. Importantly, for allthree synthetic spectra, a temperature model performs verywell in terms of finding the correct number of components(last column), as would be expected.Fig.8 illustrates the results obtained from modelling the MNRAS000 , 1–14 (2020)
I-VPFIT ScratchGeneration 1Generation 2Generation 3Generation 5 Generation 10Generation 9Generation 8Generation 7Generation 6
Figure 1.
Model evolution through Stage 1 of the AI-VPFIT process. The black histogram illustrates the synthetic temperature spectrum.The signal-to-noise per pixel is 100. The continuous red line in the panel labelled “Scratch” shows the initial condition i.e. the randomlyplaced first-guess line. No VPFIT generation has taken place. The continuous red line in subsequent panels (i.e. “Generation ”) showsthe best-fit temperature model after each generation of VPFIT. The ± σ normalised residuals are plotted at the top of each panel.MNRAS , 1–14 (2020) C. C. Lee et al
Stage 1 (End)Stage 2 (End) Stage 5 (End)Stage 3 (End) Stage 6 (Initialised)Final ModelFiducialStage 4 (End)
Figure 2.
Model evolution at each Stage in AI-VPFIT. The black histogram shows the synthetic spectrum (temperature line broadening).The signal-to-noise per pixel is 100. The red continuous line illustrates the best-fit temperature model. For Stages 1 through 5, the end-model is shown. For Stage 6, we illustrate the model at the commencement of the stage, to show how the current velocity structure ofthe primary species (Mg in this case) is re-introduced to the secondary species (FeII), for the reasons explained in Section 2.6.MNRAS000
Model evolution at each Stage in AI-VPFIT. The black histogram shows the synthetic spectrum (temperature line broadening).The signal-to-noise per pixel is 100. The red continuous line illustrates the best-fit temperature model. For Stages 1 through 5, the end-model is shown. For Stage 6, we illustrate the model at the commencement of the stage, to show how the current velocity structure ofthe primary species (Mg in this case) is re-introduced to the secondary species (FeII), for the reasons explained in Section 2.6.MNRAS000 , 1–14 (2020)
I-VPFIT Temperature Turbulent Thermal
No. log N ( MgII ) log N ( FeII ) Redshift b turb T /K b th ( Mg ) b th ( Fe ) b ( Mg ) = b ( Fe ) b ( Mg ) b ( Fe ) Interloper log N Redshift b
Table 1.
Model parameters used to generate the three synthetic spectra discussed in Sections 3 and 4. The model comprises eight heavyelement components and one interloper. Columns 5-8 correspond to the temperature synthetic spectrum, column 9 to the turbulent one,and columns 10, 11 to the thermal one. Columns 7 and 8 are the thermal b -parameters for MgII and FeII, i.e (cid:112) kT / m ( Mg, Fe ) . MgIIand FeII have the same b -parameter for turbulent broadening, as shown in column 9. Columns 10 and 11 give the thermal b -parametersfor MgII and FeII with b ( FeII ) = b ( MgII ) (cid:112) m ( Mg )/ m ( Fe ) . All b -parameters have units of km/s. -60 -40 -20 0 20 4000.51 MgII 2796.36
FeII 2382.76
Figure 3.
Upper panel: FeII 2383 for turbulent (red continuouscurve), thermal (blue dotted curve) and temperature (thin blackcurve) broadening. Lower panel: MgII 2786, where the line broad-ening is b th (Mg) from Table 1. The figure illustrates the varia-tions between each broadening model and hence highlights theimportance of using the correct fitting model. For example, line2 has a low turbulent b -parameter and a relatively high tempera-ture, so the thermal and temperature line centres are similar andmuch deeper than the turbulent line. Line 8 has a larger turbu-lent b -parameter and a lower temperature, so the turbulent andtemperature models are closer, the deeper line centre being thethermal model. synthetic spectra described above. The figure is divided intothree panels, corresponding (top to bottom) to the syntheticdata being turbulent , thermal , and temperature respectively.Four ∆ α / α measurements are illustrated within each panel.The highest point illustrates the fitted value of ∆ α / α usingVPFIT only, where the starting parameter guesses were theinput parameters i.e. the parameters used to generate thesynthetic spectrum. The second point down illustrates theresult of AI-VPFIT fitting the synthetic data using a turbu-lent model. The third point down illustrates a thermal modeland the lowest point a temperature model. The numerical re- sults are listed in Table 2. The results can be summarisedas follows:(i) When fitting the synthetic data with the correct model,the results are well-behaved. As expected, the fitted ∆ α / α agrees well with the true value (illustrated by thevertical dashed line in Fig.8).(ii) When the wrong model is used, the ∆ α / α uncertainty es-timate is magnified by a factor of 3 or 4. This is becauseimposing the wrong broadening mechanism necessarilyresults in additional velocity components (heavy ele-ments or interlopers or both) being introduced to achievea satisfactory fit to the data. The consequence of theadditional fitting parameters is that ∆ α / α is more por-rly constrained. In the last case ( temperature syntheticdata), the thermal model produces a ∆ α / α estimate thatis almost 2 σ away from the correct value.(iii) For all three synthetic spectra, when a temperature model is used, the results are good, as would be expected,since the model incorporates the full range of possibili-ties.The important conclusion of the above that using thewrong broadening mechanism has highly undesirable conse-quences. The solution is that wherever possible (i.e. whendifferent species with sufficiently different atomic massesare modelled simultaneously), one should use a temperature model. ∆ α / α It is important to understand that the modelling processoutlined in this paper, in the context of solving for ∆ α / α ,is not the same as the interactive process followed by a hu-man being. In previous studies, where no AI methodologyhas been used, and where a human has interactively con-structed an absorption line model, this is done by assumingthroughout that ∆ α / α = . Once a good model has beenfound, the procedure generally adopted in previous analyseshas been to only then include ∆ α / α as a free parameter. MNRAS , 1–14 (2020) C. C. Lee et al FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1 FeII 2382.761 +1 MgII 2796.361 +1
50 0 50
Velocity relative to z abs =1.147179 (km/s) MgII 2803.531 +1
Figure 4.
Synthetic spectra and lowest AICc AI-VPFIT models. Tick marks indicate component positions. The signal to noise per pixelis 100. Other parameters used to generate the synthetic spectra are given in Section 3. The left hand column is for the turbulent syntheticspectrum, modelled using turbulent broadening (top), thermal broadening (middle), and temperature broadening (bottom). The middlecolumn illustrates the thermal synthetic spectrum. The right hand column illustrates the temperature synthetic spectrum. Normalisedresiduals are plotted above each spectrum. The horizontal parallel lines indicate the 1 σ ranges.Spectrum Model ∆ α / α σ ( ∆ α / α ) Turbulent Turb (vp)
Turb
Therm
Temp
Thermal Therm (vp)
Turb -2.86E-07 1.43E-05 11 + 3
Therm
Temp
Temperature Temp (vp)
Turb
Therm
Temp
Table 2.
The fitted ∆ α / α for each permutation of synthetic spec-trum and model type - illustrated graphically in Fig.8. See thecaption to that figure for further details. The last column showsthe numbers of heavy element velocity components plus the num-ber of interlopers in the model, within the velocity range of ± km/s, as illustrated in Fig.4. We show in this Section that in fact very good modelscan be (and are) found with disparate values of ∆ α / α . Givenenough computations, different starting points (equivalentto using different random seeds in the AI-VPFIT code) mayreach different solutions. Whilst this is a well-known char-acteristic of non-linear least-squares procedures applied tocomplex datasets, it has not been studied in any detail inthe context of modelling high-resolution quasar absorptionspectra.Fig.9 shows an AI-VPFIT model, obtained using the thermal synthetic spectrum fitted using a turbulent model.The model illustrated gives good normalised residuals andappears to be a good representation of the data. However,the synthetic spectrum was generated using ∆ α / α = × − ,yet the best-fit value found for this particular turbulent fitis ∆ α / α = ± × − . The model is thus 4 σ away fromthe “true” value in this case. Inspecting the model detailsreveals how this occurs. Referring to Fig.9, the FeII columndensity of velocity component AI has decreased relative to MNRAS000
The fitted ∆ α / α for each permutation of synthetic spec-trum and model type - illustrated graphically in Fig.8. See thecaption to that figure for further details. The last column showsthe numbers of heavy element velocity components plus the num-ber of interlopers in the model, within the velocity range of ± km/s, as illustrated in Fig.4. We show in this Section that in fact very good modelscan be (and are) found with disparate values of ∆ α / α . Givenenough computations, different starting points (equivalentto using different random seeds in the AI-VPFIT code) mayreach different solutions. Whilst this is a well-known char-acteristic of non-linear least-squares procedures applied tocomplex datasets, it has not been studied in any detail inthe context of modelling high-resolution quasar absorptionspectra.Fig.9 shows an AI-VPFIT model, obtained using the thermal synthetic spectrum fitted using a turbulent model.The model illustrated gives good normalised residuals andappears to be a good representation of the data. However,the synthetic spectrum was generated using ∆ α / α = × − ,yet the best-fit value found for this particular turbulent fitis ∆ α / α = ± × − . The model is thus 4 σ away fromthe “true” value in this case. Inspecting the model detailsreveals how this occurs. Referring to Fig.9, the FeII columndensity of velocity component AI has decreased relative to MNRAS000 , 1–14 (2020)
I-VPFIT -40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40051015 Figure 5.
Results of modelling the synthetic turbulent spectrum. Top row: the fitted model is a turbulent model. Middle row: the fittedmodel is thermal . Bottom row: the fitted model is temperature . Blue points (circles) show the best-fit parameters for MgII. Red points(triangles) show the best-fit parameters for FeII. VPFIT error bars for each parameter are shown in all cases. Left: log N versus velocity.Centre: Redshift uncertainty (from VPFIT) versus velocity. Right: b -parameter versus velocity. -40 -20 0 20 408101214 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40051015 Figure 6.
Same as Fig.5 except using the thermal synthetic spectrum.MNRAS , 1–14 (2020) C. C. Lee et al -40 -20 0 20 408101214 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 408101214 -40 -20 0 20 40-2-1012 -40 -20 0 20 40024681012-40 -20 0 20 401011121314 -40 -20 0 20 40-2-1012 -40 -20 0 20 40051015
Figure 7.
Same as Fig.5 except using the temperature synthetic spectrum. -2 -1 0 1 2 310 -5 TurbulentThermalTemperature
TurbTurb (vp)ThermTempTurbTherm (vp)ThermTempTurbTemp (vp)ThermTemp
Synthetic Spectrum
Model Best-fit
Figure 8.
Comparing the fitted ∆ α / α for each permutation ofsynthetic spectrum and model type. The left-hand column indi-cates the synthetic spectrum used. The second column indicatesthe line broadening mechanism used in modelling. “Turb (vp)”means that VPFIT was used with the true model parameters asstarting guesses; this point is shown in black in the third column,which also shows the best-fit ∆ α / α for a turbulent , thermal , and temperature model, with corresponding parameter 1 σ uncertain-ties (from VPFIT). its “true” value, to permit the interloper 24. The same effectoccurs at velocity component AE. The interloper 25 is actu-ally present in the data. The strong line under the AM andAA components is actually a single component. However,AI-VPFIT has introduced two components here, such that
01 <> 0.99663
FeII 2382.76 f = 0.2936 = 171.6400npix = 210 = 1.5463ndf = 111 q = 1505
AW AS AMAA AU AQ AO AE AC AI
21 22 23 24
1 +1
01 <> 1.01232
MgII 2796.36 f = 0.4862 = 175.3800npix = 218 = 1.3388ndf = 131 q = 212 aw as amaa au aq ao ae ac ai
1 +1
50 0 50
Velocity relative to z abs = 1.147179 (km/s)
01 <> 0.99670
MgII 2803.53 f = 0.2416 = 191.3800npix = 218 = 1.4176ndf = 135 q = 121 aw as amaa au aq ao ae ac ai
1 +1 /Users/leechungchi/Documents/VPFIT/velplot/vpfit/data/synthetic/tbtb.13
Figure 9.
Illustration of the model non-uniqueness problem. Theinput spectrum is thermal and the fitted model is turbulent . Thismodel passed all acceptance tests described in Section 2 and is aminimum AICc result with χ n = . . However, it gives ∆ α / α = . ± . × − , compared to the “true” (synthetic spectrum)input value of ∆ α / α = . × − . The discrepancy is caused by alocal χ minimum. the column densities of the FeII AA and MgII AM compo-nent dominate the line centroid positions. The “true” modelin fact comprises 8 heavy element velocity components andone interloper (25). However, Fig.9 comprises 10 heavy ele-ment components and 5 interlopers.In this example, clearly a false minimum in χ has beenfound. There is no a priori reason to reject this model. The MNRAS000
Illustration of the model non-uniqueness problem. Theinput spectrum is thermal and the fitted model is turbulent . Thismodel passed all acceptance tests described in Section 2 and is aminimum AICc result with χ n = . . However, it gives ∆ α / α = . ± . × − , compared to the “true” (synthetic spectrum)input value of ∆ α / α = . × − . The discrepancy is caused by alocal χ minimum. the column densities of the FeII AA and MgII AM compo-nent dominate the line centroid positions. The “true” modelin fact comprises 8 heavy element velocity components andone interloper (25). However, Fig.9 comprises 10 heavy ele-ment components and 5 interlopers.In this example, clearly a false minimum in χ has beenfound. There is no a priori reason to reject this model. The MNRAS000 , 1–14 (2020)
I-VPFIT only potentially tell-tale signs of there being a problem withthis model are (i) an apparent excess of close blends (bothheavy element-heavy element and heavy element-interloper)in the stronger components, and (ii) the value of ∆ α / α ob-tained is inconsistent with that obtained using the turbulent and thermal models Fig.4.The effect described above, i.e. finding a statistically ac-ceptable model associated with a false minimum, is a naturalconsequence of an unbiased modelling the process. (An inter-active i.e. human method may artificially avoid the problemby fixing ∆ α / α = throughout the model-building process).The correct solution is to repeat the whole fitting process atleast twice (and preferably more), searching for additionalsolutions with fewer velocity components and a lower AICc.We also note that, from calculations done so far, the numberof acceptable models (i.e. the degree of non-uniqueness) islikely to be reduced when fitting a temperature model (aswould be expected). In Bainbridge & Webb (2017a,b) and in the present paperwe have made extensive use of the AICc statistic (Hurvich& Tsai 1989) to select candidate models. In the course ofthis work, we also compared AICc results with those usingthe Bayesian Information Criterion, BIC (Bozdogan 1987).Both statistics add a penalty to the usual χ statistic, thevalue of which increases with increasing number of modelparameters. The purpose is to end up with a final modelhaving the “right” number of parameters. In other words,both AICc and BIC try to limit the model complexity toavoid “over-fitting”. The AICc (Eq.(1)) and BIC penaltiesare nk /( n − k − ) and k log n , (4)where n , k are the number of data points and free parametersrespectively.Fig.10 illustrates the contribution to the penalty fromeach free parameter in the model. Placing the spectral fittingboundaries (i.e. defining n ) is to some extent arbitrary. Thusin this application to absorption line spectroscopy, a penaltythat is relatively insensitive to n is preferable. However, in-cluding pixels far from line centre is required if we want tofit continuum parameters. These two considerations aloneshow that AICc is preferable to BIC in our application.Nevertheless, even though we have adopted AICc inthis paper, it is not entirely suitable. The spectral fittingboundaries influence the AICc penalty and the former areill-defined. Moreover, for an unsaturated absorption line (forexample), pixels far from the line centre impact little on theline profile whilst the converse is true, yet each pixel acrossthe absorption line or complex carries the same weight inthe AICc penalty term. Another way of expressing this isthat each parameter can influence specific regions within thedataset (but have little or no impact elsewhere); the param-eters describing an absorption component at −
50 km/s inFig.3, for example, have virtually no influence on the modelintensity at +50 km/s. The AICc (and BIC) fail to takesuch an effect into account. AICc is thus not optimal forabsorption line spectroscopy and some other (new) statisticis needed.
10 20 30 40 50246810
Figure 10.
Penalty per free parameter (Eq.4 divided by k ) vs.the number of pixels per parameter, n / k . The blue (solid), green(dotted), and red (dashed) curves illustrate the BIC penalty termfor different k . At fixed n / k , BIC increases relatively rapidly asmore model parameters are introduced. At fixed n / k , BIC also in-creases as the spectral fitting region n is increased. The mageneta(continuous) and black (dot-dash) curves shows the AICc penaltyterm. The AICc penalty is insensitive to k and relatively insensi-tive to n / k . The vertical grey (dashed) line indicates n / k = / for the synthetic model of Fig.3. Finally, the details of model non-uniqueness (Section5) are likely to depend on the statistic used i.e. AICc orsomething else. These considerations are beyond the scopeof the present paper and will be addressed in subsequentwork.
In this paper we have extended previous work that com-bined a genetic algorithm with non-linear least-squares toautomatically model absorption line data. “Interactive” ab-sorption line fitting is subjective and generally not repro-ducible. The method presented here brings objectivity andrepeatability, important requirements for the increasinglyhigh-quality data from new and forthcoming spectroscopicfacilities like ESPRESSO/VLT and HIRES/ELT.An AI method provides an unbiased estimate of ∆ α / α .The same is not true of any interactive method that holds ∆ α / α = during the model construction phase, only to“switch on” ∆ α / α as a free parameter at the end, once themodel is essentially complete. A procedure of this sort pref-erentially selects local minima closest to ∆ α / α = .Synthetic spectra, based on a well-known intermediateredshift absorption system, have been created and used totest the method’s performance. We have shown that modeldevelopment goes wrong i.e the wrong velocity structureis obtained, if the wrong broadening mechanism is usedto model the observational data. High quality observationsof quasar absorption systems show that in general the linebroadening does indeed contain both thermal and turbulentcontributions simultaneously. However, much of the previ-ous work on estimating ∆ α / α has been done assuming a sin- MNRAS , 1–14 (2020) C. C. Lee et al gle broadening mechanism. We have shown that there arehighly undesirable consequences of this: the wrong velocitystructure may easily be obtained, resulting in far less reliableparameter estimates. Fig. 8 shows this clearly for ∆ α / α .The extent to which this model non-uniqueness issueapplies even when the correct model is used is not yet clear.It is possible, even likely, that χ –parameter space is suf-ficiently complex, with multiple local minima, that falseminima result even when solving for b th and b turb simul-taneously. If so, model non-uniqueness limits the precisionachievable from any single measurement i.e. it is not pos-sible to reach the theoretical statistical limit using a singleobservation of ∆ α / α , no matter how precise the wavelengthcalibration. The implication is that, by nature, the problemof measuring ∆ α / α is a statistical one, requiring a large sam-ple of measurements to render the non-uniqueness system-atic negligible. This principle should perhaps underpin allfuture scientific efforts to determine whether the fine struc-ture constant varies in time or space. ACKNOWLEDGEMENTS
CCL thanks the Royal Society for a Newton InternationalFellowship during the early stages of this work. JKW thanksthe John Templeton Foundation, the Department of AppliedMathematics and Theoretical Physics and the Institute ofAstronomy at Cambridge University for hospitality and sup-port, and Clare Hall for a Visiting Fellowship during thiswork.
REFERENCES
Bainbridge M. B., Webb J. K., 2017a, Universe, 3, 34Bainbridge M. B., Webb J. K., 2017b, MNRAS, 468, 1639Bozdogan H., 1987, Psychometrika, 52, 345Carswell R. F., Webb J. K., 2014, VPFIT: Voigt profile fittingprogram, Astrophysics Source Code Library (ascl:1408.015)Carswell R. F., Webb J. K., 2020, Bob Carswell’s homepage, https://people.ast.cam.ac.uk/~rfc/
ESO ELT team 2010, The science case for the European Ex-tremely Large Telescope: The next step in mankind’s questfor the Universe. ESO online,
ESO ELT team 2011, An Expanded View of the Universe Sci-ence with the European Extremely Large Telescope. ESO on-line,
Hook I., 2009, in Moorwood A., ed., Science with the VLT in theELT Era. Springer Netherlands, Dordrecht, pp 225–232Hurvich C. M., Tsai C.-L., 1989, Biometrika, 76, 297Liske J., Spyromillio J., Tamai R., 2014, Top Level Requirementsfor ELT-HIRES,
Marconi A., et al., 2016, in Ground-based and Airborne In-strumentation for Astronomy VI, eds. Christopher J. Evansand Luc Simard and Hideki Takami. SPIE, pp 676 –687, doi:10.1117/12.2231653, https://doi.org/10.1117/12.2231653
Moscato P., 1989, On Evolution, Search, Optimization, Ge-netic Algorithms and Martial Arts - Towards Memetic Al-gorithms, Technical Report Caltech Concurrent Computa- tion Program, report 826, http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.27.9474&rep=rep1&type=pdf
Pepe F., et al., 2014, Astronomische Nachrichten, 335, 8This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000