#  Revealing quantum contextuality using a single measurement device by
generalizing measurement non-contextuality

Jaskaran Singh<EMAIL_ADDRESS>Department of Physical Sciences,
Indian Institute of Science Education and Research (IISER) Mohali, Sector 81
SAS Nagar, Manauli PO 140306 Punjab India.  Arvind<EMAIL_ADDRESS>Department of Physical Sciences, Indian Institute of Science Education and
Research (IISER) Mohali, Sector 81 SAS Nagar, Manauli PO 140306 Punjab India.

###### Abstract

In this letter we generalize Spekkens’ notion of measurement non-contextuality
(NC). We show that any non-contextual ontological model based on this notion
of contextuality fails to explain the statistics of outcomes of a single
carefully constructed POVM executed sequentially on a quantum system. The
context of measurement arises from the different configurations in which the
device can be used. We develop an inequality from the non-contextual ontic
model, and construct corresponding quantum situations where the measurement
outcomes from the device violate this NC inequality. Our work brings out the
hitherto unexplored implications of contextuality for a single measurement
device, and paves the way for further study of consequences of contextuality
for sequential measurements.

Introduction.— Quantum theory is contextual since the outcomes of measurements
depend on the context of the measurement, namely, the set of commuting
observables being measured along with the desired observable Simon Kochen
(1968); Liang  _et al._ (2011). Unlike quantumness of composite systems
revealed via Bell type inequalities BELL (1966); Brunner  _et al._ (2014);
Clauser  _et al._ (1969), quantum contextuality can be demonstrated on single
indivisible systems and the simplest such scenario involves a three
dimensional quantum system and five different projective measurements Klyachko
_et al._ (2008). While the first proof that quantum theory is contextual was
provided by Kochen and Specker Kochen and Specker (1975), involving $117$
different projectors in a three dimensional Hilbert space, over time a number
of simpler and more systematic ways of revealing quantum contextuality,
particularly the ones based on graph theory, have become available Cabello
_et al._ (2014); Acín  _et al._ (2015). The graph theoretic approach has been
successful in identifying new contextual scenarios Araújo  _et al._ (2013);
Amaral  _et al._ (2014); Kurzyński and Kaszlikowski (2012); Lisoněk  _et al._
(2014); Cabello  _et al._ (2015); Kleinmann  _et al._ (2012); Cabello  _et
al._ (2016); Yu and Oh (2012); Cabello (2013); Bharti  _et al._ (2019a),
simplifying formulations of contextuality monogamy Ramanathan  _et al._
(2012), contextuality non-locality relationship Kurzyński  _et al._ (2014),
developing robust self tests Bharti  _et al._ (2019b) and information
theoretic applications of contextuality Singh  _et al._ (2017); Cabello  _et
al._ (2011); Troupe and Farinholt (2015); Saha  _et al._ (2019).

While the majority analysis in quantum contextuality has involved projective
measurements, the approach developed by Spekkens generalizes the same to
positive-operator-valued measures (POVMs), as well as provides a notion of
contextuality for preparations and transformations Spekkens (2005). With the
use of this generalized approach it was possible to bring down the number of
measurements required to exhibit contextual behavior to just $3$ on a single
qubit Kunjwal and Ghosh (2014). Furthermore, this generalization provides a
technique for noise-robust experimental demonstration of contextuality Mazurek
_et al._ (2016); Kunjwal and Spekkens (2018, 2015); Xu  _et al._ (2016) and
leads to information theoretic applications of quantum situations involving
preparation contextuality Spekkens  _et al._ (2009); Chailloux  _et al._
(2016); Ambainis  _et al._ (2019). While the minimum number of measurements
required to reveal quantum contextuality so far is three Kunjwal and Ghosh
(2014), no physical principle prohibits a smaller number. Could a more
generalized approach involving analysis based on sequential applications of
POVMs be used to reveal contextuality of quantum situations? Is it possible to
generalize Spekkens assumptions to sharpen the contextual-noncontextual
divide? Although such generalizations have been attempted Chaturvedi and Saha
(2020); Schmid  _et al._ (2020), these questions are still not fully resolved.
We take up these questions and show how such generalizations lead to the
unearthing of quantum contextuality in more general situations.

In this work we propose a notion of measurement non-contextuality (NC) which
is a generalization of, and is motivated by the assumption of measurement NC
developed by Spekkens Spekkens (2005). Our generalized notion opens up several
new and interesting scenarios which exhibit quantum correlations that cannot
be simulated using non-contextual ontological models. We formulate a new
experimentally verifiable NC inequality based on the proposed notion of
measurement NC. The scenario involves the use of a single measurement device
executed twice sequentially to reveal the contextuality of the quantum
situation involved. We show that this inequality is violated by quantum
situations, which implies that the underlying non-contextual ontic models
cannot reproduce the statistics of even a single measurement device.

Measurement non-contextuality.— The definition of measurement NC as formulated
by Spekkens Spekkens (2005), relies on defining a notion of equivalence among
different experimental procedures. Specifically, two measurement procedures
$\mathcal{M}$ and $\mathcal{M}^{\prime}$ are deemed equivalent if they yield
the same statistics for every possible preparation procedure P, that is

$p(k|\mathcal{M},\text{P})=p(k|\mathcal{M}^{\prime},\text{P})\implies\mathcal{M}\equiv\mathcal{M}^{\prime}.$
(1)

Where $p(k|\mathcal{M},\text{P})$ is the probability of obtaining the outcome
$k$ in the measurement $\mathcal{M}$ for preparation P. Consider an ontic
model where ontic state $\lambda$ can be used to predict the measurement
outcomes. Based on the definition of equivalence classes for measurement
procedures, and as motivated by Leibniz’s principle of indiscernibles, the
definition of measurement NC is

$p(k|\mathcal{M},\text{P})=p(k|\mathcal{M}^{\prime},\text{P})\implies\xi_{\mathcal{M},k}(\lambda)=\xi_{\mathcal{M}^{\prime},k}(\lambda),$
(2)

where $\xi_{\mathcal{M},k}(\lambda)$ is the ontic response function for
assigning the probability of obtaining an outcome $k$ in a measurement
$\mathcal{M}$ given the ontic state $\lambda$ with

$\xi_{\mathcal{M},k}(\lambda)\leq
1,\quad\sum_{k}\xi_{\mathcal{M},k}(\lambda)=1.$ (3)

Thus, the assumption of measurement NC implies that the ontic response
function for different outcomes of equivalent measurement procedures is the
same.

We consider a modified assumption of measurement NC which leads to several
interesting contextual scenarios not possible under the NC assumption as
conceived by Spekkens. Firstly, we define two measurement procedures,
$\mathcal{M}$ and $\mathcal{M^{\prime}}$ to be similar if

$p(k|\mathcal{M},P)=s.p(k|\mathcal{M^{\prime}},P)\quad\forall\,P,$ (4)

where $s$ is drawn from a known probability distribution independent of
$\mathcal{M}$, $\mathcal{M^{\prime}}$ and $P$. In this case, both the
measurement procedures can be clearly distinguished. Based on the notion of
similar measurement procedures, we define our assumption of measurement NC as,

$p(k|\mathcal{M},P)=s.p(k|\mathcal{M^{\prime}},P)\implies\xi_{\mathcal{M},k}(\lambda)=s.\xi_{\mathcal{M}^{\prime},k}(\lambda).$
(5)

Which states that an ontic model is non-contextual if the ontic response
functions corresponding to different outcomes for similar measurement
procedures must be the same upto the factor $s$ involved in defining the
similarity of measurement procedures. This wider notion of measurement NC
ontic models allows us to sharpen the dichotomy between contextual and non-
contextual situations and helps us in declaring more situations as contextual.
Our notion of NC can be thought of as inducing an ontic representation on
similar measurement procedures due to the fact that the parameter $s$ is
known. We utilize this assumption to describe a situation where measurement NC
can be revealed using a single measurement device, albeit applied twice
sequentially.

Example— As an example of our generalization consider two different projective
measurements on a single spin-$1/2$ system given as

$\displaystyle\mathcal{M}_{Z}:=\\{|0\rangle\langle 0|,|1\rangle\langle 1|\\},$
(6)
$\displaystyle\mathcal{M}_{X}:=\\{|+\rangle\langle+|,|-\rangle\langle-|\\},$

where $|0\rangle$ and $|1\rangle$ are the eigenkets of Pauli spin-$Z$, and
$|+\rangle$ and $|-\rangle$ are the eigenkets of Pauli spin-$X$ operator. We
define another POVM $\mathcal{M}_{ZX}$ as,

$\mathcal{M}_{ZX}:=\\{s|0\rangle\langle 0|,s|1\rangle\langle
1|,(1-s)|+\rangle\langle+|,(1-s)|-\rangle\langle-|\\},$ (7)

where $s\in\mathcal{S}$ is sampled from some known probability distribution
$\mathcal{S}$. The above measurement can be thought of as performing the
$\mathcal{M}_{Z}$ measurement with probability $s$ and $\mathcal{M}_{X}$ with
probability $(1-s)$. Spekkens’ assumption of measurement NC cannot be applied
here because the probability for the outcomes $\ket{i}\bra{i}$,
$\forall\,i\in\\{0,1,+,-\\}$ are‘ clearly different from those for the
outcomes $s\ket{i}\bra{i}$. Therefore, $\mathcal{M}_{Z}$ and
$\mathcal{M}_{ZX}$ cannot be the different ‘contexts’ of the outcomes
$\ket{0}\bra{0}$ (or for that matter for any other outcome).

However, our generalized assumption of measurement NC for the measurements
described above can be applied to this situation. Let $\xi_{i,k}$ denote the
ontic response function assigned to the outcome $k$ for measurement
$i\in\\{Z,X,ZX\\}$, where,

$\xi_{i,k}\geq 0,\quad\sum_{k}\xi_{i,k}=1.$ (8)

As an example let us consider the outcome $|0\rangle\langle 0|$ appearing in
the measurements $\mathcal{M}_{Z}$ and $\mathcal{M}_{ZX}$. We see that for all
possible preparations $P$,

$p(|0\rangle\langle 0||\mathcal{M}_{Z},P)=s.p(|0\rangle\langle
0||\mathcal{M}_{ZX},P)\quad\forall P,$ (9)

which implies according to our assumption of measurement NC as

$\xi_{Z,|0\rangle\langle 0|}=s.\xi_{ZX,|0\rangle\langle 0|}.$ (10)

A similar analysis holds for all the outcomes for $\mathcal{M}_{Z}$,
$\mathcal{M}_{X}$ and $\mathcal{M}_{ZX}$.

According to our definition even when outcomes are perfectly distinguishable,
a notion of context can still be defined. This can be achieved only when the
statistics of the outcomes are similar to each other according to Eq. (5).
This identifies a larger set of contexts for measurements as compared to what
was possible earlier. It should be noted that for $s=1$, similarity becomes
equivalence and our condition transforms into the condition set forth by
Spekkens.

The Setup— Consider a measurement procedure able to function in two different
configurations and we assume a non-contextual ontological model explaining the
statistics of both the configurations.

Figure 1: A schematic diagram of the two configurations of the measurement
device $\mathcal{M}$. In the first configuration, the device samples the
measurements $\mathcal{M}_{i}$ from a random number generator (RNG) with
probability $p_{i}$. In this configuration the dashed outcomes are collected
as final outcome $\mathcal{K}$, while outcomes with solid lines are collected
as final outcomes $E_{i}$. In the second configuration, each measurement
$\mathcal{M}_{i}$ can be performed independently with outcomes given by solid
and dashed lines.

The schematics of our setup is shown in Fig 1. The setup is comprised of $5$
projective measurements $\\{\mathcal{M}_{i}\\}$, with three outcomes each. We
choose these projective measurements to be the ones utilized in the KCBS
scenario. Using these measurements we can enable the device to work in two
possible configurations as detailed below:

  * C1:

In this configuration the measurements $\\{\mathcal{M}_{i}\\}$ are sampled
from a probability distribution $\\{p_{i}\\}$, $i\in\\{0,1,2,3,4\\}$, and the
device implements the projective measurement $\\{\Pi_{i},\Pi_{i\oplus
1},K_{i}\\}$, where $K_{i}$ is added to complete the measurement and $\oplus$
is addition modulo $5$. The projectors $\\{\Pi_{i}\\}$ satisfy ${\rm
Tr}(\Pi_{i}\Pi_{i\oplus 1})=0$ and are the same as the ones used in the
derivation of KCBS inequality. The resultant measurement is then a POVM
$\mathcal{M}$ with outcomes,

$\mathcal{M}:\\{E_{0},E_{1},E_{2},E_{3},E_{4},\mathcal{K}\\},$ (11)

where $E_{i}=(p_{i}+p_{i\oplus 1})\Pi_{i}$ and
$\mathcal{K}=2\sum_{i}p_{i}K_{i}$ with

$\sum_{0}^{4}E_{i}=\sum_{i=0}^{4}(p_{i}+p_{i\oplus 1})\Pi_{i}\leq\mathds{1},$
(12)

which is a consequence of completeness of measurement.

  * C2:

In this configuration we choose a particular setting $i$, implementing a
projective measurement on our device

$\mathcal{M}_{i}:\\{\Pi_{i},\Pi_{i\oplus 1},K_{i}\\}.$ (13)

This is akin to blocking all measurement outcomes $j\neq i$. From completeness
we again have

$\Pi_{i}+\Pi_{i\oplus 1}\leq\mathds{1}\quad\forall i.$ (14)

Five such settings are possible corresponding to the projective measurements
detailed above and each is labeled by a value of $i$.

We now formulate an ontological description of the aforementioned measurement
device. For each measurement outcome $\Pi_{i}$ from the measurement
$\mathcal{M}_{i}$, there corresponds a response function $\xi_{i}(\lambda)$,
where $\lambda$s are the ontic variables. Therefore an ontological description
of $\mathcal{M}_{i}$ will have the form

$\mathcal{M}_{i,\lambda}:\\{\xi_{i}(\lambda),\xi_{i\oplus
1}(\lambda),\xi_{K_{i}}(\lambda)\\}.$ (15)

This ontological model of configuration C2 induces an ontological description
for the configuration C1. This is due to the fact that the projectors involved
in both are the same. From the generalized assumption of measurement NC, the
induced ontological model of the measurement device in C1 configuration will
have the following description,

$\displaystyle\mathcal{M}_{\lambda}:\\{\xi_{0}^{\prime}(\lambda),\xi_{1}^{\prime}(\lambda),\xi_{2}^{\prime}(\lambda),\xi_{3}^{\prime}(\lambda),\xi_{4}^{\prime}(\lambda),\xi_{\mathcal{K}}(\lambda)\\},$
$\displaystyle\xi_{i}^{\prime}(\lambda)=(p_{i}+p_{i\oplus
1})\xi_{i}(\lambda),$ (16)

such that
$\\{\xi_{i}(\lambda),\xi_{i}^{\prime}(\lambda),\xi_{\mathcal{K}}(\lambda)\\}\in\left[0,1\right]$,
which fits quite well with our generalized notion of measurement NC given in
Eq. (2)

Further, we have

$\sum_{i=0}^{4}(p_{i}+p_{i\oplus
1})\xi_{i}(\lambda)=\sum_{i=0}^{4}\xi_{i}^{\prime}(\lambda)\leq 1,$ (17)

which is a consequence of Eq. (12). Since, by construction, the projectors
$\Pi_{i}$ and $\Pi_{i\oplus 1}$ are orthogonal, we have,

$\int_{\lambda\in\Lambda}\xi_{i}(\lambda)\xi_{i\oplus 1}(\lambda)d\lambda=0,$
(18)

which states that the overlap between the two response functions
$\xi_{i}(\lambda)$ and $\xi_{i\oplus 1}(\lambda)$ is zero.

Contextuality from a single measurement device.— We are now ready to propose
an inequality to be tested by sequential measurements and explicitly derive
its maximum non-contextual value. We then show that a quantum description of
the same leads to a violation.

######  Proposition 1.

The sum of probabilities for obtaining outcome $E_{i}$ given the outcome
$E_{j}$ when measurement $\mathcal{M}$ is performed twice sequentially for
$p_{i}=\frac{1}{5}$ is bounded by a non-contextual ontological model as,

$\mathcal{C}=\sum_{i,j=0}^{4}p(E_{i}|E_{j})\overset{\text{NC}}{\leq}3.20,$
(19)

where $p(E_{i}|E_{j})$ denotes the conditional probability of obtaining
outcome $E_{i}$ given $E_{j}$ when measurement $\mathcal{M}$ is performed
sequentially and NC represents the maximum non-contextual bound.

###### Proof.

In a non-contextual ontological model the probability of obtaining sequential
outcome $E_{i}$ given $E_{j}$ represented by response functions
$\xi_{i}(\lambda)$ and $\xi_{j}(\lambda)$ is given as

$p(E_{i}|E_{j})=2\int_{\lambda\in\Lambda}p_{i}\xi_{i}(\lambda)\xi_{j}(\lambda)d\lambda,$
(20)

where the ontic state after first measurement is represented by the response
function $\xi_{j}(\lambda)$. Using Eqn. (18) and (20) we have,

$\displaystyle\mathcal{C}$
$\displaystyle=\int_{\lambda\in\Lambda}\left[2\sum_{i=0}^{4}\left(p_{i}\xi_{i}^{2}(\lambda)+2\left(p_{i}\xi_{i}(\lambda)\xi_{i\oplus
2}(\lambda)\right)\right)\right]d\lambda$ (21)
$\displaystyle\leq\int_{\lambda\in\Lambda}\left[\frac{2}{5}\sum_{i=0}^{4}\left(\xi_{i}(\lambda)+2\left(\xi_{i}(\lambda)\xi_{i\oplus
2}(\lambda)\right)\right)\right]d\lambda$ $\displaystyle\leq
2+\frac{4}{5}\int_{\lambda\in\Lambda}\left[\left(\begin{aligned}
&\xi_{0}(\lambda)\left[\xi_{2}(\lambda)+\xi_{3}(\lambda)\right]\\\
&+\xi_{1}(\lambda)\left[\xi_{3}(\lambda)+\xi_{4}(\lambda)\right]\\\
&+\xi_{2}(\lambda)\xi_{4}(\lambda)\end{aligned}\right)\right]d\lambda$
$\displaystyle\leq
2+\frac{4}{5}\int_{\lambda\in\Lambda}\left(\xi_{0}(\lambda)+\xi_{1}(\lambda)+\xi_{2}(\lambda)\xi_{4}(\lambda)\right)d\lambda$
$\displaystyle=2+\frac{4}{5}\int_{\lambda\in\Lambda}\left(\xi_{i}(\lambda)+\xi_{i\oplus
1}(\lambda)+\xi_{i\oplus 2}(\lambda)\xi_{i\oplus
3}(\lambda)\right)d\lambda\quad\forall i$
$\displaystyle=2+\frac{4}{5}\left(1+\int_{\lambda\in\Lambda}\xi_{i}(\lambda)\xi_{i\oplus
2}(\lambda)d\lambda\right)$ $\displaystyle\leq 3.20,$

where for the third inequality we have used $\xi_{i}(\lambda)+\xi_{i\oplus
1}(\lambda)=1$ and for last inequality we have used the fact that overlap
between two non-orthogonal states as appearing in the scenario cannot be more
than $\frac{1}{2}$. This can be motivated as follows: The overlap cannot be
$0$ as it would imply the response functions are orthogonal and can be
discriminated in a single shot experiment. If the overlap is $1$, then
$\xi_{i}(\lambda)=\xi_{i\oplus 2}(\lambda)$, which further implies
$\xi_{i}(\lambda)$ is orthogonal to $\xi_{i\oplus 3}(\lambda)$, which again is
not allowed. Furthermore, due to symmetry of the problem $\xi_{i}(\lambda)$
must have an equal overlap with $\xi_{i\oplus 2}(\lambda)$ and $\xi_{i\oplus
3}(\lambda)$ $\forall i$. If the overlap of $\xi_{i}(\lambda)$ with
$\xi_{i\oplus 2}(\lambda)$ is greater than $\frac{1}{2}$, then the overlap
between $\xi_{i}(\lambda)$ with $\xi_{i\oplus 3}(\lambda)$ has to be less than
$\frac{1}{2}$ in order to avoid any overlap between $\xi_{i\oplus 2}(\lambda)$
and $\xi_{i\oplus 3}(\lambda)$, which is by definition zero. The only allowed
possibility which maximizes the quantity of interest is then for the overlap
to be $\frac{1}{2}$. This concludes the proof. ∎

In order to show that the aforementioned non-contextual ontological model is
unable to describe quantum statistics, it is enough to provide a single
counter-example which violates the inequality (19). We explicitly provide the
projectors $\Pi_{i}=|v_{i}\rangle\langle v_{i}|$ which will maximally violate
the inequality (19) as follows:

$\displaystyle\langle v_{0}|$
$\displaystyle=\left(1,0,\sqrt{\cos\frac{\pi}{5}}\right)$ (22)
$\displaystyle\langle v_{1}|$
$\displaystyle=\left(\cos\frac{4\pi}{5},-\sin\frac{4\pi}{5},\sqrt{\cos\frac{\pi}{5}}\right),$
$\displaystyle\langle v_{2}|$
$\displaystyle=\left(\cos\frac{2\pi}{5},\sin\frac{2\pi}{5},\sqrt{\cos\frac{\pi}{5}}\right),$
$\displaystyle\langle v_{3}|$
$\displaystyle=\left(\cos\frac{2\pi}{5},-\sin\frac{2\pi}{5},\sqrt{\cos\frac{\pi}{5}}\right),$
$\displaystyle\langle v_{4}|$
$\displaystyle=\left(\cos\frac{4\pi}{5},\sin\frac{4\pi}{5},\sqrt{\cos\frac{\pi}{5}}\right).$

The inequality $\mathcal{C}$ then takes on the value,

$\displaystyle\mathcal{C}$ $\displaystyle=2\sum_{i=0}^{4}\left(p_{i}|\langle
v_{i}|v_{i}\rangle|^{2}+2p_{i}|\langle v_{i}|v_{i\oplus 2}\rangle|^{2}\right)$
(23) $\displaystyle=2(1+3-\sqrt{5})$ $\displaystyle=3.52,$

which is greater than the predicted value for an ontologically non-contextual
model. Therefore, such models are unable to predict statistics even for single
measurement devices.

Conclusion.— In this letter we considered a generalization of the notion of
measurement NC proposed by Spekkens. Such a generalization opens up several
new interesting scenarios which exhibit a quantum superiority over classical
correlations by violating the non-contextual bound on the ontic model. We have
detailed one such scenario in which the statistics of outcomes from a single
measurement device is not reproducible from an underlying non-contextual ontic
model. Such a scenario has been constructed for the first time to the best of
our knowledge. Our results pave the way for future theoretical as well as
experimental work to unearth the contextuality of new quantum situations. A
straightforward extension of our work would entail generalization of the
single measurement device scenario to incorporate multiple sequential
measurements of the same or different non-contextual scenarios involving
$n$-outcomes and will be taken up elsewhere.

Acknowledgements.— JS acknowledges UGC India for financial support. JS and A
acknowledge the financial support from DST/ICPS/QuST/Theme-1/2019/General
Project number Q-68.

## References

  * Simon Kochen (1968) E. S. Simon Kochen, Indiana Univ. Math. J. 17, 59 (1968).
  * Liang  _et al._ (2011) Y.-C. Liang, R. W. Spekkens, and H. M. Wiseman, Physics Reports 506, 1 (2011).
  * BELL (1966) J. S. BELL, Rev. Mod. Phys. 38, 447 (1966).
  * Brunner  _et al._ (2014) N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Rev. Mod. Phys. 86, 419 (2014).
  * Clauser  _et al._ (1969) J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Phys. Rev. Lett. 23, 880 (1969).
  * Klyachko  _et al._ (2008) A. A. Klyachko, M. A. Can, S. Binicioğlu, and A. S. Shumovsky, Phys. Rev. Lett. 101, 020403 (2008).
  * Kochen and Specker (1975) S. Kochen and E. P. Specker, “The problem of hidden variables in quantum mechanics,” in  _The Logico-Algebraic Approach to Quantum Mechanics: Volume I: Historical Evolution_, edited by C. A. Hooker (Springer Netherlands, Dordrecht, 1975) pp. 293–328.
  * Cabello  _et al._ (2014) A. Cabello, S. Severini, and A. Winter, Phys. Rev. Lett. 112, 040401 (2014).
  * Acín  _et al._ (2015) A. Acín, T. Fritz, A. Leverrier, and A. B. Sainz, Communications in Mathematical Physics 334, 533 (2015).
  * Araújo  _et al._ (2013) M. Araújo, M. T. Quintino, C. Budroni, M. T. Cunha, and A. Cabello, Phys. Rev. A 88, 022118 (2013).
  * Amaral  _et al._ (2014) B. Amaral, M. T. Cunha, and A. Cabello, Phys. Rev. A 89, 030101 (2014).
  * Kurzyński and Kaszlikowski (2012) P. Kurzyński and D. Kaszlikowski, Phys. Rev. A 86, 042125 (2012).
  * Lisoněk  _et al._ (2014) P. Lisoněk, P. Badzia¸g, J. R. Portillo, and A. Cabello, Phys. Rev. A 89, 042101 (2014).
  * Cabello  _et al._ (2015) A. Cabello, M. Kleinmann, and C. Budroni, Phys. Rev. Lett. 114, 250402 (2015).
  * Kleinmann  _et al._ (2012) M. Kleinmann, C. Budroni, J.-A. Larsson, O. Gühne, and A. Cabello, Phys. Rev. Lett. 109, 250402 (2012).
  * Cabello  _et al._ (2016) A. Cabello, M. Kleinmann, and J. R. Portillo, Journal of Physics A: Mathematical and Theoretical 49, 38LT01 (2016).
  * Yu and Oh (2012) S. Yu and C. H. Oh, Phys. Rev. Lett. 108, 030402 (2012).
  * Cabello (2013) A. Cabello, Phys. Rev. A 87, 010104 (2013).
  * Bharti  _et al._ (2019a) K. Bharti, M. Ray, and L.-C. Kwek, Entropy 21 (2019a), 10.3390/e21020134.
  * Ramanathan  _et al._ (2012) R. Ramanathan, A. Soeda, P. Kurzyński, and D. Kaszlikowski, Phys. Rev. Lett. 109, 050404 (2012).
  * Kurzyński  _et al._ (2014) P. Kurzyński, A. Cabello, and D. Kaszlikowski, Phys. Rev. Lett. 112, 100401 (2014).
  * Bharti  _et al._ (2019b) K. Bharti, M. Ray, A. Varvitsiotis, N. A. Warsi, A. Cabello, and L.-C. Kwek, Phys. Rev. Lett. 122, 250403 (2019b).
  * Singh  _et al._ (2017) J. Singh, K. Bharti, and Arvind, Phys. Rev. A 95, 062333 (2017).
  * Cabello  _et al._ (2011) A. Cabello, V. D’Ambrosio, E. Nagali, and F. Sciarrino, Phys. Rev. A 84, 030302 (2011).
  * Troupe and Farinholt (2015) J. E. Troupe and J. M. Farinholt, (2015), arXiv:1512.02256 [quant-ph] .
  * Saha  _et al._ (2019) D. Saha, P. Horodecki, and M. Pawłowski, New Journal of Physics 21, 093057 (2019).
  * Spekkens (2005) R. W. Spekkens, Phys. Rev. A 71, 052108 (2005).
  * Kunjwal and Ghosh (2014) R. Kunjwal and S. Ghosh, Phys. Rev. A 89, 042118 (2014).
  * Mazurek  _et al._ (2016) M. D. Mazurek, M. F. Pusey, R. Kunjwal, K. J. Resch, and R. W. Spekkens, Nature Communications 7, ncomms11780 EP (2016).
  * Kunjwal and Spekkens (2018) R. Kunjwal and R. W. Spekkens, Phys. Rev. A 97, 052110 (2018).
  * Kunjwal and Spekkens (2015) R. Kunjwal and R. W. Spekkens, Phys. Rev. Lett. 115, 110403 (2015).
  * Xu  _et al._ (2016) Z.-P. Xu, D. Saha, H.-Y. Su, M. Pawłowski, and J.-L. Chen, Phys. Rev. A 94, 062103 (2016).
  * Spekkens  _et al._ (2009) R. W. Spekkens, D. H. Buzacott, A. J. Keehn, B. Toner, and G. J. Pryde, Phys. Rev. Lett. 102, 010401 (2009).
  * Chailloux  _et al._ (2016) A. Chailloux, I. Kerenidis, S. Kundu, and J. Sikora, New Journal of Physics 18, 045003 (2016).
  * Ambainis  _et al._ (2019) A. Ambainis, M. Banik, A. Chaturvedi, D. Kravchenko, and A. Rai, Quantum Information Processing 18, 111 (2019).
  * Chaturvedi and Saha (2020) A. Chaturvedi and D. Saha, Quantum 4, 345 (2020).
  * Schmid  _et al._ (2020) D. Schmid, J. H. Selby, and R. W. Spekkens, (2020), arXiv:2009.03297 [quant-ph] .


# The Voigt and complex error function:
Humlíček’s rational approximation generalized

Franz Schreier
DLR — Deutsches Zentrum für Luft- und Raumfahrt, Institut für Methodik der
Fernerkundung, 82234 Oberpfaffenhofen, Germany E-mail<EMAIL_ADDRESS>
(Accepted 2018 June 20. Received 2018 June 18; in original form 2018 May 2)

###### Abstract

Accurate yet efficient computation of the Voigt and complex error function is
a challenge since decades in astrophysics and other areas of physics. Rational
approximations have attracted considerable attention and are used in many
codes, often in combination with other techniques. The 12-term code “cpf12” of
Humlíček (1979) achieves an accuracy of five to six significant digits
throughout the entire complex plane. Here we generalize this algorithm to a
larger (even) number of terms. The $n=16$ approximation has a relative
accuracy better than $10^{-5}$ for almost the entire complex plane except for
very small imaginary values of the argument even without the correction term
required for the cpf12 algorithm. With 20 terms the accuracy is better than
$10^{-6}$. In addition to the accuracy assessment we discuss methods for
optimization and propose a combination of the 16-term approximation with the
asymptotic approximation of Humlíček (1982) for high efficiency.

###### keywords:

Techniques: spectroscopic – Methods: numerical – Line: profiles

††pubyear: 2018††pagerange: The Voigt and complex error function:  Humlíček’s
rational approximation generalized–References

##  1 Introduction

The Voigt profile (Olver et al., 2010, print companion to the NIST Digital
Library of Mathematical Functions, http://dlmf.nist.gov/), i.e. the
convolution of a Lorentzian and Gaussian profile, is ubiquitous in many
branches of physics including astrophysics, e.g., Peraiah (2002); Heng (2017).
From a computational point it is more convenient to consider the Voigt
function $K(x,y)$ depending on two variables $x$ and $y$ (essentially the
distance from the center peak and the ratio of the Lorentz to Gauss width) and
the complex error function $w(x+\mathrm{i}y)$, whose real part gives the Voigt
function. Because there is no closed-form solution for the convolution
integral, numerous algorithms have been developed in the past for efficient
and/or accurate evaluation utilizing series or asymptotic expansions,
continued fractions, rational approximations, Gauss-Hermite quadrature etc.
(for an old review, still worth while to read, see Armstrong, 1967).

For an assessment of Voigt and complex error function algorithms the intended
application has to be considered. Highly accurate codes such as Poppe & Wijers
(1990a, b) (14 significant digits stated accuracy) and Zaghloul & Ali (2011)
(20 significant digits) or arbitrary precision codes (Boyer & Lynas-Gray,
2014; Molin, 2011) are indispensable for tests of an algorithm’s accuracy, but
not necessarily fast. High resolution line-by-line (lbl) radiative transfer
modeling (Bailey & Kedziora-Chudczer, 2012; Schreier et al., 2014) presents a
“million- to billion-line challenge” (Grimm & Heng, 2015; Heng, 2017) where
numerous Voigt functions have to be evaluated (Rothman et al., 2010; Tennyson
& Yurchenko, 2012) at thousands to millions of frequency grid points and
computational speed becomes a prime concern.

In view of the few digits used for line strengths and broadening parameters in
spectroscopic data bases such as HITRAN (Gordon et al., 2017), HITEMP (Rothman
et al., 2010), or GEISA (Jacquinet-Husson et al., 2016), four significant
digits for the Voigt function can be regarded as appropriate for many lbl
calculations (Schreier, 2011). On the other hand, the increasing quality of
spectroscopic measurements has indicated deficiencies of the Voigt profile.
More sophisticated line shapes accounting for collisional Dicke narrowing,
line mixing, or speed-dependence are required (Varghese & Hanson, 1984;
Tennyson et al., 2014) that can often be expressed in terms of the complex
error function, where the imaginary parts of the arguments can differ by more
than ten orders of magnitude (Tran et al., 2013; Schreier, 2017). For
“standard” Voigt lbl modeling in astrophysical and planetary atmospheric
spectroscopy $y$ can be as large as $10^{5}$ and as small as $10^{-8}$ (Lynas-
Gray, 1993; Wells, 1999; Tepper-García, 2006, 2007; Schreier, 2011).

Rational approximations are particularly attractive because they can be
implemented efficiently (i.e. only one division plus additions and
multiplications) and allow for high accuracy. Hui et al. (1978) presented a
single rational approximation with a fifth-degree numerator polynomial and a
sixth-degree denominator polynomial that should be “sufficiently accurate for
most applications over the whole complex plane”. Humlíček (1979) proposed a
12-term rational approximation “cpf12”; including a modification for small $y$
the maximum relative error is less than $5\cdot 10^{-6}$. Humlíček (1982)
suggested a division of the complex plane into four subdomains; using
appropriate rational approximations the “w4” code evaluates the complex error
function accurate to four significant digits. Weideman (1994) developed a
single rational approximation applicable in the entire complex plane whose
accuracy can be adjusted by selection of the number $N$ of terms. Kochanov
(2011) presented two sums with four and six fractions with complex
coefficients.

Problems of the Hui et al. (1978) code for small $y$ have been reported early
by Karp (1978) and Humlíček (1982) (see also Schreier, 2011, Fig. 4). The
Kochanov (2011) approximation shows significant errors for small $y<10^{-2}$.
Shippony & Read (1993, 2003) use the Hui et al. approach (for small $|x|$ and
intermediate $y$) along with other techniques. Humlíček’s codes appear to
belong to the most popular complex error function codes (according to Scopus
(March 2018), there are 218 and 313 citing articles for the 1979 and 1982
paper, respectively). In particular the w4 code has been further developed by
Kuntz (1997); Ruyten (2004); Imai et al. (2010) and the cpf12 code has been
implemented by Tran et al. (2013) for evaluation of the “Hartmann-Tran”
profile (Tennyson et al., 2014). Wells (1999) combined approximations of both
Humlíček codes with slightly modified region bounds. Zaghloul (2017) allows to
select the accuracy and uses Humlíček (1982) for low accuracy.

Except for Weideman (1994) all codes achieve moderate to high accuracy
throughout the whole complex plane only by a combination of several methods.
Unfortunately this can make the code difficult to optimize; this is especially
true if the conditional branches are more than a simple “either-or” as for
example in the Humlíček (1979) cpf12 algorithm. To avoid complicated (e.g.
nested) if structures we have combined the asymptotic rational approximation
of Humlíček (1982) for $|x|+y>15$ with the Weideman (1994) approximation
(Schreier, 2011).

Whereas the four rational approximations of Humlíček (1982) have been utilized
by several refinements, the somewhat older and slightly more accurate Humlíček
(1979) cpf12 code has rarely been used in other complex error function codes
(except for Wells (1999)). In particular we are not aware of any algorithm
using the Humlíček (1979) rational approximation with a higher (or smaller)
number of terms (i.e. $n\neq 12$).

The objective of this paper is an assessment of the Humlíček (1979) rational
approximation for an “arbitrary” number of terms. We briefly review some basic
facts in the following section, followed by a presentation of the Humlíček
algorithm, and a discussion of our optimization strategies. In Section 3 we
evaluate the accuracy and efficiency of the “generalized” Humlíček (1979)
code, and provide our conclusions in Section 4.

##  2 Theory and Methods

###  2.1 The Voigt and complex error function

The Voigt (or Hjerting) profile and the closely related Voigt function are
defined by

$\displaystyle
g_{\text{V}}(\nu-\hat{\nu},\gamma_{\text{L}},\gamma_{\text{G}})~{}$
$\displaystyle=~{}\int_{-\infty}^{\infty}\mathrm{d}\nu^{\prime}\;g_{\text{L}}(\nu-\nu^{\prime},\gamma_{\text{L}})\,\times\,g_{\text{G}}(\nu^{\prime}-\hat{\nu},\gamma_{\text{G}})$
(1) $\displaystyle=~{}{\sqrt{\ln 2/\pi}\over\gamma_{\text{G}}}~{}K(x,y)$
$\displaystyle K(x,y)~{}$
$\displaystyle=~{}{y\over\pi}~{}\int_{-\infty}^{\infty}{\mathrm{e}^{-t^{2}}\over(x-t)^{2}+y^{2}}~{}\mathrm{d}t~{},$
(2)

where $\gamma_{\text{L}}$ and $\gamma_{\text{G}}$ are the half widths at half
maximum (HWHM) of the Lorentzian $g_{\text{L}}$ and Gaussian $g_{\text{G}}$,
respectively, and $\hat{\nu}$ is the center wavenumber or frequency (i.e.
corresponding to the energy difference of an atomic or molecular transition).
The dimensionless arguments of the Voigt function are defined as ratios

$x~{}=~{}\sqrt{\ln 2}~{}{\nu-\hat{\nu}\over\gamma_{\text{G}}}\qquad\hbox{ and
}\qquad y~{}=~{}\sqrt{\ln 2}~{}{\gamma_{\text{L}}\over\gamma_{\text{G}}}~{}.$
(3)

Note that the profiles in Eq. (1) are normalized to one, $\int
g(\nu)\,\mathrm{d}\nu=1$, whereas the Voigt function is normalized to
$\sqrt{\pi}$.

The Voigt function is closely related to the complex error function (a.k.a.
complex probability function, Fadde(y)eva function, or plasma dispersion
function)

$w(z)~{}\equiv~{}K(x,y)\,+\,\mathrm{i}L(x,y)~{}=~{}{\mathrm{i}\over\pi}~{}\int_{-\infty}^{\infty}~{}{e^{-t^{2}}\over
z-t}~{}\mathrm{d}t$ (4)

with $z=x+\mathrm{i}y$. Asymptotically the complex error function behaves as
$\mathrm{i}/(\sqrt{\pi}z)$, corresponding to a Lorentz-like Voigt function for
large $x$. For $y=0$ the real part becomes a Gauss function, i.e.
$K(x,0)=\text{Re}\bigl{(}w(x)\bigr{)}=\exp(-x^{2})$. The half width at half
maximum of the Voigt function can be estimated as (Olivero & Longbothum, 1977;
Schreier, 2011)

$x_{1/2}~{}=~{}\dfrac{1}{2}\>\left(y\,+\,\sqrt{y^{2}+4\ln 2}\right)~{}.$ (5)

###  2.2 The rational approximation by Humlíček (1979)

Humlíček presented a rational approximation in the form

$w(z)~{}=~{}\sum_{\begin{subarray}{c}k=-n/2\\\ k\neq
0\end{subarray}}^{n/2}{\alpha_{k}+\mathrm{i}\beta_{k}\over
z-x_{k}+\mathrm{i}\delta}\\\
~{}=~{}\sum_{k=1}^{n/2}\left({\alpha_{k}+\mathrm{i}\beta_{k}\over
z-x_{k}+\mathrm{i}\delta}~{}-~{}{\alpha_{k}-\mathrm{i}\beta_{k}\over
z+x_{k}+\mathrm{i}\delta}\right)$ (6)

with real constants

$\displaystyle\alpha_{k}~{}$
$\displaystyle=~{}-\alpha_{-k}~{}=~{}-{1\over\pi}\omega_{k}\mathrm{e}^{\delta^{2}}\sin(2x_{k}\delta)$
(7) $\displaystyle\beta_{k}~{}$
$\displaystyle=~{}+\beta_{-k}~{}=~{}+{1\over\pi}\omega_{k}\mathrm{e}^{\delta^{2}}\cos(2x_{k}\delta)~{}.$
(8)

Here $x_{k}$ and $\omega_{k}$ are the roots and weights of the $n$-point
Gauss-Hermite formula (with $n$ even and the nodes and weights antisymmetric
and symmetric in $k$ (e.g. Abramowitz & Stegun, 1964, Table 25.10)), and
$\delta$ is a positive constant to be chosen appropriately. Note that the
constants $\alpha_{k}$, $\beta_{k}$, $x_{k}$, $\omega_{k}$, and $\delta$
depend on the number $n$ of nodes, but the superscript $(n)$ used in the
original paper has been omitted here; furthermore, Humlíček uses the notation
$\delta_{k}^{(n)}\equiv y_{0}^{(n)}$ synonymously. As emphasized by Humlíček,
“the choice of $n$ and $\delta$ involves a compromise between several
demands”, i.e. a larger number of terms increases the quality of the
approximation, but also the computational effort. Table 1 of Humlíček (1979)
gives some hints on the proper choice of $\delta$ for $8\leq n\leq 20$.

To circumvent accuracy problems near the real axis where (6) fails to
approximate an almost Gaussian function (see subsection 3.1), Humlíček
introduced a modification for the real part

$\displaystyle K(x,y)~{}=~{}\mathrm{e}^{-x^{2}}+$
$\displaystyle\sum_{\begin{subarray}{c}k=-n/2\\\ k\neq
0\end{subarray}}^{n/2}{y\over(x-x_{k})^{2}+\delta^{2}}\,\times$ (9)
$\displaystyle{\beta_{k}[(x-x_{k})^{2}-\delta(y+\delta)]-\alpha_{k}(x-x_{k})(y+2\delta)\over(x-x_{k})^{2}+(y+\delta)^{2}}$

for $y<0.85$ and $|x|>18.1y+1.65$. The “cpf12” code with $n=12$ and
$\delta=1.5$ achieves a maximum relative error less than $2\cdot 10^{-6}$ and
$5\cdot 10^{-6}$ for the real and imaginary part, respectively.

###  2.3 Optimization and a combination with the Humlíček (1982) asymptotic
approximation

Equations (6) and (9) are suboptimal w.r.t. computational efficiency because
of the numerous divisions: According to Ueberhuber (1997) divisions are
significantly slower than additions or multiplications, and Goedecker & Hoisie
(2001) noted that “the calculation of special functions, such as divisions,
square roots, exponentials and logarithms requires anywhere from a few dozen
cycles up to hundreds of cycles … these calculations have to be decomposed
into a sequence of elementary instructions such as multiplies and adds.”
However, using elementary calculus (6) can be rewritten to the “standard form”
of a rational approximation

$w(z)~{}=~{}{P(z)\over
Q(z)}~{}=~{}{\sum\limits_{k=0}^{n-1}a_{k}z^{k}\over\sum\limits_{l=0}^{n}b_{l}z^{l}}$
(10)

Note that $a_{n-1}=\mathrm{i}/\sqrt{\pi}$ and $b_{n}=1$ in accordance with the
asymptotic expansion of the complex error function. The numerator coefficients
$a_{k}$ are real for even $k$ and purely imaginary for odd $k$; furthermore
$b_{l}=0$ for all odd $l$. Finally it is worth to mention that the
coefficients depend on the parameter $\delta$.

For the 16-term approximation transformed into the form (10) the number of
adds and multiplies is similar to the $N=32$ rational approximation of
Weideman (1994). However, the benchmark tests in Schreier (2011) have
indicated that this approximation is significantly slower than the asymptotic
rational approximation

$R_{1,2}(z)={\mathrm{i}z/\sqrt{\pi}\over z^{2}-{1\over 2}}$ (11)

used by Humlíček (1982) for $s=|x|+y>15$. (The subscript on the lefthand side
indicates the degree of the numerator and denominator polynomial.) Note that
for lbl calculations function values for only a few grid points have to be
evaluated in the line center, and most values have to be evaluated in the line
wings with large $|x|$. Accordingly we had suggested to use the Weideman
approximation ($N=24$) near the line center only and approximation (11)
otherwise. Similarly we propose a combination of the two Humlíček (1979, 1982)
approximations

$w(z)~{}=~{}\begin{cases}{\mathrm{i}z/\sqrt{\pi}\over z^{2}-{1\over
2}}&|x|+y>15\\\
{\sum\limits_{k=0}^{n-1}a_{k}z^{k}\over\sum\limits_{l=0}^{n}b_{l}z^{l}}\qquad&\text{otherwise.}\end{cases}$
(12)

##  3 Results

In the first subsection we provide an assessment of the accuracy of the
Humlíček rational approximation (6) for $n\geq 12$. As an accuracy reference
SciPy’s wofz function is used, a combination of the Poppe & Wijers (1990a, b)
and Zaghloul & Ali (2011) codes with a stated accuracy of at least 13
significant digits (Scientific Python scipy.special.wofz implementation based
on S.G. Johnson’s “Faddeeva package”, http://ab-initio.mit.edu/Faddeeva).
Preliminary timing tests are presented in Subsection 3.2, and the performance
of lbl cross section modeling is discussed in the third subsection. The
combination of (6) and (9) with $n=12$ is denoted as cpf12, whereas the
optimized form (10) with $n=12$ or $n=16$ is denoted as zpf12 and zpf16,
respectively.

Figure 1: Comparison of the Humlíček cpf12 approximation for $\delta=1.5$ with
the wofz reference. Left: approximation (6); Right: (6) and (9) combined.
White areas indicate a relative error better than $10^{-12}$. The numbers in
the title indicate the maximum relative error of $K$ and $L$ (in parentheses).

###  3.1 Accuracy

Fig. 1 depicts the relative error
$\epsilon(x,y)=|K_{\text{cpf}}-K_{\text{wofz}}|/K_{\text{wofz}}$ of the real
part of the $n=12$ approximation. Fig. 1a clearly confirms Humlíček’s
statement w.r.t. accuracy problems of the approximation (6) for small $y$ (the
imaginary part $L(x,y)$ is computed correctly with a relative accuracy better
than $4\cdot 10^{-6}$). Including the correction term (9) drastically reduces
the errors, and the maximum relative error of about $10^{-6}$ is in accordance
with Humlíček’s claims. However, the drawback of the correction term is the
difficulty to implement this efficiently.

Figure 2: Comparison of Humlíček zpf16 approximation for $\delta=1.3118$
(left) and $\delta=1.35$ (right) with the wofz reference.

The $n=16$ rational approximation is compared to the wofz reference code in
Fig. 2. Variation of $\delta$ indicates that for $\delta=1.3118$ the
approximation is optimal w.r.t. overall accuracy, i.e. the relative error is
$\leq 7.86\cdot 10^{-5}$; for $y>10^{-6}$ the error is less than about
$10^{-5}$ for all $x,y$ and the largest error is observed at $x\approx 4.8$
and $y=10^{-8}$. Interestingly, the imaginary part $L(x,y)$ is slightly less
accurate compared to the $n=12$ approximation. For $\delta\approx 1.33$ (and
very small $y$) the location of the maximum error moves into the line wings.
For $\delta=1.35$ the maximum of the relative error is slightly larger, about
$2\cdot 10^{-4}$, but this maximum is achieved for large $x$ and very small
$y$ where the function values are already tiny ($K(5,10^{-6})\approx
10^{-8}$). For the imaginary part the maximum relative deviation is $2.6\cdot
10^{-6}$.

In Fig. 3a we compare the maximum relative error $\max_{x}\epsilon(x,y)$ as a
function of $y$. For $y>10^{-4}$ the approximation (6) or (10) is accurate to
four significant digits for all $\delta$ considered.

Figure 3: The maximum relative error of the Humlíček $n=16$ and $n=20$
approximations for various $\delta$ compared to the wofz reference.

Increasing the number of terms in the rational approximation (6) further
improves the accuracy. Fig. 3b shows that six significant digits can be
achieved with $n=20$ and $\delta\approx 1.55$. Contour plots of the relative
error of the $n=18$ and $n=20$ approximations with optimal $\delta$ are shown
in Fig. 4.

Figure 4: Comparison of Humlíček $n=18$ ($\delta=1.45$, left) and $n=20$
($\delta=1.55$, right) approximation with the wofz reference.

For a further appraisal of the generalized Humlíček (1979) approximation (6)
we evaluate the Voigt function (i.e. $\text{Re}(w))$ for selected values of
$z=x+\mathrm{i}y$ and compare these with the results of Lether & Wenston
(1991) who use a corrected midpoint quadrature rule. The values of $K(x,y)$
compiled in their Table 2 “are believed to be correct to twenty-five
significant digits” and have also been cross-checked by Zaghloul (2007, his
table 2) and Boyer & Lynas-Gray (2014, table 1). Table 1 shows that for $n=16$
and $n=20$ the expansion (6) is correct to at least five significant digits
for all $x,y$ pairs except for $x=y=10^{-3}$. Note that for $y=10^{-20}$
(first data row) the Voigt function is essentially a Gaussian, and the correct
value $\exp(-1.0)=0.36787\;94411\;71442$ is reproduced with six, eight, and
nine digits for $n=16$, $n=20$, and $n=24$, respectively. In the last three
rows with moderate $x$ and $y$ the cpf values are correct to ten digits.

Table 1: Comparison of Voigt function values for the $n=16$, $n=20$, and $n=24$ approximation (6) with the precise values of Lether & Wenston (1991, their table 2 (given with 25 digits)). $x$ | $y$ | Lether & Wenston | $n=16~{}(\delta=1.35)$ | $n=20~{}(\delta=1.55)$ | $n=24~{}(\delta=1.4)$
---|---|---|---|---|---
1.0 | $10^{-20}$ | 0.36787 94411 71442 | 0.36787 93403 59605 | 0.36787 94396 18318 | 0.36787 94414 03711
10. | $10^{-4}$ | 0.57287 17561 64533 $\cdot 10^{-6}$ | 0.57287 17507 11831 $\cdot 10^{-6}$ | 0.57287 17561 87614 $\cdot 10^{-6}$ | 0.57287 17561 66040 $\cdot 10^{-6}$
$10^{-3}$ | $10^{-3}$ | 0.99887 16233 35411 | 0.99885 13083 62977 | 0.99886 16184 44897 | 0.99887 16170 22162
0.0 | 0.25 | 0.77034 65477 30996 | 0.77034 65303 12796 | 0.77034 65475 77724 | 0.77034 65475 71460
1.0 | 0.50 | 0.35490 03328 67577 | 0.35490 03328 53826 | 0.35490 03328 64360 | 0.35490 03328 66028
5.0 | 5.0 | 0.56965 43988 71769 $\cdot 10^{-6}$ | 0.56965 43988 81711 $\cdot 10^{-1}$ | 0.56965 43988 81771 $\cdot 10^{-1}$ | 0.56965 43988 81769 $\cdot 10^{-1}$
1.0 | 10.0 | 0.55598 31964 10553 $\cdot 10^{-6}$ | 0.55598 31964 10505 $\cdot 10^{-1}$ | 0.55598 31964 10555 $\cdot 10^{-1}$ | 0.55598 31964 10553 $\cdot 10^{-1}$

###  3.2 Preliminary Speed Tests with IPython

To get a first impression about the performance of the various approximations,
we have used IPython’s builtin “magic” function %timeit. A typical
“experiment” looks like

    In [1]: from cpfX import zpf16p
    In [2]: x=numpy.linspace(0.,100.,10001);  y=0.001;
            z=x+1j*y
    In [3]: %timeit zpf16p(z)
    1000 loops, best of 3: 1.04 ms per loop

Table 2 summarizes the results of these experiments. Note that the grid point
spacing $\delta x=0.1$ (first three columns) is more than adequate to sample
the line center region as the half width (5) is approximately $x_{1/2}\approx
1$ for small $y$.

For the Python translation of the original Fortran cpf12 code two versions
have been implemented: A scalar function with an if statement that is
transformed in an array function using NumPy’s vectorize function class and an
array function exploiting a mask and NumPy’s where. The first two rows of the
table clearly indicate that the vectorized version is significantly slower.

The polynomials in (10) can be implemented either with the Horner scheme coded
by hand (zpf12h, zpf16h) or exploiting NumPy’s poly1d class (zpf12p, zpf16p),
and the %timeit tests demonstrate that the hand-coded Horner scheme is
slightly faster for all cases. It should also be noted that implementing the
denominator as a polynomial of $z^{2}$ (the odd coefficients are zero) is more
efficient than the polynomial of $z$.

For comparison we also report times for several Weideman (1994)
approximations. His $N$-term approximation requires $N+4$ multiplications plus
one division, so the computational effort for weideman32 is similar to zpf16.
However, as mentioned above, the $16$-degree denominator polynomial of zpf16
can be implemented more efficiently because of the vanishing odd coefficients,
i.e. with 24 add-multiplies for zpf16 similar to weideman24.

The wofz implementation documented in the very last row appears to be faster
than all rational approximations. Note, however, that this special function
imported from SciPy is actually a C code linked to Python in contrast to our
pure Numeric Python codes.

Note that the execution time is approximately proportional to the number of
function evaluations only for the cpf12v case (first data row) that is by far
the slowest implementation. Furthermore, times in the fourth column
($n_{x}=10\,000$) are roughly ten times larger compared to the third column
($n_{x}=1000$). The deviations from strict proportionality might indicate that
%timeit estimates also include overhead for function call etc.

Table 2: Execution time measured by the % timeit function in the IPython interpreter. For all columns $y=0.001$ has been used. $n_{x}$ is the number of $x$ grid points and function evaluations. The “zpf” indicates implementations of the optimized rational approximation (10). The tests have been performed on a laptop with an Intel x86_64 CPU running at 2.7 GHz and cache size 3072 KB. | $0\leq x\leq 25$ | $0\leq x\leq 50$ | $0\leq x\leq 100$ | $0\leq x\leq 100$
---|---|---|---|---
$n_{x}$ | 251 | 501 | 1001 | 10001
cpf12v | 12.1 ms | 24.8 ms | 49.4 ms | 492 ms
cpf12w | 877 $\mu$s | 1.01 ms | 1.17 ms | 4.39 ms
zpf12h | 64.2 $\mu$s | 84 $\mu$s | 119 $\mu$s | 829 $\mu$s
zpf12p | 116 $\mu$s | 145 $\mu$s | 185 $\mu$s | 1.01 ms
zpf16h | 90.1 $\mu$s | 110 $\mu$s | 147 $\mu$s | 904 $\mu$s
zpf16p | 117 $\mu$s | 146 $\mu$s | 185 $\mu$s | 1.04 ms
weideman16 | 78.5 $\mu$s | 100 $\mu$s | 134 $\mu$s | 776 $\mu$s
weideman24 | 105 $\mu$s | 138 $\mu$s | 172 $\mu$s | 1.04 ms
weideman32 | 132 $\mu$s | 162 $\mu$s | 216 $\mu$s | 1.31 ms
wofz | 52.6 $\mu$s | 66.8 $\mu$s | 91.7 $\mu$s | 886 $\mu$s

###  3.3 Preformance of lbl cross section modeling

For realistic benchmarks we test the generalized Humlíček approximation (10)
in an atmospheric radiative transfer lbl modeling context. Radiative transfer
in Earth’s atmosphere is clearly important in the geosciences, esp. for
climate modeling and remote sensing, but is also relevant for astronomy, i.e.
for corrections of telluric absorption and emission in ground-based
observations (e.g. Seifahrt et al., 2010; Bertaux et al., 2014; Smette et al.,
2015) or radiative transfer modeling of Earth-like exoplanets (e.g. Des Marais
et al., 2002; Robinson & Reinhard, 2018). To make comparison with our previous
studies (Schreier & Kohlert, 2008; Schreier, 2011) easier, we compute high
resolution molecular absorption cross sections for the SubMillimeter
Radiometer (SMR) aboard the Swedish small satellite ODIN (Murtagh et al.,
2002; Nordh et al., 2003), an aeronomy and astronomy mission launched in 2001
(the astronomy mission was successfully concluded in spring 2007, see also
Hjalmarson et al. (2003)).

The cross sections

$k(\nu,p,T)~{}=~{}\sum\limits_{l}S_{l}(T)\>g_{\text{V}}\Bigl{(}\nu-\hat{\nu}_{l},\gamma_{l}^{\text{(L)}}(p,T),\gamma_{l}^{\text{(G)}}(T)\Bigr{)}~{}.$
(13)

(with line position $\hat{\nu}_{l}$ and line strength $S_{l}$) are computed in
the $16\,\text{--}\,17\rm\,cm^{-1}$ interval (480 – 510 GHz) for thirteen
pressure-temperature pairs (corresponding to the $0\,\text{--}\,120\rm\,km$
altitude range with $\delta z=10\rm\,km$ steps). Because of the importance of
line wing contributions all lines in the extended
$6\,\text{--}\,27\rm\,cm^{-1}$ interval are considered: For nitric acid
($\mathrm{HNO_{3}}$) HITRAN 00 (Rothman et al., 2003) lists 2376 lines (HITRAN
16 (Gordon et al., 2017) has about 250 000 lines), and there are 21 565 lines
of ozone ($\mathrm{O_{3}}$) in all versions of HITRAN since 1996. The mean
Lorentz to Doppler (Gauss) width ratio $y$ decreases from $7\cdot 10^{3}$ (at
bottom-of-atmosphere, BoA) to $1.4\cdot 10^{-4}$ (at top-of-atmosphere, ToA)
for $\mathrm{HNO_{3}}$; the Doppler width is slightly larger for
$\mathrm{O_{3}}$ because of the smaller mass, hence $y$ is somewhat smaller.

Figure 5: Execution time (ns) per function evaluation. The tests have been
performed on an idle node of a Linux cluster with an Intel XEON CPU at 2.36
GHz.

####  3.3.1 $\mathrm{HNO_{3}}$ cross sections — Python

The $\mathrm{HNO_{3}}$ cross sections are modeled using the lbl2xs.py function
script of Py4CAtS — Python for Computational Atmospheric Spectroscopy
(available at https://atmos.eoc.dlr.de/tools/Py4CAtS/) using a “brute force”
approach, i.e. without further approximations in the line wings (e.g. coarse
grid). Note that the spectral resolution is depending on line width
(essentially pressure, hence altitude), i.e. the grid point spacing $\delta x$
is set to a fraction of the mean half width, typically $\delta
x=\bar{x}_{1/2}/4$ or $\delta\nu=\bar{\gamma}_{\text{V}}/4$. The number of
wavenumber grid points increases from less than one hundred at BoA to almost
half a million at an altitude of $90\rm\,km$.

The execution times for a single function value (i.e. the total elapsed CPU
time measured with the clock function of Python’s time module divided by the
number of lines and the number of $x$ grid points) shown in Fig. 5 essentially
confirm the IPython %timeit tests reported in the previous subsection. Note
that pre- and postprocessing steps (e.g. reading the line data) have not been
considered for the timing.

Both versions of the cpf12 code are significantly slower than all other
functions: the cpf12w code needs about 600 ns for $p\leq 10\rm\,mb$ and is not
shown. Except for large pressures the various Weideman functions and the two
implementations of the approximation (10) with $n=16$ require about $100$ to
$175\rm\,ns$, and the use of an asymptotic approximation for large $x$ is
tempting. The performance of the combination (12) as well as the Humlíček
(1982)–Weideman (1994) combination (with $N=24$) are almost identical, and
somewhat worse to the wofz C function linked to Python. The combination of the
$R_{3,4}$ rational approximation of Humlíček (1982) with the Weideman $N=32$
approximation (with a cut at $s=8$, cf. Schreier (2017)) has a relative error
better than $2\cdot 10^{-6}$ and is slightly slower.

The good performance of the wofz C function reported in the previous
subsection suggests that implementation of compute intensive code segments in
a compiled language might be advantageous. Accordingly we have used the f2py
Fortran-Python interface generator (Peterson, 2009) to make calls to Fortran
subroutines possible. Both combinations of the Humlíček (1982) $R_{1,2}$
asymptotic approximation (11) with the Humlíček (1979) zpf16 approximation
(10) or with the 24-term Weideman approximation (denoted hum1zpf16_f2py and
hum1wei24_f2py) are more than a factor two faster than the corresponding pure
Python implementation and clearly superior with respect to computational
efficiency.

Compared to the tests reported in Schreier (2011, Fig. 11) the
Humlíček–Weideman as well as the pure Weideman functions are about a factor 10
faster now. Note that the minimum times seen here are comparable to the
Fortran execution times shown in Fig. 10 of Schreier (2011). This considerable
speed-up is presumably a combination of better hardware and better software
(improved performance of NumPy).

####  3.3.2 $\mathrm{O_{3}}$ cross sections — Fortran

The Humlíček (1979) $n=16$ approximation (10) along with several other complex
error function algorithms has also be implemented in Fortran 90. The execution
time for the lbl modeling of $\mathrm{O_{3}}$ cross sections is measured by
means of the cpu_time intrinsic function.

In contrast to the Python tests the runtime appears to be fairly constant for
pressures smaller than about $10\rm\,mb$ where the number of wavenumber grid
points is already large ($n_{x}>10^{4}$). Hence we only report the average run
time in Table 3. cpf12 roughly corresponds to the original code documented in
the Humlíček (1979) paper and Fortran can execute this code significantly
faster compared to NumPy. The optimized implementation (10) of the 16-term
approximation zpf16 is executed in Fortran with roughly double speed. The
combinations of the zpf16 approximation (10) or Weideman’s approximation with
the Humlíček (1982) $R_{12}$ or $R_{34}$ approximation for large $|x|$ have
approximately identical run times. Comparison with Figure 5 reveals that the
speed of the Fortran and NumPy–f2py codes is roughly similar, i.e. $15\rm\,ns$
per function value. The two versions of the GNU gfortran compiler produce code
with approximately the same performance, whereas the NAG compiled code is
somewhat slower (similar to Schreier (2011, Fig. 9)).

There appears to be no pure Fortran array implementation of wofz. The libcerf
numeric library (http://apps.jcns.fz-juelich.de/libcerf) implementation of the
Faddeeva package also provides Fortran bindings, but only for scalar
arguments. However, the execution time of the hum1zpf16 and hum1wei24 pure
Fortran subroutines and the functions wrapped to Python via f2py are almost
identical, so we can expect that a Fortran array version of wofz would perform
similar to the SciPy-C version that was about a factor 1.5 slower compared to
hum1zpf16 and hum1wei24 (compare Figure 5 and Table 3).

Table 3: Execution time [ns] per function evaluation for the Fortran 90 implementation. GNU Fortran or NAG compiler with optimization flag -O3 (the numbers indicate the version). Linux cluster node as in Fig. 5. | gfortran 4.8 | gfortran 6.2 | nagfor 5.3
---|---|---|---
cpf12 | 97 | 102.6 | 125
zpf16 | 64.6 | 61.2 | 66.24
hum1zpf16 | 14.52 | 14.61 | 16.98
hum1wei24 | 14.55 | 14.49 | 16.39
hum2wei32 | 14.75 | 14.86 | 18.15
weideman16 | 70.34 |  |
weideman32 | 131.1 |  |

##  4 Summary and Conclusions

Generalizations of the Humlíček (1979) cpf12 rational approximation (with 12
terms) for the Voigt and complex error function have been presented and
optimized implementations have been developed. The rational approximation with
a degree 16 denominator polynomial is accurate to at least five significant
digits except for very small $y$, and somewhat larger deviations show up in
the wings where the function value is already extremely small. For $n\geq 16$
the correction term needed to compensate for the problems of the $n=12$
approximation (6) for large $|x|$ and small $y$ is not required anymore, i.e.
a single rational approximation can be used to evaluate the complex error
function over almost the entire complex plane.

The sum of $n$ fractions of the original approximation of Humlíček (1979) are
computationally inefficient, but can be readily transformed into a single
fraction with an $n-1$ degree numerator polynomial and an $n$ degree
denominator polynomial. The performance of the original and generalized
Humlíček algorithms implemented as Numeric Python functions has been tested
and compared to other rational approximations. For a more realistic scenario
molecular absorption cross sections have been computed line-by-line using
Python and Fortran implementations. The combination of the $n=16$
approximation with the asymptotic approximation also developed by Humlíček
(1982) further improved the performance. The combination of the Humlíček
(1979, 1982) methods as well as the Humlíček–Weideman combination need about
$15\rm\,ns$ for a single function value both with Fortran and Python, in
strong contrast to our findings seven years ago (Schreier, 2011) where Fortran
was about a factor ten faster. This impressive speed of the Python
implementation on a par with the Fortran code explains the increasing
popularity of the “Python-based ecosystem” (www.scipy.org) and packages such
as Astropy (Astropy Collaboration, 2013).

In conclusion, the combinations of the Humlíček (1979) or Weideman (1994)
approximation with the Humlíček (1982) asymptotic approximation can be
recommended for both accurate (five or more digits) and efficient Voigt and
complex error function evaluations. For high performance with Python, a
Fortran (or C) implementation linked to Python is advantageous. When speed is
not an issue, the generalization of the original cpf12 code to 16 or more
terms can be used for all $x,y$ without any conditional branches.

## Supplementary material

Two Python source files are provided as supplementary material with the online
version of this paper. Additionally we provide the Fortran code used for the
tests of subsection 3.3.2 including the Humlíček (1979, 1982) combination
(12).

cpfX.py Various implementations of the Humlíček (1979) algorithm:

  * •

Python/NumPy implementations of the original Fortran 77 source code (with two
versions to combine the two regions);

  * •

A generalization of the region I rational approximation to arbitrary (even)
number of terms;

  * •

An optimized implementation of the region I approximation for 16 terms using a
single fraction;

  * •

The combination (12) for NumPy.

fig124.py Execution of this script inside the IPython interpreter should
produce the Figures 1, 2, and 4 of the manuscript.

Test_voigt_speed.f90 The Fortran 90 program including subroutines of various
complex error function algorithms.

## Acknowledgements

Financial support by the DFG project SCHR 1125/3-1 is greatly appreciated. I
would also like to thank Thomas Trautmann for critical reading of the
manuscript. The transformation of the sum of quotients (6) to the optimized
form (10) has been computed with SymPy (http://www.sympy.org/).

## References

  * Abramowitz & Stegun (1964) Abramowitz M., Stegun I., 1964, Handbook of Mathematical Functions. National Bureau of Standards, AMS55, New York
  * Armstrong (1967) Armstrong B., 1967, JQSRT, 7, 61
  * Astropy Collaboration (2013) Astropy Collaboration T., Robitaille T., Tollerud E., et al., 2013, A&A, 558, A33
  * Bailey & Kedziora-Chudczer (2012) Bailey J., Kedziora-Chudczer L., 2012, MNRAS, 419, 1913
  * Bertaux et al. (2014) Bertaux L., Lallement R., Ferron S., Boonne C., Bodichon R., 2014, A&A, 564, A46
  * Boyer & Lynas-Gray (2014) Boyer W., Lynas-Gray A., 2014, MNRAS, 444, 2555
  * Des Marais et al. (2002) Des Marais D., et al., 2002, Astrobiology, 2, 153
  * Goedecker & Hoisie (2001) Goedecker S., Hoisie A., 2001, Performance Optimization of Numerically Intensive Codes. SIAM, Philadelphia, PA
  * Gordon et al. (2017) Gordon I., et al., 2017, JQSRT, 203, 3
  * Grimm & Heng (2015) Grimm S. L., Heng K., 2015, ApJ, 808, 182
  * Heng (2017) Heng K., 2017, Exoplanetary Atmospheres — Theoretical Concepts and Foundations. Princeton University Press
  * Hjalmarson et al. (2003) Hjalmarson A., et al., 2003, A&A, 402, L39
  * Hui et al. (1978) Hui A., Armstrong B., Wray A., 1978, JQSRT, 19, 509
  * Humlíček (1979) Humlíček J., 1979, JQSRT, 21, 309
  * Humlíček (1982) Humlíček J., 1982, JQSRT, 27, 437
  * Imai et al. (2010) Imai K., Suzuki M., Takahashi C., 2010, Adv. Space Res., 45, 669
  * Jacquinet-Husson et al. (2016) Jacquinet-Husson N., et al., 2016, J. Mol. Spectrosc., 327, 31
  * Karp (1978) Karp A., 1978, JQSRT, 20, 379
  * Kochanov (2011) Kochanov V., 2011, Atmospheric and Oceanic Optics, 24, 432
  * Kuntz (1997) Kuntz M., 1997, JQSRT, 57, 819
  * Lether & Wenston (1991) Lether F., Wenston P., 1991, J. Comp. Appl. Math., 34, 75
  * Lynas-Gray (1993) Lynas-Gray A., 1993, Comp. Phys. Comm., 75, 135
  * Molin (2011) Molin P., 2011, Multi-precision computation of the complex error function, Preprint at https://hal.archives-ouvertes.fr/hal-00580855
  * Murtagh et al. (2002) Murtagh D., et al., 2002, Can. J. Phys., 80, 309
  * Nordh et al. (2003) Nordh H. L., et al., 2003, A&A, 402, L21
  * Olivero & Longbothum (1977) Olivero J., Longbothum R., 1977, JQSRT, 17, 233
  * Olver et al. (2010) Olver F., Lozier D., Boisvert R., Clark C., eds, 2010, NIST Handbook of Mathematical Functions. Cambridge University Press, New York, NY
  * Peraiah (2002) Peraiah A., 2002, An Introduction to Radiative Transfer: Methods and Applications in Astrophysics. Cambridge University Press
  * Peterson (2009) Peterson P., 2009, Int. J. Comp. Sci. Eng., 4, 296
  * Poppe & Wijers (1990a) Poppe G., Wijers C., 1990a, ACM TOMS, 16, 38
  * Poppe & Wijers (1990b) Poppe G., Wijers C., 1990b, ACM TOMS, 16, 47
  * Robinson & Reinhard (2018) Robinson T., Reinhard C., 2018, preprint,  (arXiv:1804.04138)
  * Rothman et al. (2003) Rothman L., et al., 2003, JQSRT, 82, 5
  * Rothman et al. (2010) Rothman L., et al., 2010, JQSRT, 111, 2139
  * Ruyten (2004) Ruyten W., 2004, JQSRT, 86, 231
  * Schreier (2011) Schreier F., 2011, JQSRT, 112, 1010
  * Schreier (2017) Schreier F., 2017, JQSRT, 187, 44
  * Schreier & Kohlert (2008) Schreier F., Kohlert D., 2008, Comp. Phys. Comm., 179, 457
  * Schreier et al. (2014) Schreier F., Gimeno García S., Hedelt P., Hess M., Mendrok J., Vasquez M., Xu J., 2014, JQSRT, 137, 29
  * Seifahrt et al. (2010) Seifahrt A., Käufl H., Zängl G., Bean J., Richter M., Siebenmorgen R., 2010, A&A, 524, A11
  * Shippony & Read (1993) Shippony Z., Read W., 1993, JQSRT, 50, 635
  * Shippony & Read (2003) Shippony Z., Read W., 2003, JQSRT, 78, 255
  * Smette et al. (2015) Smette A., et al., 2015, A&A, 576, A77
  * Tennyson & Yurchenko (2012) Tennyson J., Yurchenko S. N., 2012, MNRAS, 425, 21
  * Tennyson et al. (2014) Tennyson J., et al., 2014, Pure Appl. Chem., 86, 1931
  * Tepper-García (2006) Tepper-García T., 2006, MNRAS, 369, 2025
  * Tepper-García (2007) Tepper-García T., 2007, MNRAS, 382, 1375
  * Tran et al. (2013) Tran H., Ngo N., Hartmann J.-M., 2013, JQSRT, 129, 199
  * Ueberhuber (1997) Ueberhuber C., 1997, Numerical Computation. Springer
  * Varghese & Hanson (1984) Varghese P., Hanson R., 1984, Appl. Opt., 23, 2376
  * Weideman (1994) Weideman J., 1994, SIAM J. Num. Anal., 31, 1497
  * Wells (1999) Wells R., 1999, JQSRT, 62, 29
  * Zaghloul (2007) Zaghloul M. R., 2007, MNRAS, 375, 1043
  * Zaghloul (2017) Zaghloul M., 2017, ACM TOMS, 44, 22:1
  * Zaghloul & Ali (2011) Zaghloul M., Ali A., 2011, ACM TOMS, 38, 15:1