Inspecting spectra with sound: proof-of-concept & extension to
datacubes–Inspecting spectra with sound: proof-of-concept & extension to
datacubes

# Inspecting spectra with sound: proof-of-concept & extension to datacubes

James W. Trayford  1 C.M. Harrison E-mail<EMAIL_ADDRESS>2 R.C.
Hinz E-mail<EMAIL_ADDRESS>2 M. Kavanagh Blatt  1 S.
Dougherty  2 & A. Girdhar2,3
1Institute of Cosmology and Gravitation  University of Portsmouth  Dennis
Sciama Building  Burnaby Road  Portsmouth PO1 3FX  UK
2School of Mathematics  Statistics and Physics  Newcastle University  NE1 7RU
UK
3Ludwig-Maximilians-Universität  Professor-Huber-Platz 2  D-80539 München
Germany

(Accepted XXX. Received YYY; in original form ZZZ; 2023)

###### Abstract

We present a novel approach to inspecting galaxy spectra using sound, via
their direct audio representation (‘spectral audification’). We discuss the
potential of this as a complement to (or stand-in for) visual approaches. We
surveyed 58 respondents who use the audio representation alone to rate 30
optical galaxy spectra with strong emission lines. Across three tests, each
focusing on different quantities measured from the spectra (signal-to-noise
ratio, emission-line width, & flux ratios), we find that user ratings are well
correlated with measured quantities. This demonstrates that physical
information can be independently gleaned from listening to spectral
audifications. We note the importance of context when rating these
sonifications, where the order examples are heard can influence responses.
Finally, we adapt the method used in this promising pilot study to spectral
datacubes. We suggest that audification allows efficient exploration of
complex, spatially-resolved spectral data.

###### keywords:

sonification – integral field spectroscopy – galaxies

##  1 Introduction

Sonification is the process of conveying data using non-verbal sound (Kramer
et al., 1999). Sound has the potential to provide unique advantages when
inspecting data. For example, the ear is particularly good at parsing evolving
or transient signals, with a larger dynamic range in frequency than the eye
(Heffner & Heffner, 2007; Sliney, 2016). Audio also provides an additional
modality to inspect multi-dimensional data, mitigating the limits of purely
visual methods. This is promising in the ‘big data’ era of astronomy, where
innovation is needed to explore complex, multidimensional data. Indeed,
sonification has long existed, finding application in science and engineering
Hermann et al. (2011), and seeing a surge of interest in astronomy in recent
years (Zanella et al., 2022).

Sound has become established in astronomy as a means to provide educational
resources and accessible channels for visually impaired (VI) individuals,
either in pure audio or multi-modal context, often by using sound to accompany
and enhance visualisations (Alexander et al., 2011; Harrison et al., 2022;
Bardelli et al., 2022; García-Benito & Pérez-Montero, 2022; Zanella et al.,
2022). However, work is on going to demonstrate the utility of sonification as
a research tool. Diaz-Merced (2013) and Tucker Brown et al. (2022) provide
precedent, demonstrating how signals can be discerned from sonified
lightcurves in both multimodal (sound & visual) and standalone contexts,
through controlled user tests.

In this pilot study, we further explore the research utility of sonification.
Particularly, we prototype applications to (hyper-)spectral datacubes, a form
of data that is transforming observational astronomy with the rise of resolved
spectroscopy and interferometry (e.g. Bacon et al., 2010; Böker et al., 2022),
and across other fields from biomedical engineering to geophysics (e.g.
Cassidy et al., 2004; Bernhardt et al., 2007; Li et al., 2013). While the
detail and dimensionality (i.e. 3D) of these datacubes makes purely visual
inspection difficult, we posit that the spectral dimension can naturally be
represented with sound, accompanying the two spatial dimensions typically
projected for visual display. This could enable efficient datacube exploration
& build intuition for otherwise unwieldy data.

We utilise the STRAUSS111Sonification Tools and Resources for Astronomers
Using Sound Synthesis. code to sonify each galaxy spectrum in our experiments
(Trayford, 2021). STRAUSS is a modular python package providing diverse
sonification options for data analysis and presentation, allowing parameters
of the data to be directly mapped to expressive properties of sound (e.g.
pitch, stereo field placement, timbre, etc), via synthesis or manipulation of
recorded sounds, or ‘samples’ (applied in Harrison et al., 2022; Tucker Brown
et al., 2022).

In astronomy-associated tools for sonifying 1D data, i.e. where $y=f(x)$,
typically the $y$-variable (e.g. brightness) is mapped to a sound property
such as pitch, and its variation with the $x$-variable (e.g. date) is
represented by the modulation of a carrier tone over time (see e.g.
astronify222astronify.readthedocs.io). However, a limitation of the parameter
modulation approach is the need to ‘play out’ the entire $x$-domain to fully
inspect the data. Playback must therefore be slow enough for subtle or small-
scale variations in the data to be perceptible. This is particularly pertinent
for narrow, highly informative features, like spectral lines. Audio inspection
of large data sets (or even a single modern datacube) can therefore be a slow
process. An alternative is Spectral Audification (SA, Sturm, 2005; Newbold et
al., 2015), harnessing the analogy between electromagnetic and audio spectra.
By scaling the flux as a function of electromagnetic wavelength to an audible
frequency range, an audio signal can be synthesised, representing this
distribution of relative power across perceptible frequencies. This provides
an ‘instantaneous’ representation of the spectrum where various frequency
components can be heard simultaneously. This can be considered a form of
‘spectral compression’.

We note that despite it’s direct or intuitive relation to the data, such
‘audification’ isn’t always necessarily the most effective, with more
abstracted sonified representations sometimes superior at conveying data (e.g.
Franinovic & Serafin, 2013). Data sonification innovations have emerged across
many data fields, both in analytic contexts, as well as for more outreach-
based or artistic applications. Sonification of multidimensional medical data
for diagnostic purposes exemplify analytic applications, including
electroencephalographic data (Väljamäe et al., 2013) and optical coherence
tomography data (Ahmad et al., 2010). Audio representations of spectral data
have been demonstrated through a similar SA approach (e.g. nuclear magnetic
resonances Morawitz 2019, X-ray scattering data Schedel & Yager 2012, the CMB
angular power spectrum McGee et al. 2011) but also abstractly as e.g. feature
sets via autoencoder training on stellar spectral libraries (Riber & Garcia,
2023). Such approaches may find use in more analytic or quantitative contexts.
Here, we focus on a spectral audification approach in particular.

We propose this SA approach as a promising method (in both multimodal and
standalone contexts) for inspecting spectroscopic data, particularly in the
context of spectral datacubes. In this study, we sonify individual galaxy
spectra representing a range of physical properties, and evaluate how well
their relative properties can be gleaned from SA sonification via a series of
rating tasks, comprising a survey. § 2 details the sonification approach, § 3
details the spectra selection and survey design, § 4 then presents the results
from the 58 respondents (4.1), and the broader applications to spectral
datacubes (4.2). We conclude in § 5. For clarity, we refer to light spectra
using their wavelength, and sound using frequency.

##  2 Spectral Audification

We synthesise the audio representing each spectrum additively, with a
superposition of sinusoidal waves. A sine wave represents each wavelength,
$\lambda_{i}$, of the galaxy spectra, with amplitude corresponding to the
observed flux density, $F_{i}$. These waves are summed, following

$A(t)=\sum_{i}F_{i}\sin\left(2\uppi\nu_{i}t+\phi_{i}\right),$ (1)

where $t$ is time in seconds, $\phi_{i}$ is a $[0,2\uppi]$ uniformly-
distributed random phase (important to avoid imposing spurious phase
correlations), and $\nu_{i}$ is the audio frequency mapped from each
wavelength of the galaxy spectrum, with

$\nu_{i}=\frac{\lambda^{-}\left(\lambda^{+}-\lambda_{i}\right)\left(\nu^{+}-\nu^{-}\right)}{\lambda_{i}\left(\lambda^{+}-\lambda^{-}\right)}+\nu^{-},$
(2)

where $+$ and $-$ superscripts indicate maxima and minima of value ranges. We
chose a minimum audio frequency of $\nu^{-}$ = 250 Hz with a maximum of
$\nu^{+}$ = 1250 Hz. This range was motivated by: (1) a range typical human
hearing has an almost constant level of frequency resolution (i.e. the ability
to audibly distinguish between two different frequencies, Zwicker et al.,
1957) and; (2) there is not dramatic changes in perceived loudness at a fixed
pressure level across the chosen range for pure (i.e. monochromatic) tones
(Suzuki & Takeshima, 2004). Our complex tones, could exhibit differing
behaviour, however this is taken as a useful guide. These considerations also
helped to ensure the sounds were in a comfortable frequency range, avoiding
the highest frequencies that can be particularly uncomfortable or tiring to
listen to. We produced SAs of 3 s duration, allowing participants to replay
the audio (see § 3.1).

We note that this additive synthesis approach affords maximal flexibility,
allowing arbitrary wavelength ranges and spacings to be represented,
particularly useful for the explorative prototyping of this pilot study.
However this approach is typically much slower than using an Inverse Fast
Fourier Transform (IFFT). In § 4.2 we discuss employing IFFT to produce SAs
orders-of-magnitude faster.

##  3 Participant testing and spectra used for spectral audification

Our study was designed to assess if participants ranked spectra on what they
heard from hearing a SA of the spectra, then these ratings would correspond to
the emission-line properties. Across three separate tests, detailed below, we
assessed the participants’ ability to hear variations in three emission-line
properties: (a) signal-to-noise ratio (SNR); (b) emission-line width and; (c)
emission-line flux ratios. The survey design is presented in § 3.1 and the
galaxy spectra used for each test is presented in § 3.2. A transcript of the
survey form and the sounds presented to participants are made available333See
doi.org/10.25405/data.ncl.22816442.

Figure 1: Visual representations of galaxy spectra that were played as audio
files to the participants during the initial training (the visuals were not
shown). The three columns correspond to the three tests (SNR, emission-line
width, and emission-line ratio). In each case, the top panel shows an example
of Rating 1 and the bottom, an example of Rating 5. The audio files can be
found in the supplementary material.

###  3.1 Survey design

We performed an anonymous participant survey via Google Forms, using embedded
video (with blank visuals) to play sounds. Participants were not informed that
the sounds corresponded to galaxy spectra, but simply that they survey was
part of “exploring new ways to represent astronomy with sound.” We also
collected data on demographics (musical and data analysis experience,
optionally age) and user experience, but forego its analysis in this pilot
study, given our limited sample size.

We ran three separate tests, each testing variation in a different quantity of
interest (QoI, see § 3.2). At the beginning of each of the tests, six example
sounds were presented to participants, functioning as a tutorial: three
examples of the very low QoI values (’Rating 1’) and three examples of very
high values (’Rating 5’). For each test, the participants were then asked to
rate a set of 10 sounds (presented in pseudo-random order) with an integer
value of 1–5, inclusive. Of these 10 sounds, one was a duplicate (i.e., each
test contained 9 unique sounds, of which 1 was presented to the participants
twice), where we avoid duplicates being presented consecutively. The survey
was run from 8th – 29th March 2022, with 58 responses. Of these respondents,
26 were aged $<$ 30, 19 aged 30-50 and 10 aged $>$ 50\. Users also exhibited
diversity in their self-assessed experience as musicians (data analysts) with
20 (37) Experienced or Very Experienced, 30 (14) with some level of experience
and 8 (7) with No Experience.

###  3.2 Galaxy spectra

We opted to use spectra of real galaxies (as opposed to synthetic or idealised
data). All spectra and corresponding automated emission-line properties, are
taken from the Sloan Digital Sky Survey (SDSS) Data Release 7 (Abazajian et
al., 2009). We focused on emission-line properties, using continuum-subtracted
spectra. Additionally, we select a narrow redshift range of $z=0.09$–0.11 and
required a signal-to-noise ratio (SNR) of $>$3 for each emission line of
interest (H$\beta$, [O iii], H$\alpha$, and [N ii]). We shifted all spectra to
the rest-frame, using the spectroscopic redshifts provided by SDSS.

Test A: Variations in SNR: For this sub-sample we were interested if
participants could successfully rate the SNR of spectral emission lines from
their SAs alone. To minimise other factors (than SNR) changing the sound, we
excluded Type 1 active galactic nuclei (AGN) and used a relatively narrow
emission-line flux ratio selection, with
$\log_{10}({\rm[O~{}III]/H\beta})\in[0.3,0.6]$ and
$\log_{10}({\rm[N~{}II]/H\alpha})\in[-0.7,-0.4]$. We calculated an average,
indicative, SNR for each spectra by taking an average over the four emission
lines: H$\beta$, [O iii], H$\alpha$, and [N ii]. For the six examples
presented during the participant training, we used three galaxy spectra with
low SNRs of 9–13, which we described as ‘Very Noisy Signal’ (Rating 1) and
three galaxy spectra with high SNRs of 84–122, described as ‘Very Clear
Signal’ (Rating 5, see first column of Fig. 1). For the nine cases used in
participant testing, we randomly selected spectra covering a SNR range 11–115.
The wavelength range used for their sonification was 400-800 nm, presented in
the left column of Fig. 1.

Test B: Variations in emission-line width: For this sub-sample we were
interested if participants could successfully rate the H$\alpha$ emission-line
width of spectra, from SA alone. We used a wavelength range of 650–665 nm
(i.e., centred around H$\alpha$ and [N ii]), presented in the central column
of Fig. 1. To minimise the influence of varying emission-line flux ratios, we
only used spectra in a narrow range of
$\log_{10}({\rm[N~{}II]/H\alpha})\in[-0.7,-0.4]$. For the six examples
presented during the participant training, we used three spectra with low
values of full-width at half maximum (FWHM) for H$\alpha$ (FWHM = 220–240 km
s-1), which we label as ‘Type B’ (Rating 5) and three examples with large
widths of H$\alpha$ FWHM=2400–2900 km s-1, labelled as ‘Type A’ (Rating 1). We
note that this labelling was chosen in lieu of clear shared vocabulary to
describe the distinct qualities of the sounds. This results in a dichotimised
rating system for this particular test, where 1 and 5 representing ‘high
confidence’ of Type A and B, respectively, with 3 indeterminate (revisited in
§ 4).

Test C: Variations in emission-line flux ratio: For this sub-sample we were
interested if participants could successfully rate [O iii]/H$\beta$ emission
line flux ratio in spectra, from their SA alone. To minimise the influence of
emission-line width variations on the sound, we excluded Type 1 AGN and
required that the FWHM of the [O iii] line be less than 250 km s-1. We also
required the H$\beta$ SNR to be $>$50\. For the training examples, we used
three galaxy spectra with high [O iii]/H$\beta$ ratios ($\log$[O
iii]/H$\beta\in[0.69,0.79]$), described as ‘Low Pitch Dominated’ (Rating 1)
and three examples with low ratios ($\log$[O iii]/H$\beta\in[-0.48,-0.27]$),
which we described as ‘High Pitch Dominated’ (Rating 5). The wavelength range
used in the sonification was 480–504 nm, with examples shown in the right
column of Fig. 1.

Figure 2: Mean participant rating (error bars show standard error on the mean)
for each of the galaxy SAs, as a function of the rank ordering of the
principle physical quantify varied in each test. Panels correspond to
different test (top: SNR; middle: emission-line width; bottom: emission-line
ratio). The physical quantity is indicated in the top-left of each panel, with
precise values inset along the $x$-axis (with units where appropriate). The
background violin plots represent the full distribution of participant
rankings. Orange dotted lines highlight the pair of identical examples for
each experiment. As their values are identical, relative ranking instead
represents the order they appear in the survey (earlier to later). The text
annotating data points indicates their presentation order to respondents.

##  4 Results and Discussion

###  4.1 Rating results

Our primary survey results are summarised in Fig. 2, with panels a) to c) (top
to bottom) corresponding to the three tests detailed in § 3.2. The mean
participant rating for each galaxy SA is plotted in rank order by QoI, such
that user rating and physical ranking should positively correlate if ratings
are representative. In all three tests, the average ratings show a clear
positive correlation, particularly between the extrema. In all cases the
highest and lowest ranked examples are rated as such by the plurality of
users. We regard this outcome as a very encouraging sign for the technique,
particularly given our rudimentary rating tasks over just 58 respondents with
limited training and knowledge about how the SA approach works.

Given the nascent stage of this approach, we wish to avoid over-interpretation
of these first results. However, we perform some basic numerical analysis,
focusing on rank correlation (monotonicity) for its insensitivity to the
property spacing and perceived scaling of the QoI from the SAs. We also
exclude the repeated SAs (right-hand distributions in orange boxes) from this
correlation for now, as their relative rank is unrelated to the QoI. The mean
ratings correlate well with physical rank for all tests, with Spearman’s rank
coefficients of 0.86, 0.9 and 0.95 for tests a) to c), respectively. When
applied to individual ratings (constituting the underlying distributions in
Fig. 2), these correlations are weaker but still positive at 0.61, 0.72 and
0.69, respectively. However, rank correlation of individual responses are
affected by the limited sample size and heavy quantisation of ratings and
rankings. To quantify this uncertainty, we compute the rank correlation
coefficient bounds for quantised data, proposed by Couso et al. 2018 (based on
Kendall’s $\tau$, a similar metric to Spearman’s). We recover [0.13, 0.71],
[0.41, 0.79] and [0.2, 0.76] for tests a to c, respectively. Definitively
positive correlations, but indicating that correlation strength is poorly
constrained given survey limitations.

Clearly there is some underlying user-user variance to investigate further.
For example, while the line ratio test (c) is most accurate in terms of
average rating, it is less precise than the line width test (b), where the
three highest and lowest rank extrema are reliably rated at high and low
extremes by users, respectively. This could partly be due to how questions are
posed; test c) asks users to assess the relative strength of high or low
pitches in the sound, while test b) asks their confidence in sorting examples
between Types A and B (broad to narrow H$\alpha$). Test b) respondents
correctly rate the 3 broadest and narrowest H$\alpha$ galaxies with low and
high ratings, respectively. The 4 intermediate H$\alpha$ galaxies are found to
be relatively indeterminate. We suggest respondents are treating test b) as
more of a sorting task, yielding more ‘stepped’ mean ratings corresponding to
A, B and indeterminate types, contrasting the continuous variation for test
c).

Another aspect that may influence respondents is the (pseudo-random) order SAs
are presented in. As a simple test, we rank examples by mean user rating
(including repeats), and measure how these ranks are offset from their
physical property ranking. We then compare this offset to the difference in
physical rank between a SA from the one heard directly before. We find that
these quantities show a significant (2$\sigma$) rank correlation for test a) -
users are more likely to (under-) over-rate the clarity of the signal if they
heard a relatively (clear) noisy example before. Test b) only correlates at
the 1$\sigma$ level, and test c) shows no evidence of correlation444For test
c), concordance in the ordering by physical rank and mean rating (i.e. their
good correlation) means there is insufficient data to measure this.. These
findings are reflected by the repeat examples; for tests a) and b),
respondents are more likely to rate a repeat example higher if the preceding
SA is ranked relatively lower, excepting test c (which incidentally shows the
smallest absolute difference between the mean rating of repeats). While more
investigation is needed, it is evident that context can influence how users
perceive SAs, with users often sensitive to contrast, reacting against the SA
they heard before. This sensitivity may have advantages, e.g. when inspecting
continuously evolving sound or where differential changes are interesting.
This is particularly relevant to our prototype spectral datacube ‘listener’,
next. Particular sensitivity to subtle relative differences in auditory
stimuli has been established in psychoacoustics, conceptualised as a ‘just
noticeable difference’ (Zwicker & Fastl, 2007; Ahmad et al., 2010).

Respondents also rated the tolerability of the SA examples, from 1 (Very
intolerable) to 5 (Tolerable). These ratings were chosen 1, 7, 11, 26, and 13
times respectively, suggesting users generally found these on the tolerable
side, though there was a minority finding the sounds difficult to listen to.

###  4.2 Future applications to datacube exploration

We illustrate our prototype datacube explorer in Fig. 3 (linking to sonified
animation in caption), using a MUSE (Bacon et al., 2010) datacube of galaxy
J1316+1753, a $z=0.15$ quasar host with defined spiral arms (studied in
Girdhar et al., 2022). Here, the SA is generated from the aggregate spectrum
in an aperture, freely moved around the 2 projected spatial dimensions of the
datacube. We render a figure-eight path about the galaxy center, sonifying the
wavelength region around the rest-frame H$\alpha$ complex. Doppler shifted
emission can be heard as a distinct pitch-shift between the emission either
side of the disc (sites 1 to 3). In the central region (site 2) the dominance
of broad lines leads to a ‘noisier’ tone. Further detail can be heard, e.g.
the narrow emission and higher recession velocities around site 3.

Ideally, a spectral datacube SA tool would be ‘real-time’, allowing aperture
size, and wavelength parameters to be adjusted on the fly. This requires low-
latency SA generation. To this end we have implemented a dedicated
‘spectraliser’ generator type in STRAUSS, via an IFFT algorithm555python
notebooks using this new spectraliser approach in STRAUSS can be found at
github.com/james-trayford/AudibleUniverseWorkbooks/tree/group3. To generate
the 3 s clips of our study (which could be looped and faded between based on
position), this approach takes an average of 2.3 ms, 4 orders of magnitude
faster than the additive approach, and crucially, safely faster than the
maximum tolerable latency for interactive applications (Lago & Kon, 2004).
This augurs well for real-time, responsive SA of spectral datacubes in the
future666Such a tool would complement existing experiments and tools for
resolved astronomical data, see e.g. the Soniverse and VoxMagellan projects.

Figure 3: Prototype audification of a spectral datacube. Top panel shows
H${\rm\alpha}$ emisson, overlaid with three aperture positions along a figure-
eight track through the galaxy center (grey dashed). Bottom panel shows the
continuum-subtracted and normalised flux for each aperture position, over a
wavelength range capturing the H$\alpha$ line in the overall galaxy rest-
frame. The upper $x$-axis illustrate how this spectrum is directly mapped to
audio frequency. A full sonified animation can be found at
https://doi.org/10.25405/data.ncl.22816442.

##  5 Conclusions & Future Directions

We present spectral audification (SA) as a novel, auditory approach for
inspecting astronomical spectra, to complement and enhance traditionally
visual approaches in research, while also providing an accessible channel for
VI people to inspect data. Using galaxy spectra (Fig. 1), we demonstrate how
key properties (e.g. SNR, linewidth, line ratios) can be gleaned from their
audio representation alone, via a 58-respondent survey. We show that user
ratings correlate appropriately with their physical ranking in Fig. 2; very
encouraging in the context of our small sample and minimal training of users.
User-user variance is noted, and partly attributed to the role of context (SA
order), finding that mean user ranking offsets can correlate with how the QoI
has changed from the spectrum heard immediately prior (particularly for SNR, §
3.1). We note that this sensitivity to contrast may be useful in the context
of our prototype datacube application. We provide an illustration and sonified
animation of a prototype SA explorer for spectral datacubes (Fig. 3), noting
that live sonification of spectral datacubes are feasible given our modified
IFFT approach.

In future surveys, improving our perceptual evaluation methods (Misdariis et
al., 2022), and larger sample sizes could give quantitative insight into the
reliability of SA inspection, beyond our proof-of-concept (see § 4.1). The
role of context could probed further, for example, by randomising SA order by
user (as e.g. a Zooniverse project). Tucker Brown et al. (2022) note the
importance of training when introducing sonification approaches for research.
In this work, we provided limited training, also deliberately giving users no
insight or explanation of the technique (§3.1). Complementary, extended
evaluation studies where users can develop more knowledge and experience
around the approach (particularly for the datacube explorer) may better
reflect the potential of audio modalities; after all, most researcher’s visual
inspection skills are honed over decades.

More extensive testing is also important to understand how the training rate
and overall effectiveness of SA inspection may vary statistically with
demographic factors, including age, experience and level of vision. Similarly,
comparisons between the effectiveness of visual, auditory and hybrid
inspection modes (i.e. additional testing including visual representations of
spectra in 1D or as spectrograms) could be informative. Developing effective
training approaches and materials for SA (as e.g. Alexander et al., 2011),
alongside improved sonification techniques, appears key to better outcomes
from auditory inspection of data. We also identify a number of avenues for
investigation to improve the aesthetics of future SA representations.

In the short-term, we will work on the spectral datacube explorer (§4.2),
seeing great potential for data intuition, accessibility and discovery. With
the proliferation of (hyper)spectral datacubes, enhanced inspection modes
could improve efficiency in research. The discovery potential of SA, where
certain aspects may be clearer or enhanced through sound, could offer
invaluable unforeseen benefits (e.g. finding signals or artefacts that could
be otherwise overlooked in the confounding detail of modern spectral
datacubes). Longer term, integration of datacube listeners into online data
archives (e.g. MAST777https://mast.stsci.edu/viz/ui/) is a goal. By
demonstrating the efficacy of analytical sonification in astronomical data, we
could also generalise our methods to other fields involving spectral
(datacube) analysis (e.g. medicine & engineering Cassidy et al., 2004;
Bernhardt et al., 2007; Li et al., 2013), complementing prior work on audio
representations and their practical efficacy in those data fields (e.g. Ahmad
et al., 2010; Lee et al., 2022).

## Data Availability

A survey transcript and all of the audio files used during the tutorials and
testing, along with our animated prototype of a sonified IFU datacube, can be
found here: https://doi.org/10.25405/data.ncl.22816442.

## Acknowledgements

We thank the anonymous referees for a their constructive feedback which helped
to improve the manuscript, and suggesting the term ‘spectral compression’. We
also thank all of the anonymous, volunteer participants who carried out the
surveys. We thank the Lorentz Centre for hosting the Audible Universe workshop
that partly inspired this work and Scott Fleming for encouraging and practical
discussion of how to integrate an audio IFU explorer into MAST. We also thank
Stephen Barrass for useful discussions. JWT acknowledges support via the STFC
Early Stage Research & Development Grant, reference ST/X004651/1. CMH
acknowledges funding from an United Kingdom Research and Innovation grant
(code: MR/V022830/1).

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Ahmad et al. (2010) Ahmad A., Adie S. G., Wang M., Boppart S. A., 2010, Optics Express, 18, 9934
  * Alexander et al. (2011) Alexander R., Gilbert J., Landi E., Simoni M., Zurbuchen T., Roberts D., 2011, in 17th International Conference on Auditory Display (ICAD2011).
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508 (arXiv:2211.16795), doi:10.1117/12.856027
  * Bardelli et al. (2022) Bardelli S., Ferretti C., Ludovico L. A., Presti G., Rinaldi M., 2022, in Proceedings of "Interactive Cultural Heritage and Arts", 23rd HCI International Conference. (arXiv:2202.05539), doi:10.1007/978-3-030-77411-0_12
  * Bernhardt et al. (2007) Bernhardt M., Cowell C., Oxford W., 2007, Proceedings of SPIE - The International Society for Optical Engineering
  * Böker et al. (2022) Böker T., et al., 2022, A&A, 661, A82
  * Cassidy et al. (2004) Cassidy R., Berger J., Lee K., Maggioni M., Coifman R., 2004, in Conference Record of the Thirty-Eighth Asilomar Conference on Signals, Systems and Computers, 2004.. pp 1611–1615 Vol.2, doi:10.1109/ACSSC.2004.1399429
  * Couso et al. (2018) Couso I., Strauss O., Saulnier H., 2018, Fuzzy Sets and Systems, 343, 50
  * Diaz-Merced (2013) Diaz-Merced W., 2013, PhD thesis, School of Computing Science, University of Glasgow
  * Franinovic & Serafin (2013) Franinovic K., Serafin S., 2013, Sonic interaction design
  * García-Benito & Pérez-Montero (2022) García-Benito R., Pérez-Montero E., 2022, in Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 28–33 (arXiv:2205.12984), doi:10.22201/ia.14052059p.2022.54.06
  * Girdhar et al. (2022) Girdhar A., et al., 2022, MNRAS, 512, 1608
  * Harrison et al. (2022) Harrison C., Trayford J., Harrison L., Bonne N., 2022, Astronomy and Geophysics, 63, 2.38
  * Heffner & Heffner (2007) Heffner H. E., Heffner R. S., 2007, Journal of the American Association for Laboratory Animal Science, 46, 20
  * Hermann et al. (2011) Hermann T., Hunt A., Neuhoff J. G., 2011, The Sonification Handbook. Logos Verlag, Berlin
  * Kramer et al. (1999) Kramer G., et al., 1999, in The sonification report: Status of the field and research agenda. Report prepared for the National Science Foundation by members of the International Community for Auditory Display.
  * Lago & Kon (2004) Lago N. P., Kon F., 2004, in vol. 2004 of ICMC. http://hdl.handle.net/2027/spo.bbp2372.2004.142
  * Lee et al. (2022) Lee S.-W., Lee S. H., Cheng Z., Yeo W. S., 2022, Diagnostics, 12, 1728
  * Li et al. (2013) Li Q., He X., Wang Y., Liu H., Xu D., Guo F., 2013, Journal of biomedical optics, 18, 100901
  * McGee et al. (2011) McGee R., Van der Veen J., Wright M., Kuchera-Morin J., Alper B., Lubin P., 2011, International Conference for Auditory Display
  * Misdariis et al. (2022) Misdariis N., Özcan E., Grassi M., Pauletto S., Barrass S., Bresin R., Susini P., 2022, Nature Astronomy, 6, 1249
  * Morawitz (2019) Morawitz F., 2019, International Conference on Auditory Display, pp 169–175
  * Newbold et al. (2015) Newbold J. W., Hunt A., Brereton J., 2015, in 21st International Conference on Auditory Display (ICAD2015). pp 329–330
  * Riber & Garcia (2023) Riber A. G., Garcia F. S., EasyChair, 2023, AI-Rmonies of the Spheres, EasyChair Preprint no. 9890
  * Schedel & Yager (2012) Schedel M., Yager K. G., 2012, International Conference on Auditory Display, pp 95–98
  * Sliney (2016) Sliney D., 2016, Eye, 30, 222
  * Sturm (2005) Sturm B. L., 2005, Leonardo, 38, 143
  * Suzuki & Takeshima (2004) Suzuki Y., Takeshima H., 2004, The Journal of the Acoustical Society of America, 116, 918
  * Trayford (2021) Trayford J., 2021, james-trayford/strauss: v0.1.0 Pre-release, doi:10.5281/zenodo.5776280
  * Tucker Brown et al. (2022) Tucker Brown J., Harrison C. M., Zanella A., Trayford J., 2022, MNRAS, 516, 5674
  * Väljamäe et al. (2013) Väljamäe A., Steffert T., Holland S., Marimon X., Benitez R., Mealla S., Oliveira A., Jordà S., 2013, International Conference on Auditory Display, pp 85–93
  * Zanella et al. (2022) Zanella A., Harrison C. M., Lenzi S., Cooke J., Damsma P., Fleming S. W., 2022, Nature Astronomy, 6, 1241
  * Zwicker & Fastl (2007) Zwicker E., Fastl H., 2007, Psychoacoustics: Facts and models. Vol. 22, Springer Science & Business Media
  * Zwicker et al. (1957) Zwicker E., Flottorp G., Stevens S. S., 1957, The Journal of the Acoustical Society of America, 29, 548


#  Even order contributions to relative energies vanish for antisymmetric
perturbations

O. Anatole von Lilienfeld<EMAIL_ADDRESS>Vector Institute
for Artificial Intelligence, Toronto, ON, M5S 1M1, Canada  Departments of
Chemistry, Materials Science and Engineering, and Physics, University of
Toronto, St. George Campus, Toronto, ON, Canada  Machine Learning Group,
Technische Universität Berlin and Institute for the Foundations of Learning
and Data, 10587 Berlin, Germany  Giorgio Domenichini  University of Vienna,
Faculty of Physics, Kolingasse 14-16, AT-1090 Vienna, Austria  University of
Vienna, Vienna Doctoral School in Physics, Boltzmanngasse 5, AT-1090 Vienna,
Austria

###### Abstract

We show that even order contributions to energy differences between any two
iso-electronic compounds vanish when using perturbation theory around an
averaged electronic reference Hamiltonian. This finding generalizes the
previously introduced alchemical chirality concept [von Rudorff, von
Lilienfeld, Science Advances, 7 2021] by lifting the symmetry requirements for
transmutating atoms in the iso-electronic reference system. The leading order
term corresponds to twice the Hellmann-Feynman derivative evaluated using the
electron density of the averaged Hamiltonian. Analogous analysis reveals Mel
Levy’s formula for relative energies [J. Chem. Phys. 70, 1573 (1979)] to
include the first order contribution while overestimating the higher odd order
energy contributions by a factor linearly increasing in order. Using density
functional theory, we illustrate the predictive power of the leading order
term for estimating relative energies among diatomics in the charge-neutral
iso-electronic 14 proton series N2, CO, BF, BeNe, LiNa, HeMg, HAl, and the
united atom, Si. The framework’s potential for the simultaneous exploration of
multiple dimensions in chemical space is demonstrated for toluene by
evaluating relative energies between all the possible 35 antisymmetric BN
doped isomers (dubbed “alchemical diastereomers”). Based solely on toluene’s
electron density, necessary to evaluate all the respective Hellmann-Feynman
derivatives, mean absolute errors of predicted total potential energy
differences between the alchemical diastereomers are on the scale of mHa.

##  I Introduction

Chemical compound space (CCS) is vast Kirkpatrick and Ellis (2004), and
deepening a physics based understanding of its structure is of fundamental
importance to the chemical and materials sciences. It can also be beneficial
for accelerating the discovery and design of materials and molecules. The
calculation of approximate solutions to the electronic Schrödinger equation,
most frequently obtained by numerically solving the variational problem for an
approximated expectation value of the electronic Hamiltonian, constitutes one
of the major bottleneck when pursuing this goal. A possible use case of such
solutions involves the estimation of relative energies. Alchemical
perturbation density functional theory (APDFT) von Rudorff and von Lilienfeld
(2020) represents a computationally less demanding approach to estimate
relative energies. APDFT relies on the continuous interpolation of the
external potential as already previously discussed and studied, for example by
Foldy Foldy (1951), Wilson E. B. Wilson, Jr. (1962), Levy Levy (1978a, 1979),
or Politzer and Parr Politzer and Parr (1974). Related research dealing with
alchemical changes includes, among others, Refs. Weigend  _et al._ (2004); von
Lilienfeld  _et al._ (2005); Wang  _et al._ (2006); Beste  _et al._ (2006);
Sheppard  _et al._ (2010), and more recently Refs. Lesiuk  _et al._ (2012);
Balawender  _et al._ (2013); Weigend (2014); Chang  _et al._ (2016); Fias  _et
al._ (2017); Munoz and Cardenas (2017); Fias  _et al._ (2018); Saravanan  _et
al._ (2017); Al-Hamdani  _et al._ (2017); Balawender  _et al._ (2018); Chang
and von Lilienfeld (2018); Griego  _et al._ (2018); von Rudorff and von
Lilienfeld (2019); Balawender  _et al._ (2019); Shiraogawa and Ehara (2020);
Griego  _et al._ (2020a, b); von RAtomicAPDFTudorff and von Lilienfeld (2020);
Muñoz  _et al._ (2020); Griego  _et al._ (2021); Gómez  _et al._ (2021); Eikey
_et al._ (2022a); Shiraogawa and Hasegawa (2022); Shiraogawa  _et al._ (2022);
Eikey  _et al._ (2022b); Shiraogawa and Hasegawa (2023); Balawender and
Geerlings (2023).

The concept of alchemical chirality introduces approximate electronic energy
degeneracies among seemingly unrelated pairs of iso-electronic compounds
(“alchemical enantiomers”) von Rudorff and von Lilienfeld (2021) that share
the same geometry but differ in constitution or composition. For two compounds
to correspond to such alchemical enantiomers (and exhibit approximate
degeneracy in electronic energy) their external potentials have to average to
an external reference potential exhibiting such symmetry that each atom being
alchemically transmutated happens to have the same chemical environment. In
this paper, we discuss what happens if this symmetry requirement for the
transmutating atoms is lifted. In particular, we consider anti-symmetric
alchemical iso-electronic perturbations of any electronic reference
Hamiltonian. The effect of lifting the alchemical chirality symmetry
requirement is that the first order term does no longer vanish. As we will see
below, all even order contributions, however, vanish, and the alchemical
chirality case is identified as the special cases where the first order
contribution disappears due to the high symmetry in the reference system
leading to the exact cancellation of odd perturbing potential with even
electron density in the Hellmann-Feynman derivative.

##  II Theory

Consider two iso-electronic compounds $l$ and $r$ whose electronic Hamiltonian
only differs in their respective external potentials, $\Delta v=v_{r}-v_{l}$.
Corresponding changes of the non-relativistic concave electronic energy can be
obtained through $\Delta
E=E_{r}-E_{l}=\int_{-1}^{+1}d\lambda\partial_{\lambda}E$ using a single one-
dimensional coupling parameter,
$-1\leq\lambda\leq 1,\lambda\in\mathbb{R}$, in the Hamiltonian,
$\hat{H}(\lambda)$. Assuming a linear interpolation and defining the
Hamiltonian such that the extreme values of $\lambda$ correspond to the two
respective compounds,
$\hat{H}(\lambda)=(\hat{H}_{r}+\hat{H}_{l})/2+\lambda(\hat{H}_{r}-\hat{H}_{l})/2$,
we show below that the averaged mid-point Hamiltonian,
$\hat{H}(\lambda=0)=(\hat{H}_{r}+\hat{H}_{l})/2=:\hat{H}_{0}$, represents a
reference system of remarkable interest: It is a pivot point for perturbations
in chemical compound space. Consider the linear iso-electronic variation of
the external potential, $\hat{H}^{\prime}:=(\hat{H}_{r}-\hat{H}_{l})/2=\delta
v$, one can thus think of
$\hat{H}(\lambda)=\hat{H}_{0}+\lambda\hat{H}^{\prime}$ as an antisymmetric
perturbation of the average and along either direction of the path connecting
compounds $r$ and $l$, exactly at their mid-point. Setting $\lambda=-1,+1$,
note how the corresponding two potentials, $v_{l},v_{r}$ are respectively
recovered, i.e. $2\hat{H}^{\prime}=2\delta
v=(v_{r}-v_{0})+(v_{0}-v_{l})=\Delta v$.

Assuming convergence, we can now expand the respective electronic energy as a
generic perturbation series using the averaged Hamiltonian as reference,
$E_{0}=\langle\hat{H}_{0}\rangle$, and involving antisymmetric variations
towards positive or negative changes in $\lambda$, i.e. $E_{r}:=E(\lambda=+1)$
and $E_{l}:=E(\lambda=-1)$. More specifically,

$\displaystyle E_{r}$ $\displaystyle=$ $\displaystyle
E_{0}+\partial_{\lambda}E_{0}+\frac{1}{2}\partial^{2}_{\lambda}E_{0}+\frac{1}{6}\partial^{3}_{\lambda}E_{0}+\frac{1}{24}\partial^{4}_{\lambda}E_{0}+\frac{1}{120}\partial^{5}_{\lambda}E_{0}+...$
$\displaystyle E_{l}$ $\displaystyle=$ $\displaystyle
E_{0}-\partial_{\lambda}E_{0}+\frac{1}{2}\partial^{2}_{\lambda}E_{0}-\frac{1}{6}\partial^{3}_{\lambda}E_{0}+\frac{1}{24}\partial^{4}_{\lambda}E_{0}-\frac{1}{120}\partial^{5}_{\lambda}E_{0}+....$

where $\Delta\lambda$ has been set to -1 and +1 to calculate $E_{l}$ and
$E_{r}$, respectively. Subtraction yields an energy difference $\Delta
E=E_{r}-E_{l}$ for which all even order contributions have vanished,

$\displaystyle\Delta E$ $\displaystyle=$ $\displaystyle
2\partial_{\lambda}E_{0}+\frac{1}{3}\partial^{3}_{\lambda}E_{0}+\frac{1}{60}\partial^{5}_{\lambda}E_{0}+...$
(1) $\displaystyle=$ $\displaystyle\int d{\bf r}\Delta v({\bf r})\rho_{0}({\bf
r})+\frac{1}{6}\int d{\bf r}\Delta v({\bf
r})\partial^{2}_{\lambda}\rho_{0}({\bf r})+\frac{1}{120}\int d{\bf r}\Delta
v({\bf r})\partial^{4}_{\lambda}\rho_{0}({\bf r})+...$

where we are following Hellmann-Feynman’s theorem Feynman (1939); Wilson
(1962); von Lilienfeld (2009); von Rudorff and von Lilienfeld (2020),
$\partial_{\lambda}E_{0}=\int d{\bf r}\delta v({\bf r})\rho_{0}({\bf r})$,
$\partial^{3}_{\lambda}E_{0}=\int d{\bf r}\delta v({\bf
r})\partial^{2}_{\lambda}\rho_{0}({\bf r})$, $\partial^{5}_{\lambda}E_{0}=\int
d{\bf r}\delta v({\bf r})\partial^{4}_{\lambda}\rho_{0}({\bf r})$, and where
$\Delta v=2\delta v$. As already discussed in the context of alchemical
perturbation density functional theory von Rudorff and von Lilienfeld (2020),
these derivatives have clear meaning summing up to the integral over the
product of the perturbing Hamiltonian with the Taylor expansion in perturbed
electron densities. Such expansions have been shown to rapidly converge for
alchemical changes involving variation in the nuclear charge distribution and
fixed geometries von Rudorff (2021). Note that Eq. 1 can easily be adapted to
also estimate relative energies for arbitrarily distant alchemical
diastereomers, simply by increasing $\Delta\lambda$ to any other natural
number as long as it is not larger than the smallest nuclear charge of a
transmutating atom in the reference Hamiltonian. Correspondingly, energy
differences between alchemical diastereomers will grow linearly in
$|\Delta\lambda|$ as long as they happen to be situated on the same $\lambda$
dimension in chemical space. Fig. 1 illustrates this point: The energy
difference between N2 and BF is roughly half the size of the energy difference
between CO and BeNe.

Figure 1: Electronic energy as a continuous function of linear coupling
parameter $\lambda$ for some iso-electronic charge neutral diatomics with 14
protons. N2 and BF correspond to alchemical diastereomers of CO, implying that
the leading order term of their energy difference $\Delta E$ corresponds to
twice the corresponding Hellmann-Feynman derivative evaluated for CO (See Eq.
1). The relative energy between more distant alchemical diastereomers CO and
BeNe is twice that amount.

Note how this expansion recovers the case of alchemical chirality von Rudorff
and von Lilienfeld (2021) (where first order terms vanish) whenever the
perturbing Hamiltonian and the reference system are chosen such that the
parity of $\rho_{0}$ and anti-parity of $\Delta v$ result in an overlap
integral that exactly averages out. Due to the reflection plane in the
reference Hamiltonian’s external potential and due to the energetic degeneracy
up to third order, the iso-electronic compounds corresponding to $\hat{H}_{l}$
and $\hat{H}_{r}$ were dubbed ‘alchemical enantiomers’. However, the
alchemical chirality symmetry condition with vanishing Hellmann-Feynman
derivatives is not met when transmutating atoms with differing chemical
environments. And thus, in the more general case the non-vanishing values for
the Hellmann-Feynman derivative at the reference system become the leading
order contributions to energy differences of arbitrary iso-electronic compound
pairs. Correspondingly, we dub the latter ‘alchemical diastereomers’. In other
words, alchemical diastereomers become enantiomers with approximate energy
degeneracy whenever the relevant transmutating atoms in their averaged
Hamiltonian posses the same chemical environment.

It is also noteworthy that the antisymmetry of the perturbation leads to the
alternating signs in the left hand side expansion which, in turn, results in
the cancellation of the even order terms in Eq. 1. Such alternations are
frequently exploited, cf. the time reversal symmetry within molecular dynamics
simulation when using the Verlet algorithm: Velocity and higher odd order time
derivatives cancel and based on previous and current position, force based
classical Newtonian propagation is exact up to 4th order Frenkel and Smit
(2002). Accordingly if $E_{l}$ was also known, one could estimate $E_{r}$ as
follows

$\displaystyle E_{r}$ $\displaystyle=$ $\displaystyle
2E_{0}-E_{l}+\partial^{2}_{\lambda}E_{0}+\frac{1}{12}\partial^{4}_{\lambda}E_{0}+...$
(2)

where all odd order terms have vanished — in complete analogy to Verlet’s
integration. Correspondingly, this approach might be beneficial when the
electron density perturbation is available, (for example from coupled
perturbed self-consistent field calculation Lesiuk  _et al._ (2012)), and
yield $E_{r}$ exactly up to fourth order.

Perturbing the electronic density for the averaged Hamiltonian (in complete
analogy to the energy above), we can see the connection to Levy’s estimation
of iso-electronic energy differences

$\displaystyle\Delta E$ $\displaystyle\approx$ $\displaystyle\int d{\bf
r}\,\Delta v({\bf r})\,\bar{\rho}({\bf r})$ (3)

that relies on the averaged electron densities,
$\bar{\rho}=(\rho_{r}+\rho_{l})/2$. Levy (1978b) Consider

$\displaystyle\rho_{+}$ $\displaystyle=$
$\displaystyle\rho_{0}+\partial_{\lambda}\rho_{0}+\frac{1}{2}\partial^{2}_{\lambda}\rho_{0}+\frac{1}{6}\partial^{3}_{\lambda}\rho_{0}+\frac{1}{24}\partial^{4}_{\lambda}\rho_{0}...$
$\displaystyle\rho_{-}$ $\displaystyle=$
$\displaystyle\rho_{0}-\partial_{\lambda}\rho_{0}+\frac{1}{2}\partial^{2}_{\lambda}\rho_{0}-\frac{1}{6}\partial^{3}_{\lambda}\rho_{0}+\frac{1}{24}\partial^{4}_{\lambda}\rho_{0}...$
(4)

Averaging, and insertion in Levy’s formula, yields

$\displaystyle\int d{\bf r}\,\Delta v({\bf r})\,\bar{\rho}({\bf r})$
$\displaystyle=$ $\displaystyle\int d{\bf r}\Delta v({\bf
r})\left(\rho_{0}+\frac{1}{2}\partial^{2}_{\lambda}\rho_{0}+\frac{1}{24}\partial^{4}_{\lambda}\rho_{0}+...\right)$

for which the even order terms have also vanished. Term-wise comparison to Eq.
1 indicates that while Levy’s formula recovers the first order energy term,
i.e. Eq. 6, exactly, it overestimates the third, fifth, and seventh order
energy terms respectively by factors 3, 5, and 7, etc.etc. This implies that
Levy’s approximation should be less accurate than the first order term in Eq.
1.

##  III Results and discussions

Figure 2: Estimated total energy difference ($2\partial_{\lambda}E$) vs. SCF
based total energy difference ($\Delta E$) for 35 anti-symmetric BN doped
alchemical diastereomers of toluene (mean absolute error = 4.3 mHa). Figure 3:
Three randomly selected dimensions in chemical space defining three pairs of
alchemical diasteomers of toluene (out of the 35 dimensions for which relative
energies are shown in Fig. 2). Alchemical diastereomers are placed opposite
from each other in order to underscore their mutual relation in terms of
antisymmetric alchemical perturbations of toluene as the joint averaged
reference Hamiltonian, and with the lines indicating the corresponding
dimension in chemical compound space.

Truncating the expansion Eq. 1 after the leading order term, we have
numerically evaluated the predictive power when estimating the energy
difference as

$\displaystyle\Delta E$ $\displaystyle\approx$ $\displaystyle
2\partial_{\lambda}E_{0}\;=\;\int d{\bf r}\Delta v({\bf r})\rho_{0}({\bf r})$
(6)

Fig. 1 provides a graphical illustration of this idea: Due to the antisymmetry
condition of the perturbation, the slope at the reference system provides a
first order estimate of the energy difference between the diastereomers. We
note that if the energy was parabolic in $\lambda$, this estimate would become
exact.

In Tab 1, numerical predictions of Eq. 6 and of Levy’s formula are shown for
diatomics. Fair agreement, albeit not chemically accurate, is found despite
some of the perturbations being very substantial. As expected from above
analysis, Eq. 6 provides a systematically better estimate of the actual energy
difference than Levy’s formula. This is encouraging since it requires the SCF
procedure to obtain the electron density only for one (averaged reference)
system, rather than for each of the two end-points as it is necessary for
Levy.

If sufficiently accurate, the usefulness of this approximation could be
considerable: It would imply that given the electron density of a reference
systems, energy differences resulting from any of its combinatorially many
possible antisymmetric alchemical perturbations could be estimated with
negligible overhead merely by evaluating the corresponding Hellmann-Feynman
derivatives. In other words, relative energy estimates for an arbitrarily
large set of alchemical diastereomers, defined by their respective
$\\{\lambda_{i}\\}$-perturbations, can be generated for negligible additional
computational cost—as long as they all share the same averaged reference
Hamiltonian. We have numerically exemplified this point for the BN doping of
the aromatic moiety in the molecule toluene. In this case, BN doping defines
35 $\lambda$-dimensions space, along which we have estimated all the
respective energy differences between the corresponding alchemical
diastereomers — for ‘free’, i.e. via Eq. 6 and only based on the electron
density obtained for the joint averaged reference Hamiltonian of toluene.
Figs. 2 shows a scatter plot of the numerical estimates of energy differences
between 35 pairs of alchemical diastereomers of toluene. The resulting MAE is
4.3 mHa which is close to chemical accuracy ($\sim$1 kcal/mol) and similar to
the accuracy of hybrid density functional approximations. Fig. 3 exemplifies
the central role of the averaged reference Hamiltonian of toluene for three
pairs of alchemical diastereomers.

$\hat{H}_{l}$-$\hat{H}_{0}$-$\hat{H}_{r}$ |  $\Delta E$[Ha] |  $2\partial_{\lambda}E_{0}$ [Ha] | Levy [Ha]Levy (1979)
---|---|---|---
NN-CO-BF | 13.176560 | 13.196927 | 13.123976
NN-BF-LiNa | 52.201112 | 52.495903 | 51.669013
NN-BeNe-HAl | 115.963855 | 117.201259 | 112.968808
CO-BF-BeNe | 26.207973 | 26.247951 | 26.132007
CO-BeNe-HeMg | 77.811426 | 78.134173 | 77.094902
CO-LiNa-Si | 152.885958 | 155.007039 | 148.255308
BF-BeNe-LiNa | 39.024551 | 39.067086 | 38.958482
BF-LiNa-HAl | 102.787295 | 103.338026 | 101.560490
BeNe-LiNa-HeMg | 51.603454 | 51.669013 | 51.482530
BeNe-HeMg-Si | 126.677985 | 127.795948 | 124.707031
LiNa-HeMg-HAl | 63.762743 | 63.897974 | 63.490776
HeMg-HAl-Si | 75.074531 | 75.312538 | 74.768959
Table 1:  Relative electronic potential energy estimates between first and
third diatomics. $2\partial_{\lambda}E_{0}$ corresponds to Eq. 6, whereas Levy
corresponds to Eq. 3 Levy (1979). Numbers obtained at 1.1 Å interatomic
distance using the PBE0 functional in the pcX-2 basis.

##  IV Computational Details

All numerical calculations were done using PySCF Sun  _et al._ (2017); Harris
_et al._ (2020); Pérez and Granger (2007) with the PBE0 density functional
approximation Perdew  _et al._ (1996); Ernzerhof and Scuseria (1999); Adamo
and Barone (1999). The importance of basis set effects on alchemical
interpolations having been established previously Domenichini  _et al._
(2020); Domenichini and von Lilienfeld (2022), we have used Jensen’s pc2 basis
set for hydrogens Jensen (2001, 2002), and the universal pcX-2 basis by
Ambroise and Jensen for all other atoms Ambroise and Jensen (2019). The
interatomic distance of all systems in the diatomics series was set to 1.1Å.
The geometry of the 70 BN doped toluene derivatives was kept fixed to the
equilibrium geometry of toluene obtained at the same level of theory and using
Hermann’s geometry optimizer PyBerny Hermann (2020); it is shown in Tab. 2.

Atom | X [Å] | Y [Å] | Z [Å]
---|---|---|---
H | -0.00048 | 2.76025 | 1.05495
H | 0.88360 | 2.79792 | -0.46619
H | -0.88315 | 2.79793 | -0.46700
C | -0.00000 | 2.38982 | 0.02594
C | 0.00001 | 0.89122 | -0.00419
C | 1.19407 | 0.17608 | -0.00721
C | -1.19407 | 0.17608 | -0.00705
C | 1.19708 | -1.21024 | -0.00623
C | -0.00000 | -1.90984 | -0.00527
C | -1.19708 | -1.21023 | -0.00606
H | 2.13561 | 0.71445 | -0.01262
H | 2.13883 | -1.74590 | -0.01129
H | -0.00000 | -2.99283 | -0.00868
H | -2.13883 | -1.74589 | -0.01100
H | -2.13560 | 0.71446 | -0.01235
Table 2:  Optimized atomic coordinates of toluene.

##  V Conclusion

We have generalized the alchemical chirality concept towards the notion of
alchemical antisymmetric perturbations that result in vanishing even order
contributions to the relative energies of alchemical diastereomers. Our
analysis also provides an interpretation of Levy’s formula Levy (1979).
Numerical evidence suggests that the leading first order term gives
meaningful, and sometimes even accurate, estimates of energy differences
between alchemical diastereomers that are close to each other. The choice of
the averaged reference system is key, it defines the number of dimensions in
chemical space along which energy differences between alchemical diastereomers
can be estimated with negligible computational overhead. We have exemplified
this point for 35 BN doped alchemical diastereomers of toluene for which
energy differences were calculated using 35 Hellmann-Feynman derivates with
the same single electron density and yielding a mean absolute error of only
4.3 mHa. Future studies will deal with the role of the quality of the electron
density used within the Hellmann-Feynman derivative, with the correlation
between prediction error and magnitude of perturbation, and the importance of
higher order terms.

###### Acknowledgements.

The authors acknowledge discussions with M. Chaudhari, D. Khan, S. Krug, M.
Meuwly, GF von Rudorff, M. Sahre, and A. Savin. O.A.v.L. has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement No. 772834). O.A.v.L.
has received support as the Ed Clark Chair of Advanced Materials and as a
Canada CIFAR AI Chair.

## References

  * Kirkpatrick and Ellis (2004) P. Kirkpatrick and C. Ellis, Nature 432, 823 (2004).
  * von Rudorff and von Lilienfeld (2020) G. F. von Rudorff and O. A. von Lilienfeld, Phys. Rev. Research 2, 023220 (2020).
  * Foldy (1951) L. L. Foldy, Phys. Rev. 83, 397 (1951).
  * E. B. Wilson, Jr. (1962) E. B. Wilson, Jr., J. Chem. Phys. 36, 2232 (1962).
  * Levy (1978a) M. Levy, The Journal of Chemical Physics 68, 5298 (1978a).
  * Levy (1979) M. Levy, The Journal of Chemical Physics 70, 1573 (1979).
  * Politzer and Parr (1974) P. Politzer and R. G. Parr, J. Chem. Phys. 61, 4258 (1974).
  * Weigend  _et al._ (2004) F. Weigend, C. Schrodt, and R. Ahlrichs, J. Chem. Phys. 121, 10380 (2004).
  * von Lilienfeld  _et al._ (2005) O. A. von Lilienfeld, R. Lins, and U. Rothlisberger, Phys. Rev. Lett. 95, 153002 (2005).
  * Wang  _et al._ (2006) M. Wang, X. Hu, D. N. Beratan, and W. Yang, J. Am. Chem. Soc. 128, 3228 (2006).
  * Beste  _et al._ (2006) A. Beste, R. J. Harrison, and T. Yanai, J. Phys. Chem. 125, 074101 (2006).
  * Sheppard  _et al._ (2010) D. Sheppard, G. Henkelman, and O. A. von Lilienfeld, J. Chem. Phys. 133, 084104 (2010).
  * Lesiuk  _et al._ (2012) M. Lesiuk, R. Balawender, and J. Zachara, J. Chem. Phys. 136, 034104 (2012).
  * Balawender  _et al._ (2013) R. Balawender, M. A. Welearegay, M. Lesiuk, F. De Proft, and P. Geerlings, J. Chem. Theory Comput. 9, 5327 (2013).
  * Weigend (2014) F. Weigend, The Journal of chemical physics 141, 134103 (2014).
  * Chang  _et al._ (2016) K. Y. S. Chang, S. Fias, R. Ramakrishnan, and O. A. von Lilienfeld, J. Chem. Phys. 144, 174110 (2016), http://dx.doi.org/10.1063/1.4947217 .
  * Fias  _et al._ (2017) S. Fias, F. Heidar-Zadeh, P. Geerlings, and P. W. Ayers, Proceedings of the National Academy of Sciences 114, 11633 (2017).
  * Munoz and Cardenas (2017) M. Munoz and C. Cardenas, Physical Chemistry Chemical Physics 19, 16003 (2017).
  * Fias  _et al._ (2018) S. Fias, K. S. Chang, and O. A. von Lilienfeld, The journal of physical chemistry letters 10, 30 (2018).
  * Saravanan  _et al._ (2017) K. Saravanan, J. R. Kitchin, O. A. von Lilienfeld, and J. A. Keith, The Journal of Physical Chemistry Letters 8, 5002 (2017).
  * Al-Hamdani  _et al._ (2017) Y. S. Al-Hamdani, A. Michaelides, and O. A. von Lilienfeld, The Journal of chemical physics 147, 164113 (2017).
  * Balawender  _et al._ (2018) R. Balawender, M. Lesiuk, F. De Proft, and P. Geerlings, Journal of chemical theory and computation 14, 1154 (2018).
  * Chang and von Lilienfeld (2018) K. S. Chang and O. A. von Lilienfeld, Physical Review Materials 2, 073802 (2018).
  * Griego  _et al._ (2018) C. D. Griego, K. Saravanan, and J. A. Keith, Advanced Theory and Simulations , 1800142 (2018).
  * von Rudorff and von Lilienfeld (2019) G. F. von Rudorff and O. A. von Lilienfeld, The Journal of Physical Chemistry B 123, 10073 (2019).
  * Balawender  _et al._ (2019) R. Balawender, M. Lesiuk, F. De Proft, C. Van Alsenoy, and P. Geerlings, Physical Chemistry Chemical Physics 21, 23865 (2019).
  * Shiraogawa and Ehara (2020) T. Shiraogawa and M. Ehara, The Journal of Physical Chemistry C 124, 13329 (2020).
  * Griego  _et al._ (2020a) C. D. Griego, J. R. Kitchin, and J. A. Keith, International Journal of Quantum Chemistry , e26380 (2020a).
  * Griego  _et al._ (2020b) C. D. Griego, L. Zhao, K. Saravanan, and J. A. Keith, AIChE Journal 66, e17041 (2020b).
  * von RAtomicAPDFTudorff and von Lilienfeld (2020) G. F. von RAtomicAPDFTudorff and O. A. von Lilienfeld, Physical Chemistry Chemical Physics 22, 10519 (2020).
  * Muñoz  _et al._ (2020) M. Muñoz, A. Robles-Navarro, P. Fuentealba, and C. Cárdenas, The Journal of Physical Chemistry A 124, 3754 (2020).
  * Griego  _et al._ (2021) C. D. Griego, J. R. Kitchin, and J. A. Keith, International Journal of Quantum Chemistry 121, e26380 (2021).
  * Gómez  _et al._ (2021) T. Gómez, P. Fuentealba, A. Robles-Navarro, and C. Cárdenas, Journal of Computational Chemistry 42, 1681 (2021).
  * Eikey  _et al._ (2022a) E. A. Eikey, A. M. Maldonado, C. D. Griego, G. F. Von Rudorff, and J. A. Keith, The Journal of chemical physics 156, 064106 (2022a).
  * Shiraogawa and Hasegawa (2022) T. Shiraogawa and J.-y. Hasegawa, The Journal of Physical Chemistry Letters 13, 8620 (2022).
  * Shiraogawa  _et al._ (2022) T. Shiraogawa, G. Dall’Osto, R. Cammi, M. Ehara, and S. Corni, Physical Chemistry Chemical Physics 24, 22768 (2022).
  * Eikey  _et al._ (2022b) E. A. Eikey, A. M. Maldonado, C. D. Griego, G. F. Von Rudorff, and J. A. Keith, The Journal of Chemical Physics 156, 204111 (2022b).
  * Shiraogawa and Hasegawa (2023) T. Shiraogawa and J.-y. Hasegawa, (2023).
  * Balawender and Geerlings (2023) R. Balawender and P. Geerlings, in  _Chemical Reactivity_ (Elsevier, 2023) pp. 15–57.
  * von Rudorff and von Lilienfeld (2021) G. F. von Rudorff and O. A. von Lilienfeld, Science Advances 7 (2021), https://doi.org/10.1126/sciadv.abf1173.
  * Feynman (1939) R. P. Feynman, Phys. Rev. 56, 340 (1939).
  * Wilson (1962) E. B. Wilson, J. Chem. Phys. 36, 2232 (1962).
  * von Lilienfeld (2009) O. A. von Lilienfeld, J. Chem. Phys. 131, 164102 (2009).
  * von Rudorff (2021) G. F. von Rudorff, The Journal of Chemical Physics 155, 224103 (2021).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ (Academic Press, 2002).
  * Levy (1978b) M. Levy, J. Chem. Phys. 68, 5298 (1978b).
  * Sun  _et al._ (2017) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K. Chan, “Pyscf: the python‐based simulations of chemistry framework,”  (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1340 .
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Pérez and Granger (2007) F. Pérez and B. E. Granger, Computing in Science and Engineering 9, 21 (2007).
  * Perdew  _et al._ (1996) J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996).
  * Ernzerhof and Scuseria (1999) M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999).
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Domenichini  _et al._ (2020) G. Domenichini, G. F. von Rudorff, and O. A. von Lilienfeld, The Journal of Chemical Physics 153, 144118 (2020).
  * Domenichini and von Lilienfeld (2022) G. Domenichini and O. A. von Lilienfeld, The Journal of Chemical Physics 156, 184801 (2022).
  * Jensen (2001) F. Jensen, J. Chem. Phys. , 9113 (2001).
  * Jensen (2002) F. Jensen, J. Chem. Phys. , 7372 (2002).
  * Ambroise and Jensen (2019) M. A. Ambroise and F. Jensen, Journal of Chemical Theory and Computation 15, 325 (2019), https://doi.org/10.1021/acs.jctc.8b01071 .
  * Hermann (2020) J. Hermann, “Pyberny is an optimizer of molecular geometries with respect to the total energy, using nuclear gradient information.” (accessed in November 2020), Github project: https://github.com/jhrmnn/pyberny, Zenodo database: https://doi.org/10.5281/zenodo.3695038 .