# Urban Crowdsensing using Social Media: An Empirical Study on Transformer and
Recurrent Neural Networks

Jerome Heng1, Junhua Liu1 and Kwan Hui Lim1  1Information Systems Technology
and Design Pillar, Singapore University of Technology and Design
Email: {jerome_heng<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

An important aspect of urban planning is understanding crowd levels at various
locations, which typically require the use of physical sensors. Such sensors
are potentially costly and time consuming to implement on a large scale. To
address this issue, we utilize publicly available social media datasets and
use them as the basis for two urban sensing problems, namely event detection
and crowd level prediction. One main contribution of this work is our
collected dataset from Twitter and Flickr, alongside ground truth events. We
demonstrate the usefulness of this dataset with two preliminary supervised
learning approaches: firstly, a series of neural network models to determine
if a social media post is related to an event and secondly a regression model
using social media post counts to predict actual crowd levels. We discuss
preliminary results from these tasks and highlight some challenges.

††publicationid: pubid: 978-1-7281-6251-5/20/$31.00 © 2020 IEEE

##  I Introduction

With increasing urbanization in the world, it is important for urban planners
to understand crowd levels for better urban planning [1]. Many users have made
their social media posts publicly accessible and this can potentially be a
source of useful insights for urban analytics, such as identifying events,
detecting crowd levels and other environmental phenomenon [2, 3]. While
another possible solution is to employ the use of physical sensors, such
sensors are potentially costly and time consuming to implement on a large
scale.

The challenge with using social media is that it is often unstructured and
sparse, hence great care must be taken to filter the dataset such that the
information is relevant, and aggregation functions can be used in appropriate
places to reduce sparsity [4]. Existing solutions are based on social media
filtering and clustering techniques, e.g., sending the data through a pipeline
to detect clusters. Researchers at Twitter have constructed a real-time
pipeline to detect clusters of trending tweets and topics [5]. Above a certain
threshold, these clusters are identified as events/outliers.

We intend to extend and improve on such approaches for urban crowdsensing
using social media. In this paper, we make the following contributions:

  * •

We curate and introduce two datasets that can be used for various urban
sensing applications, particularly for detecting crowd levels and identifying
urban events with the corresponding ground truth labels (Section II). In
developing these datasets, we also identify the limitations of social media
data and the steps we took to overcome these limitations. This dataset is
publicly available at https://github.com/kraftedcheese/crowdsensing_datasets.

  * •

We proposed two problems that can be investigated on this dataset. The first
problem is an event detection problem where we determine if a specific social
media post is related to an urban event that is characterized by higher human
traffic at a given location and time. The second problem is a crowd level
prediction problem that is modelled as a regression problem where we predict
the relative counts of people in various areas of a city from social media
posts in that area.

  * •

For our first problem on event detection, we experiment with a set of
transformer and recurrent neural networks models and report some preliminary
results. For the second problem on crowd level prediction, we utilize a
regression model using social media post counts to predict actual crowd
levels.

##  II Dataset

Geo-tagged social media has been frequently used for solving various location-
based problems, such as predicting the location of a social media posts [6],
predicting next visit locations [7] and recommending tour itineraries [8],
among others. Similarly, we curate our dataset from mainly geo-tagged Twitter
messages (i.e., tweets) and Flickr photos, which we retrieve using the Twitter
and Flickr APIs, respectively. Tables I and II gives an outline of our curated
dataset of Twitter geo-tagged messages and Flickr geo-tagged photos.

In the following sections, we will describe them in more details as we
demonstrate their application for an event detection task and crowd prediction
task.

Date | 2014-01-12 to 2018-02-12
---|---
Number of Tweets | 1,393,561
Number of Unique Users | 77,642
Number of Events | 22,241
Number of Non-events | 1,371,320
TABLE I: Dataset for Event Detection |  | Number of | Number of | Number of
---|---|---|---|---
Name | Date | Data Points | crowd sensors | unique users
Sensor Counts | 2009-05-01 to 2020-04-30 | 3,132,346 | 66 | -
Tweets (Sensor-related) | 2010-09-12 to 2018-06-04 | 266,931 | - | 20,176
Tweets (max 100 duplicate coords.) | 2010-09-21 to 2018-02-02 | 24,638 | 42 | 5,659
Tweets (max 10 duplicate coords.) | 2010-09-21 to 2018-02-02 | 12,801 | 42 | 3,785
Flickr (near to Town Hall (West)) | 2010-07-08 to 2020-03-29 | 6,076 | - | 852
TABLE II: Dataset for Crowd Level Prediction | Weighted Average | "Event" Class | "Non-Event" Class
---|---|---|---
Model | Precision | Recall | F1 | Precision | Recall | F1 | Precision | Recall | F1
LSTM | 0.7080 | 0.7292 | 0.6863 | 0.6267 | 0.2489 | 0.3563 | 0.7431 | 0.9361 | 0.8285
GRU | 0.7230 | 0.7205 | 0.6400 | 0.7299 | 0.1141 | 0.1974 | 0.7200 | 0.9818 | 0.8308
S-LSTM | 0.7130 | 0.7078 | 0.6042 | 0.7260 | 0.0476 | 0.0894 | 0.7074 | 0.9923 | 0.8260
DistilBERT | 0.4884 | 0.6989 | 0.5750 | 0.0000 | 0.0000 | 0.0000 | 0.6989 | 1.0000 | 0.8227
TABLE III: Results for the Event Detection Task

###  II-A Events Dataset

Table I shows the summary statistics of our event detection dataset. This
dataset comprises a set of tweets and the corresponding list of events. As a
proxy for events label, we utilize a list of all events in Melbourne that had
an event permit issued111https://data.melbourne.vic.gov.au/Events/Event-
permits-2014-2018-including-film-shoots-phot/sex6-6426. As a further pre-
processing, we filtered out events that were likely to not be public, such as
film/movie shoots. Using the Twitter API, we then extracted the tweets posted
during each event and within 100m of the event location, as well as all tweets
in Melbourne from the first event date to the last event date. Following
which, this resulted in a dataset of 1.39 million tweets, of which 22k were
tweets related to an event

###  II-B Crowd Level Dataset

Table II summarizes the main statistics of our crowd level prediction dataset.
There are two components to this dataset, namely a list of pedestrian sensors
and their pedestrian counts and the geo-tagged tweets in the vicinity of these
sensors. The pedestrian sensor component comprises 3.1 million pedestrian
count readings from 66 sensors, retrieved from Melbourne Open
Data222https://data.melbourne.vic.gov.au/Transport/Pedestrian-Counting-
System-2009-to-Present-counts-/b2ak-trbp. In addition, we collected 266k geo-
tagged tweets posted within Melbourne City between 2010 to 2018. These tweets
were then mapped to individual pedestrian sensor locations if their distances
differ by less than 100m.

Upon inspecting the resultant dataset, we realized that many of the tweets had
duplicate coordinates. A likely cause of this issue is the tendency for users
to tag tweets with Twitter “places” [9], which are polygonal areas
corresponding to places of interest that have a standard set of coordinates;
in contrast to having their phone GPS report their accurate coordinates. To
address this issue, we wrote a script to check for duplicate lat/long
coordinates. We set a threshold of 100 maximum duplicates, which left us with
a filtered dataset 10% of its original size.

##  III Event Detection

For the event detection task, one concern is with the imbalance distribution
of classes. To address this issue, we first performed undersampling as our
dataset is very imbalanced with a minority of non-events class [10]. After
this undersampling process, we obtain a dataset with approximately 30% events
and 70% non-events, comprising a total of 74,136 tweets. Using this set of
tweets, we used 90% for training and 5% each for our development and test set.

In our experiments, we empirically evaluate the following neural network based
models, namely: Long Short-Term Memory (LSTM) [11], Gated Recurrent Unit (GRU)
[12], Stacked LSTM (S-LSTM) and DistilBERT [13]. In these recurrent neural
networks, we used GloVe embedding [14] to represent the individual words in
the tweets. We also empirically select the best performing values for the
learning rate, dropout rate, hidden states, batch size on the development set
and report the performance on the test set.

Figure 1: Frequency counts of sensor-based pedestrian counts (left) and tweet
counts by sensor locations (right)

For our evaluation, we use the commonly used metrics of precision, recall and
F1-score. Table III shows the results of our preliminary experiments in terms
of precision, recall and F1-score for each class as well as the weighted mean.
In terms of Weighted F1-score, LSTM is the best performer, followed by GRU,
S-LSTM and DistilBERT. We believe that the smaller training dataset after the
undersampling process contributed to a poorer performance, particularly for
DistilBERT. Moving forward, we intend to further augment this dataset with
more tweets that correspond to the "Event" class and repeat the experiments on
similar datasets.

##  IV Crowd Level Exploratory Analytics

By plotting the tweet and pedestrian count by sensor location in a specific
year, we observe that there is some correlation between the two; where there
are high pedestrian counts near a particular sensor location, the
corresponding number of tweets is likely to also be high. An example is shown
in Figure 1. The exceptions to this trend are largely due to some sensors not
having any linked tweets, especially when there are many sensors close
together (e.g in the city center).

These filtering algorithms left us with a very sparse twitter dataset compared
to our pedestrian count dataset. For instance, taking the counting sensor at
Town Hall (West) in 2017, we had a tweet range from 0-500 per hour and a
pedestrian count range from 200-1000 per hour. Based on a scatter plot of the
tweets and pedestrian hourly counts, we observed a large proportion of data
points where there were no tweets and a varying number of pedestrian counts.
This makes it challenging to find meaningful correlations of the tweet counts
to the pedestrian counts from this scatter plot. This was confirmed by an
initial pass through of the data to a few logistic regression and linear
regression algorithms which showed little to any usable results – the average
accuracy was below 10% for the test set.

A dataset used for crowd level prediction must have reduced sparsity, for
instance by including data from other sources. Flickr is one such source that
we looked at. Flickr is particularly suitable for our purposes, because geo-
location data in Flickr comes from the camera’s metadata. We used the Flickr
API to crawl posts from the same timeframe as the Twitter dataset. Thus, we
can use Flickr descriptions to augment the text data from tweets. However, a
detailed method of extracting meaningful data patterns from a combined data
source requires more investigation.

##  V Conclusion

Our main contribution is in curating two datasets, one for event detection and
one for crowd level prediction. We identify the difficulties of creating a
dataset for crowd level regression – particularly the issues relating to data
sparsity and geo-location accuracy. We also discussed some preliminary results
on the two problems of event detection and crowd level prediction. For the
event detection problem, we experimented with an initial set of neural network
models, such as LSTM, GRU, Staked LSTM and DistilBERT, and discuss preliminary
results and their limitations. For the crowd level prediction problem, we
performed exploratory data analytics to identify trends between social media
data and pedestrian sensors.

##  VI Acknowledgement

This research is funded in part by the Singapore University of Technology and
Design under grant SRG-ISTD-2018-140.

## References

  * [1] K. Al-Kodmany, “Crowd management and urban design: New scientific approaches,” _Urban Design International_ , vol. 18, no. 4, pp. 282–295, 2013.
  * [2] M. Zappatore, A. Longo, and M. A. Bochicchio, “Crowd-sensing our smart cities: a platform for noise monitoring and acoustic urban planning,” _Journal of Communications Software and Systems_ , vol. 13, no. 2, pp. 53–67, 2017.
  * [3] M. Yong and K. H. Lim, “Urban Heat Islands: Beating the Heat with Multi-Modal Spatial Analysis,” in _Proceedings of the 2020 IEEE International Conference on Big Data (BigData’20)_ , Dec 2020.
  * [4] S. B. Ranneries, M. E. Kalør, S. A. Nielsen, L. N. Dalgaard, L. D. Christensen, and N. Kanhabua, “Wisdom of the local crowd: Detecting local events using social media data,” in _Proceedings of the 8th ACM Conference on Web Science_ , 2016, pp. 352–354.
  * [5] M. Fedoryszak, B. Frederick, V. Rajaram, and C. Zhong, “Real-time event detection on social data streams,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 2774–2782.
  * [6] W.-H. Chong and E.-P. Lim, “Fine-grained geolocation of tweets in temporal proximity,” _ACM Transactions on Information Systems (TOIS)_ , vol. 37, no. 2, pp. 1–33, 2019.
  * [7] B. Altaf, L. Yu, and X. Zhang, “Spatio-temporal attention based recurrent neural network for next location prediction,” in _2018 IEEE International Conference on Big Data (Big Data)_. IEEE, 2018.
  * [8] J. Liu, K. L. Wood, and K. H. Lim, “Strategic and Crowd-Aware Itinerary Recommendation,” in _Proceedings of the 2020 European Conference on Machine Learning and Knowledge Discovery in Databases (ECML-PKDD’20)_ , Sep 2020\.
  * [9] F. Botta, H. S. Moat, and T. Preis, “Quantifying crowd size with mobile phone and twitter data,” _Royal Society open science_ , vol. 2, no. 5, p. 150162, 2015.
  * [10] A. Fernández, S. García, M. Galar, R. C. Prati, B. Krawczyk, and F. Herrera, _Learning from imbalanced data sets_. Springer, 2018.
  * [11] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [12] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using rnn encoder–decoder for statistical machine translation,” in _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , 2014, pp. 1724–1734.
  * [13] V. Sanh, L. Debut, J. Chaumond, and T. Wolf, “Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter,” _arXiv preprint arXiv:1910.01108_ , 2019.
  * [14] J. Pennington, R. Socher, and C. D. Manning, “Glove: Global vectors for word representation,” in _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_ , 2014.


# Calculation of isotope shifts and King plot nonlinearities in Ca+

A. V. Viatkina1,2 V. A. Yerokhin1,3 A. Surzhykov1,2 1Physikalisch-Technische
Bundesanstalt, 38116 Braunschweig, Germany  2Technische Universität
Braunschweig, 38106 Braunschweig, Germany  3Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany

###### Abstract

Many-body perturbation theory is implemented in order to calculate the isotope
shifts of $4s$, $4p_{1/2}$, $4p_{3/2}$, $3d_{3/2}$, and $3d_{5/2}$ energy
levels of Ca+, for even isotopes $A=$40, 42, 44, 46, 48. The results are
presented for mass shift and field shift, as well as for higher-order field
shifts, quadratic mass shift, nuclear polarization correction, and the cross
term between field and mass shifts. Additionally, we examine King-plot
nonlinearities introduced by the higher-order isotope-shift corrections to the
combinations of $3d_{3/2}\rightarrow 4s$, $3d_{5/2}\rightarrow 4s$, and
$4p_{1/2}\rightarrow 4s$ transitions. For these transitions, second-order mass
shift and nuclear polarization correction are identified as the dominant
sources of King plot nonlinearity.

## Introduction

High-precision atomic spectroscopy is currently in the spotlight of
experimental research, particularly as a means of testing physics beyond
standard model [1, 2, 3]. The recent advances in experimental technique, most
notably the application of quantum logic spectroscopy [4, 5], allows one to
determine atomic transitions with high accuracy which in certain cases reaches
$0.1$ Hz [6]. One of the precision-spectroscopy methods aimed at detecting new
interactions is King plot analysis; the basic idea was proposed in Ref. [7, 2]
and later refined in Refs. [8, 3, 9]. Roughly speaking, King plot is a plot
constructed from normalized isotope shifts of atomic transitions which are
measured in a series of isotopes of a given element; such plots can be proven
to be linear to a high precision. However, it has been shown that an
interaction between the electron cloud and a hypothetical light boson would
distort this linearity in the majority of cases [2]. Hence, if a nonlinearity
in a King plot is detected in experiment, it could be attributed to a new
boson, given that any other explanation is ruled out [2]. On the other hand,
as one may expect, there are numerous sources of King-plot nonlinearities
found already in the framework of standard model.

Recently, a significant King-plot nonlinearity has been discovered in Yb and
Yb+ transitions [3, 10, 11] and the efforts of its interpretation are still
ongoing, see e.g. Refs. [10, 12]. In contrast, King plots in a succession of
Ca+ isotopes were found to be linear within experimental uncertainties [13,
14]. Isotope shifts were measured for $4s\ ^{2}S_{1/2}\rightarrow 4p\
^{2}P_{1/2}$ and $3d\ ^{2}D_{3/2}\rightarrow 4p\ ^{2}P_{1/2}$ transitions with
accuracy below 100 kHz [15], for $4s\ ^{2}S_{1/2}\rightarrow 3d\ ^{2}D_{5/2}$
with accuracy below 10 Hz [13], and for the interval $3d\ ^{2}D_{3/2}-3d\
^{2}D_{5/2}$ with $\sim 20$ Hz accuracy [14]. One can expect that in the
future the measurement accuracy reaches 1 Hz; in such case, the King-plot
linearity might no longer hold and a careful theoretical examination of the
origins of nonlinearity will be needed. In this contribution, therefore, we
aim to lay ground for an analysis of Ca+ King plot nonlinearity.

It should be noted that, as a system intended for isotope-shift spectroscopy,
Ca+ has several advantages. First, calcium has five even-even stable
($A=40,42,44,46$) or long-lived ($A=48$) isotopes, which is a welcome fact for
building a King plot. Second, the electronic structure and atomic transitions
of Ca+ are theoretically well-understood, given its alkali-like electronic
configuration. Third, Ca+ is a convenient object for experimental study: as
was mentioned above, a series of experiments has already been conducted to
determine several of its atomic transitions with high accuracy [15, 13, 14],
with a realistic possibility for future improvement.

In the present work, we calculate the main contributions to isotope shift in
Ca+ ions and evaluate their impact on King-plot linearity. The paper is
structured as follows. In Sec. I we introduce theoretical origins of isotope
shift terms: we discuss mass shift (Sec. I.1), field shift (Sec. I.2), nuclear
polarization, and field- and mass-shift cross term (Sec. I.3). The total
isotope shift—the sum of all considered contributions—is presented in Sec.
I.4. In Sec. II the specific methods for many-body calculations are laid out,
which we subsequently use in Sec. III to evaluate each of the isotope shift
terms. Finally, we discuss Ca+ King plots and their nonlinearity in Sec. IV
and summarize our results in Sec. V.

##  I Isotope shift: theory

In the first-order approximation, the isotope shift of an energy level in the
$i$-th isotope with respect to a reference isotope $a$ can be written as

$\Delta
E_{ia}=K\left(\frac{m}{M_{i}}-\frac{m}{M_{a}}\right)+F\left(R^{2}_{i}-R^{2}_{a}\right),$
(1)

where $m$ is electron mass and $M_{j}$ denote nuclear masses ($j=a,i$). Here,
we introduce a dimensionless nuclear radius

$R_{j}^{2}=\langle r^{2}\rangle_{j}/\lambdabar^{2}_{C}\ ,$ (2)

which is a mean square nuclear radius $\langle r^{2}\rangle_{j}$ divided by
(the square of) reduced Compton wavelength $\lambdabar_{C}$. The first term in
Eq. (1) corresponds to mass shift (MS) and the second term to field shift
(FS); the former is the consequence of the change in nuclear mass and the
latter of the change in nuclear charge distribution between two isotopes. Note
that in the first-order approximation (1) the electronic structure constants
$K$ and $F$ do not depend on the nuclear parameters of the $i$-th isotope,
while they may still implicitly include the parameters of the reference
isotope $a$. Effectively, Eq. (1) assumes that electronic wavefunctions do not
yet ‘notice’ the change in nuclear mass and shape between isotopes.

In the present work, we investigate the effects beyond the first-order
approximation (1). We take into account that both electronic ‘constants’ $K$
and $F$ depend on the isotope $i$ in question and hence the energy shift can
be written as

$\Delta
E_{ia}=K_{ia}\left(\frac{m}{M_{i}}-\frac{m}{M_{a}}\right)+F_{ia}\left(R^{2}_{i}-R^{2}_{a}\right)\
.$ (3)

Below, we will examine the isotope dependence of $K_{ia}$ and $F_{ia}$,
treating mass shift and field shift separately. Moreover, we will consider
further contributions to isotope shifts which do not, strictly speaking,
belong either to mass or field shift: nuclear polarization correction and
cross term between field and mass shifts. At the end, we will show how to
extend Eq. (3) to take those additional effects into account. Note that, in
this work, we will use 40Ca isotope as our reference isotope $a$ ($A_{a}=40$).

###  I.1 Mass shift

Isotope mass shift arises from the difference in nuclear recoil effect between
two isotopes. For our purposes, it is sufficient to consider this effect
nonrelativistically.111The nonrelativistic treatment allows us to express
normal mass shift in terms of reduced mass. In a relativistic treatment,
however, one has to compute normal mass shift explicitly, which introduces
additional numerical errors. Those errors would obscure the small relativistic
effects present in Ca+ MS. Thus, let us write the Schrödinger Hamiltonian of
the atom:

$H=\frac{\vec{P}^{2}}{2M}+\sum_{k}\frac{\vec{p}_{k}^{\ 2}}{2m}+V_{C}\ .$ (4)

Here $\vec{P}$ refers to the nuclear momentum, $\vec{p}_{k}$ to the momentum
of the $k$-th electron. By $V_{C}$ we denote the Coulomb potential

$V_{C}=-\sum_{k}\frac{Z\alpha}{|\vec{r}_{k}|}+\sum_{k<l}\frac{\alpha}{|\vec{r}_{k}-\vec{r}_{l}|}\
,$ (5)

with $\vec{r}_{k}$ being the position vector of a $k$-th electron. In the
center-of-mass reference frame the nuclear momentum is given by
$\vec{P}=-\sum_{k}\vec{p}_{k}$ and we can rewrite Eq. (4) as

$H=\sum_{k}\frac{\vec{p}_{k}^{\ 2}}{2m_{r}}+V_{C}+\frac{1}{M}V_{\mathrm{SMS}}\
,$ (6)

where $m_{r}=mM/(m+M)$ is the reduced mass. In this equation, the first term
is the normal mass shift (NMS) operator, whose effect can be observed already
in hydrogenlike systems, while the last term $V_{\mathrm{SMS}}$ is the many-
electron specific mass shift (SMS) operator

$V_{\mathrm{SMS}}=\sum_{k<l}\vec{p}_{k}\cdot\vec{p}_{l}\ .$ (7)

To further evaluate (6), we introduce $\mu=m_{r}/m$ and shift the variables in
Eq. (6) as $\vec{r}\rightarrow\mu^{-1}\vec{r}\,$ and, therefore,
$\vec{p}\rightarrow\mu\vec{p}$. Thereby we obtain

$\displaystyle H$ $\displaystyle=\mu\left[\sum_{k}\frac{\vec{p}_{k}^{\
2}}{2m}+V_{C}+\frac{\mu}{M}V_{\mathrm{SMS}}\right]$
$\displaystyle\equiv\mu\left[H_{0}+\frac{\mu}{M}V_{\mathrm{SMS}}\right]\ ,$
(8)

where $H_{0}$ is the nonrelativistic atomic Hamiltonian in the infinite
nuclear mass limit. If we denote the eigenfunctions of $H_{0}$ as $\psi_{0}$,
its eigenvalues as $E_{0}$, and expand the eigenvalues of $H$ in the powers of
$m/M$, we obtain the first- and second-order nuclear recoil (or mass-shift)
corrections to an electronic energy level:

$\displaystyle\delta E_{\mathrm{MS}}^{(1)}$
$\displaystyle=\frac{m}{M}\left(-E_{0}+K^{(1)}_{\mathrm{SMS}}\right)\ ,$ (9)
$\displaystyle\delta E_{\mathrm{MS}}^{(2)}$
$\displaystyle=\left(\frac{m}{M}\right)^{2}\left(E_{0}-K^{(1)}_{\mathrm{SMS}}+K^{(2)}_{\mathrm{SMS}}\right)\
.$ (10)

Here, the first- and second-order specific mass shift constants are given by

$\displaystyle K^{(1)}_{\mathrm{SMS}}$
$\displaystyle=\langle\psi_{0}|V_{\mathrm{SMS}}|\psi_{0}\rangle\ ,$ (11)
$\displaystyle K^{(2)}_{\mathrm{SMS}}$ $\displaystyle={\sum_{n\neq
0}}\frac{|\langle\psi_{0}|V_{\mathrm{SMS}}|\psi_{n}\rangle|^{2}}{E_{0}-E_{n}}\
.$ (12)

By defining the MS constants of the first and second order

$\displaystyle K^{(1)}$ $\displaystyle\equiv-E_{0}+K^{(1)}_{\mathrm{SMS}}\ ,$
(13) $\displaystyle K^{(2)}$ $\displaystyle\equiv
E_{0}-K^{(1)}_{\mathrm{SMS}}+K^{(2)}_{\mathrm{SMS}}\ ,$ (14)

we obtain an expression for the mass shift in an isotope $j$ relative to a
hypothetical infinite-mass isotope:

$\Delta
E_{\mathrm{MS},j}=K^{(1)}\frac{m}{M_{j}}+K^{(2)}\frac{m^{2}}{M^{2}_{j}}\ .$
(15)

Accordingly, the mass shift of an isotope $i$ with respect to the reference
isotope $a$ is

$\displaystyle\Delta$ $\displaystyle E_{\mathrm{MS},ia}=$ $\displaystyle
K^{(1)}\left(\frac{m}{M_{i}}-\frac{m}{M_{a}}\right)+K^{(2)}\left(\frac{m^{2}}{M_{i}^{2}}-\frac{m^{2}}{M_{a}^{2}}\right)\
.$ (16)

It should be noted that when calculating $K^{(1)}$ and $K^{(2)}$ in Sec. III.2
and presenting the results in Table 1 below, we use the experimental values of
binding energies as our $E_{0}$. Strictly speaking, experimental energies $E$
are not equal to the eigenvalues $E_{0}$ of the infinite-mass-isotope
Hamiltonian; however, the difference $|E-E_{0}|$ is negligible for our
purposes, since the uncertainties of our many-body atomic calculations are
much larger.

###  I.2 Field shift

The nuclei of two isotopes differ not only in mass, but also in the parameters
of the nuclear charge distribution. Electronic energy shifts which result from
the difference in nuclear charge distribution—or, more precisely, from the
difference in nuclear potential—are called field shifts (FS). For the purposes
of the present investigation, it is sufficient to assume that all isotopes
have the same shape of the nuclear charge distribution and differ only by the
values of the charge radii. In this approximation, the electrostatic potential
$V$ of a nucleus depends solely on the nuclear radius, $V=V(R)$; its total
isotopic variation would be $\delta V_{ia}=V(R_{i})-V(R_{a})$. In the second-
order perturbation theory, the field shift is given by

$\Delta E_{\mathrm{FS},ia}=\langle\psi_{a}|\delta
V_{ia}|\psi_{a}\rangle+\sum_{n\neq a}\frac{|\langle\psi_{a}|\delta
V_{ia}|\psi_{n}\rangle|^{2}}{E_{a}-E_{n}}\ .$ (17)

Here $\psi_{a}$ is the electronic wave function and $E_{a}$ the electronic
energy level in the reference isotope $a$.

Let $\delta R_{ia}^{2}=R_{i}^{2}-R_{a}^{2}$ be the difference between (squares
of) dimensionless radii. In Eq. (17), we would like to isolate the dominant
first-order field shift $F^{(1)}\delta R_{ia}^{2}$ [see Eq. (1)] from other
field-shift contributions. To do this, we introduce the standard field-shift
operator $V_{\mathrm{FS}}\equiv\partial V(R)/\partial\left(R^{2}\right)$ and
rewrite $\delta V_{ia}$ as

$\delta V_{ia}=V_{\mathrm{FS}}\,\delta R_{ia}^{2}+\left(\frac{\delta
V_{ia}}{\delta R_{ia}^{2}}-V_{\mathrm{FS}}\right)\delta R_{ia}^{2}\ .$ (18)

Now we can express the first term in Eq. (17) as follows:

$\langle\psi_{a}|\delta V_{ia}|\psi_{a}\rangle=F^{(1)}(R_{a})\delta
R_{ia}^{2}+\delta_{R_{i}}F^{(1)}(R_{a})\delta R_{ia}^{2}\ .$ (19)

Here, the first coefficient

$F^{(1)}(R_{a})=\left<\psi_{a}\left|V_{\mathrm{FS}}\right|\psi_{a}\right>,$
(20)

is the standard FS constant in Eq. (1). The remaining part of the field shift
is the higher-order correction

$\delta_{R_{i}}F^{(1)}(R_{a})=\left<\psi_{a}\left|\frac{\delta V_{ia}}{\delta
R_{ia}^{2}}-V_{\mathrm{FS}}\right|\psi_{a}\right>.$ (21)

Since the second term of Eq. (17) is a small correction, one can replace
$\delta V_{ia}$ by $V_{FS}$ in it. Hence, we introduce the second-order field-
shift electronic constant

$F^{(2)}(R_{a})=\sum_{n\neq
a}\frac{|\langle\psi_{a}|V_{\mathrm{FS}}|\psi_{n}\rangle|^{2}}{E_{a}-E_{n}}\
.$ (22)

Finally, we can express the resulting field shift as

$\displaystyle\Delta
E_{\mathrm{FS},ia}=\left[F^{(1)}(R_{a})+\delta_{R_{i}}F^{(1)}(R_{a})\right]$
$\displaystyle\delta R_{ia}^{2}$ $\displaystyle+F^{(2)}(R_{a})$
$\displaystyle\left(\delta R_{ia}^{2}\right)^{2}.$ (23)

Note that the residual potential $\left(\delta V_{ia}/\delta
R^{2}_{ia}-V_{\mathrm{FS}}\right)$ is localized in the nuclear region
similarly to the dominant field-shift potential $V_{\mathrm{FS}}$. Hence, it
makes sense to represent the higher-order correction
$\delta_{R_{i}}F^{(1)}(R_{a})$ as a factor
$f_{\mathrm{ho}}\left(R_{a},R_{i}\right)$ multiplied by the field-shift
constant:

$\delta_{R_{i}}F^{(1)}\left(R_{a}\right)=-F^{(1)}\left(R_{a}\right)f_{\mathrm{ho}}\left(R_{a},R_{i}\right)\cdot
10^{-3}\ .$ (24)

This form is convenient for presenting the results of our numerical
calculation. We will report the numerical results for FS constants
$F^{(1)}\left(R_{a}\right)$, $\delta_{R_{i}}F^{(1)}\left(R_{a}\right)$, and
$F^{(2)}\left(R_{a}\right)$ in Sec. III.3.

###  I.3 Additional terms

####  I.3.1 Nuclear polarization

Not only the nuclear shape, but also the disposition of a nucleus to be
polarized by an electric field changes from isotope to isotope. Nuclear
polarization (np) manifests itself in a correction to electronic energy
levels; varying nuclear polarization between isotopes thus results in a
contribution to isotope shift. Let $\zeta$ be an atomic state and $E_{\zeta}$
its energy; the np-correction to $E_{\zeta}$ would be [16, 17]

$\displaystyle\Delta E_{\mathrm{np}}$ $\displaystyle=$ $\displaystyle-\alpha$
$\displaystyle\sum_{LM}B(EL)\sum_{n}\frac{|\left<\zeta|\mathcal{F}_{L}Y_{LM}|n\right>|^{2}}{E_{n}-E_{\zeta}+\mathrm{sgn}(E_{n})\omega_{L}}\
,$ (25)

where $n$ runs over the complete spectrum of electronic states including
positive and negative continuum, $B(EL)=B(EL;L\rightarrow 0)$ are the reduced
probabilities of nuclear transitions from the excited (“$L$”) to the ground
state (“$0$”), $\omega_{L}$ are the nuclear excitation energies with respect
to the ground state, $Y_{LM}$ the spherical harmonics, and $\mathcal{F}_{L}$
are characteristic radial functions in the sharp-edge-nucleus approximation
[17]:

$\mathcal{F}_{L=0}\left(r\right)=\frac{5\sqrt{\pi}}{2r_{0}^{3}}\left[1-\frac{r^{2}}{r^{2}_{0}}\right]\theta\left(r_{0}-r\right)\
,$ (26) $\displaystyle\mathcal{F}_{L>0}\left(r\right)$
$\displaystyle=\frac{4\pi}{(2L+1)r_{0}^{L}}\left[\frac{r^{L}}{r_{0}^{L+1}}\theta\left(r_{0}-r\right)\right.$
$\displaystyle+\left.\frac{r_{0}^{L}}{r^{L+1}}\theta\left(r-r_{0}\right)\right]\
,$ (27)

where $r_{0}=\sqrt{\langle r\rangle^{2}}$ is the radius of the nuclear sphere
and $\theta$ the step function.

The main contributions to the sum in Eq. (25) arise from two kinds of
transitions: giant resonances and lowest-lying rotational transitions. The
former are dominant because of their large transition strengths $B(EL)$, while
the latter are enhanced due to the small denominator containing the transition
frequency $\omega_{L}$. Accordingly, the dominant low-lying-level contribution
to $\Delta E_{\mathrm{np}}$ in the vast majority of even-even nuclei comes
from the electric quadrupole transition from the first rotational level to the
ground state, $2^{+}\rightarrow 0^{+}$.

To facilitate numerical calculations in many-electron systems, we introduce a
nuclear polarization potential $V_{\mathrm{np}}$ defined by its matrix
elements between single-electron atomic states:

$\displaystyle\langle\zeta|V_{\mathrm{np}}|\xi\rangle=$ $\displaystyle-$
$\displaystyle\alpha\sum_{LM}B(EL)\sum_{n}\frac{\left<\zeta|F_{L}Y_{LM}|n\right>\left<n|F_{L}Y_{LM}|\xi\right>}{E_{n}-m_{e}+\mathrm{sgn}(E_{n})\omega_{L}}\
.$ (28)

Here we replaced $E_{\zeta}$ with electron mass, thereby discarding the (very
weak) dependence of $V_{\mathrm{np}}$ on the binding energy of the state
$\zeta$.

Similarly to the field-shift operator $V_{\mathrm{FS}}$, the operator
$V_{\mathrm{np}}$ is mainly localized in the nuclear region. Hence it is
convenient to present the np correction to a given energy level in Ca+ as

$\Delta E_{\mathrm{np}}\equiv-
g_{\mathrm{np},j}\,R_{j}^{2}\,F^{(1)}\left(R_{a}\right)\cdot 10^{-3}\ ,$ (29)

where $R_{j}$ is the dimensionless nuclear radius (2), $F^{(1)}$ is defined in
Eq. (20), and $g_{\mathrm{np},j}$ is a np coefficient which depends both on
the isotope and on the electronic state in question. Accordingly, the isotope
shift due to nuclear polarization in a given calcium isotope $i$ with regards
to the isotope $a$ would be:

$\Delta
E_{\mathrm{np},ia}=-F^{(1)}\left(R_{a}\right)\left(R_{i}^{2}\,g_{\mathrm{np,}i}-R_{a}^{2}\,g_{\mathrm{np,}a}\,\right)\cdot
10^{-3}\ .$ (30)

The numerical methods for calculating the np coefficients $g_{\mathrm{np}}$
are described in Sec. III.4.2.

####  I.3.2 FS and MS cross term

Finally, let us turn to the isotope shift contribution which is a mixture of
mass shift and field shift. It is convenient to present it as a nuclear-mass-
dependent correction to the field-shift constant. The leading nonrelativistic
effect comes from the reduced mass and could be included into the isotope
shift by the substitution [18] (see also Ref. [19]):

$F\rightarrow\mu^{3}F=F\left(1-3\frac{m}{M}+\dots\right).$ (31)

In calculating the cross term below, we use $F^{(1)}(R_{a})$ in place of the
total $F$, since the difference between the two cases is minuscule in the
already small cross term:

$\Delta F=-3\frac{m}{M}F^{(1)}(R_{a})\ .$ (32)

Additionally, the contributions of the same order induced by the specific
mass-shift operator are expected to be smaller than the reduced-mass effect
and are likewise neglected.

###  I.4 General isotope shift formula

The total isotope shift of an energy level consists of mass shift (Sec. I.1),
field shift (Sec. I.2), and additional terms (Sec. I.3). Let us denote
$F^{(1)}\equiv F^{(1)}(R_{a})$ and $F^{(2)}\equiv F^{(2)}(R_{a})$; written as
a single sum, the main contributions to the isotope shift between the isotopes
$i$ and $a$ are:

$\displaystyle\Delta
E_{ia}=K^{(1)}\left(\frac{m}{M_{i}}-\frac{m}{M_{a}}\right)+K^{(2)}\left(\frac{m^{2}}{M_{i}^{2}}-\frac{m^{2}}{M_{a}^{2}}\right)$
$\displaystyle+\left[1-f_{\mathrm{ho}}\left(R_{a},R_{i}\right)\cdot
10^{-3}-3\left(\frac{m}{M_{i}}-\frac{m}{M_{a}}\right)\right]F^{(1)}\delta
R_{ia}^{2}$
$\displaystyle-F^{(1)}\left(R_{i}^{2}\,g_{\mathrm{np,}i}-R_{a}^{2}\,g_{\mathrm{np,}a}\right)\cdot
10^{-3}+F^{(2)}\left(\delta R_{ia}^{2}\right)^{2}.$ (33)

##  II Many-body perturbation theory

We describe an atom with the relativistic no-pair Dirac-Coulomb Hamiltonian
$H$ which is a sum of the zeroth-order Hamiltonian

$H_{0}=\ \sum_{i}\Big{[}\bm{\alpha}_{i}\cdot{\bm{p}}_{i}+\beta_{i}\,m+V_{\rm
nuc}(r_{i})+U(r_{i})\Big{]}\,,$ (34)

where $V_{\rm nuc}$ is the nuclear binding potential and the residual
electron-electron interaction

$V_{I}=\
\sum_{i<j}\Lambda_{++}\,I(r_{ij})\,\Lambda_{++}-\sum_{i}\Lambda_{+}\,U(r_{i})\,\Lambda_{+}\
.$ (35)

In the latter expression, $\Lambda_{++}$ and $\Lambda_{+}$ are the projection
operators to the positive-energy part of the Dirac spectrum, $I$ is the
electron-electron interaction operator

$I(r_{ij})=\
\frac{\alpha}{r_{ij}}-\frac{\alpha}{2r_{ij}}\Big{[}\bm{\alpha}_{i}\cdot\bm{\alpha}_{j}+(\bm{\alpha}_{i}\cdot\hat{{\bm{r}}}_{ij})(\bm{\alpha}_{j}\cdot\hat{{\bm{r}}}_{ij})\Big{]}\,.$
(36)

Here $\hat{{\bm{r}}}={\bm{r}}/|\bm{r}|$ and $U$ is a screening potential
introduced in the zeroth-order Hamiltonian to partially account for the
electron-electron interaction. An important instance of such a potential is
the Dirac-Hartree-Fock potential $V_{\rm HF}$, whose matrix elements are given
by

$\displaystyle\langle i|V_{\rm HF}|j\rangle\equiv(V_{\rm
HF})_{ij}=\sum_{a}I_{ai;aj}\,,$ (37)

where $I_{ab;cd}\equiv I_{abcd}-I_{abdc}$, $I_{abcd}\equiv\langle
ab|I|cd\rangle$, and $I$ is the operator of the electron-electron interaction
defined in Eq. (36). Here, we will adopt the standard notation from Ref. [20]:
the letters $a$, $b$, $c$, $\ldots$ designate occupied core orbitals; $n$,
$m$, $r$, $\ldots$ signify excited orbitals outside the core, including the
valence orbital; $i$, $j$, $k$, $\ldots$ can be either excited or occupied
orbitals. The letter $v$ stands for the valence orbital.

Within the many-body perturbation theory (MBPT), the energy of the valence
state $E_{v}$ is presented as a perturbation expansion
$E_{v}=E^{(0)}+E^{(1)}+E^{(2)}+\ldots$ When only the contributions to the
valence ionization energy (‘val’) are considered, the expansion terms can be
obtained [20]:

$\displaystyle E^{(0)}_{\rm val}=$ $\displaystyle\ \varepsilon_{v}\ ,$ (38a)
$\displaystyle E^{(1)}_{\rm val}=$ $\displaystyle\big{(}V_{\rm
HF}-U\big{)}_{vv}\,,$ (38b) $\displaystyle E^{(2)}_{\rm val}=$ $\displaystyle\
\sum_{amn}\frac{I_{vamn}\,I_{mn;va}}{\epsilon_{av}-\epsilon_{mn}}-\sum_{abm}\frac{I_{abmv}\,I_{mv;ab}}{\epsilon_{ab}-\epsilon_{vm}}$
$\displaystyle+$ $\displaystyle 2\,\sum_{am}\frac{(V_{\rm
HF}-U)_{am}\,I_{mv;av}}{\varepsilon_{a}-\varepsilon_{m}}$ $\displaystyle+$
$\displaystyle\sum_{i\neq v}\frac{(V_{\rm HF}-U)_{vi}\,(V_{\rm
HF}-U)_{iv}}{\varepsilon_{v}-\varepsilon_{i}}\,,$ (38c)

where $\epsilon_{ab}\equiv\varepsilon_{a}+\varepsilon_{b}$. The expressions
for the third-order MBPT correction $E^{(3)}$ are quite lengthy; they are
presented in Ref. [20] and the angular reduction of these formulae is
described in Ref. [21]. In practical calculations, it is typical to choose the
screening potential $U$ to be the frozen-core Dirac-Fock potential. In such
case, all terms which include the matrix elements of $(V_{\rm HF}-U)$ can be
omitted. However, these terms should be preserved if we plan to perturb the
above formulae with an additional potential; a perturbation of this kind will
be described at the end of this section.

In the present work, we aim to calculate the matrix elements of one-body
(field shift) and two-body (specific mass shift) operators. Moreover, we need
to compute the second-order iterations of such operators. The simplest way to
achieve this is to use the finite-difference approach: the perturbing
operators are first added to the Hamiltonian and then the numerical derivative
with respect to the perturbations is evaluated. Specifically, for calculating
the field-shift constants $F^{(1)}$ and $F^{(2)}$, we add the perturbing
potential $V_{\rm FS}$ with an arbitrary pre-factor $\lambda$ to the nuclear
potential, i.e., $V_{\rm nuc}\to V_{\rm nuc}+\lambda\,V_{\rm FS}$. Then the
valence energies $E_{v}(\lambda)$ of the modified Hamiltonian are calculated.
Finally, the field-shift constants are obtained by computing numerical
derivatives with respect to the parameter $\lambda$:

$F^{(1)}=\left.\frac{\partial
E_{v}(\lambda)}{\partial\lambda}\right|_{\lambda=0}\,,\ \ \
F^{(2)}=\left.\frac{1}{2}\frac{\partial^{2}E_{v}(\lambda)}{\partial\lambda^{2}}\right|_{\lambda=0}\,.$
(39)

Similarly, the specific mass-shift constants $K^{(1)}_{\rm SMS}$ and
$K^{(2)}_{\rm SMS}$ are calculated by adding the perturbing potential $V_{\rm
SMS}$ to the electron-electron interaction, $I(r_{ij})\to
I(r_{ij})+\beta\,V_{\rm SMS}(r_{ij})$ and finding the derivative of the
perturbed valence energies $E_{v}(\beta)$ with respect to $\beta$,

$\displaystyle K^{(1)}=\left.\frac{\partial
E_{v}(\beta)}{\partial\beta}\right|_{\beta=0}\,,\ \ \
K^{(2)}=\left.\frac{1}{2}\frac{\partial^{2}E_{v}(\beta)}{\partial\beta^{2}}\right|_{\beta=0}\,.$
(40)

The parameters $\lambda$ and $\beta$ are chosen in such a way that the
resulting changes in energy are much larger than the round-off errors in the
numerical calculation and, at the same time, sufficiently small for the
numerical derivative to be stable against variations of $\lambda$ and $\beta$.

Table 1: Mass shift and field shift constants in Ca+, see Eq. (33) and Eqs. (13–14); the respective units are described in Sec. III.1, Eqs. (46a–46d). Our results for $\tilde{K}^{(1)}_{\mathrm{SMS}}$ and $\tilde{F}^{(1)}(R_{a})$ are compared with the SMS and FS constants calculated in Refs. [22, 23]. The dimensionless nuclear radius $R_{a}$ [see Eq. (2)] belongs to the 40Ca isotope. The overall sign of the field shift constants from Refs. [22, 23] was reversed to conform to the definition used in the present work.  | Units |  | $4s$ | $4p_{1/2}$ | $4p_{3/2}$ | $3d_{3/2}$ | $3d_{5/2}$
---|---|---|---|---|---|---|---
$\tilde{K}^{(1)}$ | GHz$\times$amu |  | 1324 | 940 | 951 |  $-$1136 |  $-$1124
$\tilde{K}^{(1)}_{\rm SMS}$ | GHz$\times$amu |  |  $-$251 |  $-$221 |  $-$206 |  $-$2487 |  $-$2474
|  | MBPT+RPA [22] |  $-$259 |  $-$204 |  $-$200 |  $-$2601 |  $-$2595
|  | RCCSD(T) [23] |  $-$243 |  $-$208 |  $-$204 |  $-$2364 |  $-$2357
$\tilde{K}^{(2)}$ | GHz$\times$amu2 |  |  $-$3.26 |  $-$8.79 |  $-$8.78 |  $-$3.87 |  $-$3.85
$\tilde{K}^{(2)}_{\rm SMS}$ | GHz$\times$amu2 |  |  $-$2.53 |  $-$8.27 |  $-$8.26 |  $-$4.49 |  $-$4.46
$\tilde{F}^{(1)}(R_{a})$ | MHz$\times$fm-2 |  | 266.8 |  $-$19.6 |  $-$19.9 |  $-$112.2 |  $-$111.6
|  | MBPT+RPA [22] | 266.6 |  $-$19.6 |  $-$19.9 |  $-$111.8 |  $-$111.2
|  | RCCSD(T) [23] | 263.3 |  $-$15.9 |  $-$19.5 |  $-$93.9 |  $-$112.3
$\tilde{F}^{(2)}(R_{a})$ | kHz$\times$fm-4 |  |  $-$89.25 | 6.58 | 6.65 | 37.50 | 37.30

In our calculations, we choose the potential $U$ in the zeroth-order $H_{0}$
to be the Dirac-Fock potential. When a perturbation is added to the
Hamiltonian $H_{0}$ in the finite-difference approach, one can include this
perturbation into the self-consistent procedure of computing the Dirac-Fock
potential. It was demonstrated in Ref. [24] that such an inclusion is
advantageous because it accounts for the infinite sequence of diagrams known
as the random-phase-approximation (RPA) corrections, thus yielding a
substantial improvement in the accuracy of calculations. We adopt this
approach and ensure the self-consistency of the Dirac-Fock potential after the
perturbation is added to $H_{0}$, for each value of the parameters $\lambda$
and $\beta$.

Apart from the finite-difference method, we also implement a more traditional
perturbative approach for calculating the first-order matrix elements of one-
body operators—specifically, the nuclear polarization correction and the
higher-order finite nuclear size correction. Let $V$ be the one-body
potential; we consider the linear-in-$V$ perturbations of the MBPT formulas
for energy levels. Accordingly, the perturbations of single-electron energies
and wave functions are

$\displaystyle\varepsilon_{i}\to\langle i|V|i\rangle\equiv V_{ii}\,,\ \ \
|i\rangle\to|\delta i\rangle=\sum_{k\neq
i}\frac{|k\rangle\,V_{ki}}{\varepsilon_{i}-\varepsilon_{k}}\,.$ (41)

The perturbations of the first two contributions in Eq. (38a) yield

$\displaystyle\delta E^{(0)}_{V}=$ $\displaystyle\ V_{vv}\,,$ (42)
$\displaystyle\delta E^{(1)}_{V}=$ $\displaystyle\ \sum_{k\neq
v}\frac{\big{(}V_{\rm
HF}-U\big{)}_{vk}\,V_{kv}}{\varepsilon_{v}-\varepsilon_{k}}+2\,\sum_{a}I_{\delta\\!av;av}\,.$
(43)

Formulas for the next-order correction $\delta E^{(2)}_{V}$ are easily
obtained from Eq. (38c) by perturbing all single-electron energies and wave
functions according to Eq. (41) and preserving only the part which is linear
in $V$.

In order to obtain results equivalent to those delivered by the finite-
difference approach, we need to include the contributions from RPA explicitly.
This can be achieved by defining the RPA-corrected single-particle matrix
elements as [22]

$\displaystyle V_{an}^{\,\rm RPA}=V_{an}+\sum_{bm}\frac{V_{bm}^{\,\rm
RPA}\,I_{am;nb}+I_{ab;nm}\,V_{mb}^{\,\rm
RPA}}{\varepsilon_{b}-\varepsilon_{m}}\,,$ (44)

and solving the equation iteratively, obtaining the ‘dressed’ single-particle
matrix elements $V_{an}^{\,\rm RPA}$. Then, in the calculation of $\delta
E^{(1)}_{V}$, the resulting $V_{an}^{\,\rm RPA}$ are used instead of $V_{an}$.
One can demonstrate that the first-order RPA iteration accounts for a part of
$\delta E^{(2)}_{V}$, while the higher-order RPA iterations correspond to an
infinite sequence of higher-order diagrams.

Furthermore, the dressed RPA matrix elements can be implemented in
calculations of $\delta E^{(2)}_{V}$, thus accounting for additional sequence
of higher-order contributions. Of course, the double-counting terms should be
removed from the RPA-corrected $\delta E^{(1)}_{V}$ and $\delta E^{(2)}_{V}$
contributions. A very similar approach was used in Ref. [22] for the first-
order field-shift correction.

##  III Calculation of isotope shift parameters

###  III.1 Units

Before we consider the specifics of Ca+ numerical calculations, let us briefly
discuss the units of isotope shift parameters. According to our definition in
Eq. (3), the parameters $K_{ia}$ and $F_{ia}$ have units of energy, since both
$m/M_{j}$ and $R^{2}_{j}$ are dimensionless. Such a definition is
indispensable for theoretical treatment, e.g., when we consider the small-
parameter expansion in Sec. I.1. In the literature, however, the mass shift
constant is expressed in GHz$\times$amu and the field shift constant in
MHz$\times$fm-2. In this case, Eq. (3) would have the following form

$\Delta
E_{ia}=\tilde{K}_{ia}\left(\frac{1}{M_{i}}-\frac{1}{M_{a}}\right)+\tilde{F}_{ia}\delta\langle
r^{2}\rangle_{ia}\ ,$ (45)

The connection between the two definitions of the isotope-shift constants is
as follows:

$\displaystyle\hbar\,\tilde{K}^{(1)}$ $\displaystyle=mK^{(1)}\ ,$ (46a)
$\displaystyle\hbar\,\tilde{K}^{(2)}$ $\displaystyle=m^{2}K^{(2)}\ ,$ (46b)
$\displaystyle\hbar\,\tilde{F}^{(1)}$
$\displaystyle=\lambdabar_{c}^{-2}F^{(1)}\ ,$ (46c)
$\displaystyle\hbar\,\tilde{F}^{(2)}$
$\displaystyle=\lambdabar_{c}^{-4}F^{(2)}\ .$ (46d)

In order to make our results more accessible for calculating isotope shifts,
in Table 1 we present the numerical results in the form of $\tilde{K}^{(1)}$,
$\tilde{K}^{(2)}$, $\tilde{F}^{(1)}$ and $\tilde{F}^{(2)}$.

It should be mentioned that, in all calculations reported in this work, we do
not include the contributions induced solely by the core electrons: they are
the same for all states investigated here and cancel when transitions between
two states are considered.

Table 2: Higher-order field-shift and nuclear-polarization contributions for the ground and first excited states of Ca+ evaluated with different methods. The methods (i), (ii), and (iii) are described in Sec. III.4.1 and III.4.2. Higher-order field-shift constant correction $\delta_{R_{46}}F^{(1)}=\delta_{R_{46}}F^{(1)}(R_{40})$ is evaluated for the nuclear radii $r_{40}=3.4776$ fm and $r_{46}=3.4953$ fm [26]; note that $R_{j}=r_{j}/\lambdabar_{C}$. Nuclear polarization correction $\Delta E_{\rm np}$ to the electronic energy levels is found for the 40Ca+ isotope. Units for $\delta_{R_{46}}F^{(1)}$ are kHz$\times$fm-2, while $\Delta E_{\rm np}$ is given in MHz.  | $4s$ | $4p_{1/2}$ | $4p_{3/2}$ | $3d_{3/2}$ | $3d_{5/2}$
---|---|---|---|---|---
Method | $\delta_{R_{46}}F^{(1)}$ | $f_{\rm ho}$ | $\delta_{R_{46}}F^{(1)}$ | $f_{\rm ho}$ | $\delta_{R_{46}}F^{(1)}$ | $f_{\rm ho}$ | $\delta_{R_{46}}F^{(1)}$ | $f_{\rm ho}$ | $\delta_{R_{46}}F^{(1)}$ | $f_{\rm ho}$
(i) | -4.810 | 0.0238 | 0.642 | 0.0238 | 0.645 | 0.0237 | 3.814 | 0.0238 | 3.794 | 0.0238
(ii) | -5.354 | 0.0238 | 0.452 | 0.0238 | 0.456 | 0.0237 | 2.011 | 0.0238 | 1.998 | 0.0238
(iii) | -6.202 | 0.0238 | 0.447 | 0.0238 | 0.453 | 0.0237 | 2.519 | 0.0238 | 2.506 | 0.0238
| $\Delta E_{\rm np}$ | $g_{\rm np}$ | $\Delta E_{\rm np}$ | $g_{\rm np}$ | $\Delta E_{\rm np}$ | $g_{\rm np}$ | $\Delta E_{\rm np}$ | $g_{\rm np}$ | $\Delta E_{\rm np}$ | $g_{\rm np}$
(i) | -0.891 | 0.364 | 0.119 | 0.365 | 0.119 | 0.363 | 0.709 | 0.366 | 0.705 | 0.366
(ii) | -0.991 | 0.364 | 0.084 | 0.365 | 0.084 | 0.363 | 0.374 | 0.366 | 0.372 | 0.366
(iii) | -1.148 | 0.364 | 0.083 | 0.365 | 0.084 | 0.363 | 0.469 | 0.366 | 0.466 | 0.366
Table 3: The higher-order field shift parameter $f_{\mathrm{ho}}\left(R_{a},R_{i}\right)$ [see Eq. (24), $A_{a}=40$] for atomic states of Ca+.  $A_{i}$ | $f_{\mathrm{ho}}\left(R_{a},R_{i}\right)$
---|---
| $4s$ | $4p_{1/2}$ | $4p_{3/2}$ | $3d_{3/2}$ | $3d_{5/2}$
42 | 0.0410 | 0.0411 | 0.0410 | 0.0410 | 0.0410
44 | 0.0542 | 0.0543 | 0.0541 | 0.0541 | 0.0541
46 | 0.0238 | 0.0238 | 0.0238 | 0.0238 | 0.0238
48 |  $-$0.0007 |  $-$0.0007 |  $-$0.0007 |  $-$0.0007 |  $-$0.0007
Table 4: Nuclear parameters and the nuclear polarization correction coefficient $g_{\rm np}$ [see Eqs. (29–30)] for various isotopes of Ca+. The energy of the first nuclear rotational state is $\omega$, $B(E2){\uparrow}$ is the corresponding nuclear transition strength, and $\sqrt{\langle r^{2}\rangle}$ is the root-mean-square nuclear charge radius.  $A$ |  $\omega$ [25] |  $B(E2){\uparrow}$ [25] |  $\sqrt{\langle r^{2}\rangle}$ [26] | $g_{\rm np}$
---|---|---|---|---
| [MeV] | $[e^{2}{\rm b}^{2}]$ | [fm] | $4s$ | $4p_{1/2}$ | $4p_{3/2}$ | $3d_{3/2}$ | $3d_{5/2}$
40 | 3.904 | 0.010 | 3.4776 | 0.3643 | 0.3646 | 0.3628 | 0.3658 | 0.3660
42 | 1.525 | 0.042 | 3.5081 | 0.3961 | 0.3965 | 0.3945 | 0.3977 | 0.3979
44 | 1.157 | 0.047 | 3.5179 | 0.4122 | 0.4126 | 0.4106 | 0.4139 | 0.4141
46 | 1.346 | 0.018 | 3.4953 | 0.4120 | 0.4123 | 0.4102 | 0.4138 | 0.4140
48 | 3.832 | 0.010 | 3.4771 | 0.4242 | 0.4245 | 0.4224 | 0.4262 | 0.4263

###  III.2 Mass shift

Our calculations of the first- and second-order specific-mass-shift constants
$K^{(1)}_{\rm SMS}$ and $K^{(2)}_{\rm SMS}$ were performed with the finite-
difference approach described in Sec. II which reduces the problem at hand to
the computation of energies. We included all MBPT corrections to energy up to
the third order, i.e., up to the three-photon exchange. This approach accounts
for both the two-photon-exchange corrections to the matrix element of the SMS
operator and an infinite sequence of higher-order corrections delivered by the
random-phase approximation. The numerical results for the mass shift constants
are presented in Table 1. Our values for the first-order constant are in a
reasonable agreement with previous calculations [22, 23]. For the second-order
constants, there have been no results reported in the literature.

###  III.3 Field shift

The calculations of the first- and second-order field-shift constants
$F^{(1)}$ and $F^{(2)}$ were performed with the finite-difference approach as
described in Sec. II. We included MBPT corrections to energy up to the second
order, which accounts for all two-photon-exchange corrections to the matrix
element of the field shift operator.

The numerical results for $F^{(1)}$ and $F^{(2)}$ are presented in Table 1.
Our values for the first-order constant $F^{(1)}$ are in a very good agreement
with results by Safronova and Johnson [22]. There is, however, a significant
difference with results of Ref. [23], especially for $4p_{1/2}$ and $3d_{3/2}$
states. The probable reason could be that the numerical approach of Ref. [23]
is not particularly suitable for computation of the field shift.

###  III.4 Additional terms

####  III.4.1 Higher order field shift

In the present work, we consider three different approaches for evaluating the
field-shift correction $\delta_{R_{i}}F^{(1)}\left(R_{a}\right)$ [see Eq.
(21)]. In the order of the improvement of approximation, they are: (i) Dirac-
Fock including core relaxation, based on Eqs. (42) and (43); (ii) random phase
approximation (RPA) which relies on adding Eq. (44) to the method (i); and,
finally, (iii) MBPT2+RPA, in which we added two-photon exchange correction to
the method (ii). Table 2 compares the results for
$\delta_{R_{46}}F^{(1)}(R_{40})$ in Ca+ obtained with different methods. In
this table, together with the $\delta_{R_{46}}F^{(1)}(R_{40})$ correction, we
display the factor $f_{\rm ho}$ which was introduced in Eq. (24). We conclude
that, to the level of $10^{-3}$ relative accuracy, the factor $f_{\rm ho}$ is
both method-independent and the same for all states considered. One can also
note that the results obtained are very close to the hydrogenic $1s$ values.
For instance, the hydrogenic result is $f_{\rm ho}(R_{46},R_{40})=0.02378$ for
$A_{i}=46$ in Ca+. Our final results for the higher-order field-shift
correction are summarized in Table 3.

####  III.4.2 Nuclear polarization

The nuclear polarization correction is induced by the operator $V_{\rm np}$
defined in Eq. (28). The single-electron matrix elements of $V_{\rm np}$ were
calculated in the same way as in our previous work [27]. Specifically, we
included the dominant $E2$ nuclear rotational transition and the giant
resonance transitions with $L=0,1,2,3$. The sum over the Dirac spectrum was
performed using the finite-basis-set $B$-spline method [28]. For hydrogenic
matrix elements we reproduce the results of Ref. [17]. A similar calculation
has been recently presented in Ref. [29]; the difference is that approximate
empirical formulas for $B(E2)$ were used and only the dominant $L=1$ giant
resonance was included. Based on Eq. (29), the np correction is conveniently
expressed in terms of a $g_{\rm np}$ coefficient which is defined by the ratio
of nuclear-polarization and field-shift contributions. While calculating this
ratio, it is important to use the same method for both contributions. In such
case, the ratio will not depend on the method of accounting for electron
correlations.

The electron-structure corrections to the np effect are calculated with the
same approaches as the higher-order field shift in Sec. III.4.1: (i) Dirac-
Fock with core relaxation, (ii) RPA, and (iii) MBPT2$+$RPA. The numerical
results for both quantities $\Delta E_{\mathrm{np}}$ and $g_{\rm np}$
delivered by these three methods for the singly-ionized calcium ion 40Ca+ are
listed in Table 2. We find that the numerical results expressed in terms of
$g_{\rm np}$ do not depend on the method of calculation and only weakly depend
on the electronic valence state. This is not surprising since it is known that
for light atoms both the nuclear polarization and the finite nuclear size
corrections are roughly proportional to the expectation value of the Dirac
$\delta$-function [30].

We also find that the results for $g_{\rm np}$ in 40Ca+ from Table 2 are very
close to the hydrogenic $1s$ value $g_{\rm np}(1s)=0.360$. The deviation is
within 2%, which is much smaller than the uncertainty associated with the
approximate treatment of the nuclear polarization effect. Therefore, for many
practical purposes, it is sufficient just to use the hydrogenic values of
$g_{\rm np}$.

In Table 4 we present our results for the nuclear polarization contribution in
various isotopes of Ca+. The experimental values of the nuclear quadrupole
transition strengths $B(E2)$ and the excitation energies $\omega_{L}$
originate from Ref. [25], while nuclear charge radii are from Ref. [26].

###  III.5 Numerical results

Our numerical results for the isotope shifts of Ca+ isotopes with $A=42$, 44,
46, and 48, relative to the reference isotope $A_{0}=40$ are summarized in
Table 5. Note that the numerical results listed in the table do not include a
contribution of the electron core that is the same for all states considered
and cancels for the transition energies.

The uncertainty of our theoretical values of the mass-shift and field-shift
constants stems predominantly from the higher-order electron-correlation
effects, because all computational errors (basis-set truncation, errors of
numerical differentiation) are small. There is no safe way to estimate the
omitted electron-correlation effects. To some extent, this can be done by
comparing calculations performed by different methods. The specific mass shift
is known to be particularly difficult to calculate reliably for many-electron
atoms. The agreement of results obtained by different method summarized in
Table 1 for $K^{(1)}_{\mathrm{SMS}}$ is on the level of 10%, which could be
taken as an estimation of the uncertainty. Calculations of the field shift can
generally be performed to a higher accuracy than for the mass shift. We assume
that the uncertainty of our results for $F^{(1)}$ should be within 5%.

There are no independent calculations to compare with for the second-order
isotope shift constants. We assume here that their uncertainties should be
comparable with those of the first-order constants $K^{(1)}_{\mathrm{SMS}}$
and $F^{(1)}$. Furthermore, there are uncertainties originating from the
nuclear model employed in the calculation of the nuclear polarization effect.
We assume these uncertainties to be on the level of 10%.

Table 5: Isotope-shift contributions in MHz and kHz, according to Eq. (33), with reference isotope $A_{0}=40$. Note that we do not include the large contributions of core electrons to isotope shifts which would be the same for every single-electron state and cancel out when the difference is considered. MS and FS stand for mass shift and field shift, respectively; ‘sec.’ indicates ‘second order’ and ‘h.o.’ means ‘higher order’. Nuclear polarization contribution is denoted by ‘np’, while ‘cross term’ means field- and mass-shift cross term. $A$ | isotope shift | units | $4s$ | $4p_{1/2}$ | $4p_{3/2}$ | $3d_{3/2}$ | $3d_{5/2}$
---|---|---|---|---|---|---|---
$42$ | total | MHz |  $-$1519.6 |  $-$1122.9 |  $-$1136.5 | 1329.8 | 1315.6
| MS | MHz |  $-$1576.5 |  $-$1119.2 |  $-$1132.8 | 1353.4 | 1339.1
| MS, sec. | MHz | 0.190 | 0.511 | 0.511 | 0.225 | 0.224
| FS | MHz | 56.8 |  $-$4.18 |  $-$4.24 |  $-$23.9 |  $-$23.8
| FS, h.o. | kHz |  $-$2.33 | 0.172 | 0.174 | 0.981 | 0.975
| FS, sec. | kHz |  $-$4.05 | 0.299 | 0.302 | 1.70 | 1.69
| np | kHz |  $-$125.1 | 9.24 | 9.31 | 52.8 | 52.5
| cross term | kHz | 0.111 |  $-$0.008 |  $-$0.008 |  $-$0.047 |  $-$0.047
$44$ | total | MHz |  $-$2935.1 |  $-$2141.7 |  $-$2167.7 | 2553.2 | 2526.0
| MS | MHz |  $-$3010.4 |  $-$2137.1 |  $-$2163.1 | 2584.3 | 2557.0
| MS, sec. | MHz | 0.354 | 0.955 | 0.954 | 0.420 | 0.418
| FS | MHz | 75.2 |  $-$5.54 |  $-$5.62 |  $-$31.6 |  $-$31.5
| FS, h.o. | kHz |  $-$4.08 | 0.301 | 0.304 | 1.71 | 1.70
| FS, sec. | kHz |  $-$7.09 | 0.523 | 0.529 | 2.98 | 2.96
| np | kHz |  $-$185.5 | 13.7 | 13.8 | 78.4 | 78.0
| cross term | kHz | 0.282 |  $-$0.021 |  $-$0.021 |  $-$0.118 |  $-$0.118
$46$ | total | MHz |  $-$4287.2 |  $-$3068.2 |  $-$3105.5 | 3695.8 | 3656.6
| MS | MHz |  $-$4320.5 |  $-$3067.1 |  $-$3104.4 | 3709.0 | 3669.7
| MS, sec. | MHz | 0.498 | 1.34 | 1.34 | 0.590 | 0.587
| FS | MHz | 32.9 |  $-$2.42 |  $-$2.46 |  $-$13.9 |  $-$13.8
| FS, h.o. | kHz |  $-$0.784 | 0.058 | 0.059 | 0.330 | 0.328
| FS, sec. | kHz |  $-$1.36 | 0.100 | 0.101 | 0.571 | 0.568
| np | kHz |  $-$167.5 | 12.3 | 12.4 | 70.9 | 70.5
| cross term | kHz | 0.177 |  $-$0.013 |  $-$0.013 |  $-$0.074 |  $-$0.074
$48$ | total | MHz |  $-$5528.5 |  $-$3922.6 |  $-$3970.3 | 4746.9 | 4696.6
| MS | MHz |  $-$5528.0 |  $-$3924.4 |  $-$3972.1 | 4745.6 | 4695.4
| MS, sec. | MHz | 0.624 | 1.68 | 1.68 | 0.741 | 0.737
| FS | MHz |  $-$0.928 | 0.068 | 0.069 | 0.390 | 0.388
| FS, h.o. | kHz |  $-$0.00065 | 0.00005 | 0.00005 | 0.00027 | 0.00027
| FS, sec. | kHz |  $-$0.0011 | 0.0001 | 0.0001 | 0.0005 | 0.0005
| np | kHz |  $-$192.9 | 14.2 | 14.3 | 81.8 | 81.2
| cross term | kHz |  $-$0.0064 | 0.0005 | 0.0005 | 0.0027 | 0.0027

##  IV King plots

###  IV.1 Theoretical introduction

The method of King plots [31] is a popular way to analyze the experimentally
measured isotope shifts. There are different ways to construct King plots from
spectroscopy data to serve similar purposes, see Refs. [31, 3]; we will
consider the most widely used version. The core idea of a King plot is that,
to the leading order, isotope shifts depend linearly on $m/M_{i}$ and
$R^{2}_{i}$, see Eq. (1). Hence, if one considers isotope shifts of two
different transitions, it is possible to eliminate the poorly known
$R^{2}_{i}$ from a system of linear equations. To illustrate this approach,
let $\Delta\nu=\Delta E^{(e)}-\Delta E^{(g)}$ be a transition frequency
between an excited and ground atomic energy levels. Then, to the first order

$\displaystyle\begin{cases}\Delta\nu_{1,ia}=K_{\nu_{1}}\mu_{ia}+F_{\nu_{1}}\delta
R^{2}_{ia}\ ,\\\ \Delta\nu_{2,ia}=K_{\nu_{2}}\mu_{ia}+F_{\nu_{2}}\delta
R^{2}_{ia}\ ,\end{cases}$ (47)

with $\mu_{ia}=\left(m/M_{i}-m/M_{a}\right)$. Consequently, we can write a
linear relation

$n_{2,ia}=\frac{F_{\nu_{2}}}{F_{\nu_{1}}}\
n_{1,ia}+\left(K_{\nu_{2}}-\frac{F_{\nu_{2}}}{F_{\nu_{2}}}K_{\nu_{1}}\right)\
,$ (48)

where $n_{k,ia}=\Delta\nu_{k,ia}/\mu_{ia}$ are modified frequencies ($k=1,2$).
The plot of $n_{2,ia}$ against $n_{1,ia}$ for different isotope pairs $(i,a)$
is called King plot (see Fig. 1) and, to the first order in isotope shift, it
is linear.

The higher-order corrections described in Sec. I and summarized in Eq. (33)
distort the linear relation (48). To quantify this deviation from linearity,
we will apply two different methods shown schematically in Fig. 1. The first
method relies on averaging the difference between the outlier point and the
line defined by two other points:

$\Delta_{\mathrm{line}}=(\Delta_{1}+\Delta_{2})/2,\quad\Delta_{1,2}>0\ .$ (49)

To scale $\Delta_{\mathrm{line}}$ back to frequency units, we multiply its
value by $|m/M_{j}-m/M_{a}|$, where ($j,a$) is the isotope pair of the outlier
point.

Figure 1: Two ways to evaluate the nonlinearity of a King plot. Here $n_{1,2}$
are modified isotope shift frequencies and $(i,a)$ indicates the isotope pair
which corresponds to the given point; $i=1,2,3$ and we assume $a$ to be the
reference isotope 40Ca. Method a) characterizes the average difference
$\Delta_{\mathrm{line}}=\left(\Delta_{1}+\Delta_{2}\right)/2$ between the
outlier point and the line defined by the two other points. Method b)
quantifies the area of the triangle defined by the three points, see Eq. (51)
and Ref. [2].

The second method considers the area of a triangle defined by three points of
a King plot [2]. To formalize this approach, we need to define ‘isotope shift
vectors’ whose coordinates correspond to modified frequency shifts in
different isotope pairs:

$\vec{n}_{k}=(n_{k,1a},\ n_{k,2a},\ n_{k,3a}).$ (50)

The area of the King-plot triangle can be calculated as [2, 8]

$\Delta_{V_{2}}=\frac{1}{2}\,\mathrm{det}\left(\vec{n}_{1},\ \vec{n}_{2},\
\vec{\mathbb{1}}_{3}\right),$ (51)

where $\vec{\mathbb{1}}_{3}\equiv\left(1,1,1\right)$ is a unit vector. Again,
in order to scale $\Delta_{\mathrm{V_{2}}}$ to frequency units, we multiply it
by the second power of the mean value $\overline{|m/M_{i}-m/M_{a}|}^{2}$,
where ($i,a$) are all four isotope pairs.

The latter approach is easily generalized to include more isotopes and more
transitions [8]. For instance, the isotope shifts in Ca+ can be measured in
four pairs of isotopes: the reference $A=$40 and each of the $A=$42, 44, 46,
48 isotopes. In this way, we would obtain four-dimensional isotope-shift
vectors:

$\vec{n}_{k}=\left(n_{k,1a},\ n_{k,2a},\ n_{k,3a},\ n_{k,4a}\right).\\\ $ (52)

The nonlinearity is defined by

$\Delta_{V_{3}}=\frac{1}{6}\,\mathrm{det}\left(\vec{n}_{1},\ \vec{n}_{2},\
\vec{n}_{3},\ \vec{\mathbb{1}}_{4}\right),$ (53)

where $\vec{\mathbb{1}}_{4}=\left(1,1,1,1\right)$ is a unit vector. We scale
$\Delta_{\mathrm{V_{3}}}$ back to frequency units by introducing the mean-
value factor $\overline{|m/M_{i}-m/M_{a}|}^{3}$, where ($i,a$) are all four
isotope pairs. This case corresponds to a three-dimensional King plot.

In the present work we analyze nonlinearities both for the two- and three-
dimentional King plots. In former case, we consider two sets of isotopes:
$(A_{i},40)$ with $A_{i}=$42, 44, 46 and $A_{i}=$42, 44, 48. Note that the
charge radii of 40Ca and 48Ca are almost equal, while the radii of the $A=$42,
44, 46 isotopes are larger than both, as can be seen from Table 4; this fact
may influence the pattern of King plot nonlinearity in a way that is specific
for calcium.

###  IV.2 King-plot nonlinearity in Ca+

We calculate two-dimensional King-plot nonlinearities (NLs) with two methods
illustrated in Fig. 1, for two pairs of transitions: ($3d_{3/2}\rightarrow
4s$; $3d_{5/2}\rightarrow 4s$) and ($4p_{1/2}\rightarrow 4s$;
$3d_{3/2}\rightarrow 4s$). The first pair consists of two narrow transitions
and is suitable for the search of possible new-physics effects [14]. However,
this transition pair poses a challenge for theoretical predictions of NLs, as
it involves significant cancellations between the isotope shifts. On the other
hand, the second pair of transitions does not have such strong cancellations
and, consequently, is better suited for detecting King-plot nonlinearities
originating from the standard model, as will be discussed below.

Since there are four isotope pairs available for Ca+—four points on a 2D King
plot—we calculate NLs both for the ‘third point’ ($A=$42, 44, 46 vs. $A=$40)
and the ‘fourth point’ ($A=$42, 44, 48 vs. $A=$40). Additionally, we determine
the 3D King plot NL as defined in Eq. (53), for the following set of
transitions: $3d_{3/2}\rightarrow 4s$, $3d_{5/2}\rightarrow 4s$, and
$4p_{1/2}\rightarrow 4s$.

Our numerical results are presented in Tables 6 and 7 for the NL of the 2D and
3D King plots, respectively. Individual contributions to 2D NL are roughly
additive; in Table 6 we present the total NL results as well as the individual
contributions expressed in percentage of the total values. On the other hand,
individual contributions to the 3D NL turn out to be non-additive. So, we
successively exclude each of the terms and list the resulting change of NL
values in Table 7.

Examining individual contributions to NL, we observe that the second-order and
higher-order FS corrections and the cross term are essentially negligible for
Ca+ and can be omitted in future studies. This was expected, since the
nuclear-size effects tend to be small for light atoms. The main source of NL
in Ca+ is the second-order mass shift, which is consistent with our previous
findings for argon [18]. Quite surprisingly, however, we find a large
contribution from the nuclear polarization for the 4th-point
($3d_{3/2},3d_{5/2}\rightarrow 4s$) 2D King-plot and for the 3D King plot. A
possible reason is the irregular behavior of the nuclear radii: the charge
radii of 40Ca and 48Ca are almost equal, whereas the radii of 42,44,46Ca are
larger than both of them (see Table 4).

One might notice that the methods a) and b) produce NLs that are, first, of
different units and, second, very different numerically. This is not
surprising since the method a) measures a linear distance whereas the method
b) measures the area of a triangle. Naturally, both methods describe the same
phenomenon and thus the resulting NLs correspond to the same experimental
uncertainty which would allow to detect them. We checked this numerically by a
Monte Carlo simulation and confirmed that the experimental errors at which the
NLs become visible is the same for a) and b) methods of the NL determination.

We now address the question of how accurately the present theory can predict
the NLs of King plots within the standard model. Theoretical accuracy of
isotope-shift (particularly, the specific mass shift) constants of many-
electron systems achievable in modern calculations is not very high—on the
level of few percent. Initially, it was assumed that, because the NLs are very
small, even an approximate theoretical result would provide valuable
information for searching for new physics with Kings plots. However, as has
been shown in our previous work [19], the numerical values of NLs in
hydrogenlike ions turn out to be very sensitive to the experimental
uncertainties of the nuclear radii, which places limits on theoretical
predictions of NLs. In the present work, we studied the sensitivity of the
numerical values of NLs in Ca+ to the theoretical uncertainties of isotope-
shift parameters and to the experimental errors of nuclear radii. Since the
standard methods of error propagation do not work in this case, we used a
Monte Carlo simulation. Each of the computed isotope-shift constants was
represented by a set of quasirandom normally-distributed numbers, with
parameters of the normal distribution defined by the central value and the
uncertainty of the corresponding isotope-shift constant. The NL errors were
obtained by analyzing the resulting distributions of the NL values. The
uncertainty estimates described in Sec. III.5 were employed for the isotope-
shift constants. For the nuclear radii, we used the errors from Ref. [26],
specifically the systematic uncertainty of 0.0020 fm for all isotopes and the
relative uncertainty of 0.0009 fm of nuclear radii of isotopes $A=$42, 44, 46,
and 48 relative to $A=40$.

The numerical results are summarized in Table 8. The two methods of the NL
determination yield very similar results, therefore it is sufficient to
present the values only for the method a). We find that for the
($3d_{3/2},\,3d_{5/2}\rightarrow 4s$) pair the theoretical uncertainty
$\epsilon_{\mathrm{const}}$ is an order of magnitude larger than the central
value. This is not surprising, given the strong cancellations between the
isotope shifts of the two transitions. As is evident from Table 5, we cannot
reliably predict the isotope shift of the $3d_{5/2}$–$3d_{3/2}$ fine-structure
difference, and this fact leads to a large uncertainty in the corresponding
NLs. Hence, for the ($3d_{3/2},\,3d_{5/2}\rightarrow 4s$) pair only an upper
limit of NL can be obtained: we conclude that the standard-model NLs for this
transition pair in Ca+ might be observed when the experimental accuracy is
below 200 Hz. As a matter of fact, such accuracy is already achievable in
modern experiments with Ca+ [13, 14], although no confirmed NLs have been
reported so far.

The second pair of transitions considered in this work, ($4p_{1/2}\rightarrow
4s$; $3d_{3/2}\rightarrow 4s$), involves transitions with different principal
quantum numbers $n$, hence the cancellation between the isotope shifts is much
weaker. Accordingly, the NLs predicted for this KP transition pair are
significantly larger and can be seen already at the experimental accuracy of
about 1.5 kHz. The purely theoretical uncertainty—which stems from the
uncertainty of isotope-shift constants—amounts to $\sim$20% of the NL value,
which means that, in this case, our theory can quantitatively predict a NL.
Additionally, the nuclear-radii uncertainty turns out to be about twice as
large as the purely theoretical one. This means that if observed, the
($4p_{1/2}\rightarrow 4s$; $3d_{3/2}\rightarrow 4s$) King-plot can be used to
improve our knowledge of nuclear radii.

Table 6: Two-dimensional King-plot NLs. Methods a) and b) are depicted in Fig. 1. The NL units are Hz for method a) and kHz2 for method b). The ‘3rd point’ column corresponds to the isotope shifts in $A=42,44,46$ isotopes and the ‘4th point’ to the isotope shifts in $A=42,44,48$ isotopes with respect to $A_{a}=40$. Approximate percentages of each higher-order term contribution to the total NLs are given in respective columns. The abbreviations for the individual terms are the same as in Table 5.  transitions | contribution | 3rd point | 4th point
---|---|---|---
|  | a) | b) | a) | b)
$3d_{3/2}\rightarrow 4s;$ | total | 44 | 1.1$\times 10^{3}$ | 180 | 3.3$\times 10^{3}$
$3d_{5/2}\rightarrow 4s.$ | MS, sec. | 99.2% | 99.2% | 24.7% | 24.7%
| FS, h.o. | 0.1% | 0.1% | 0.0% | 0.0%
| FS, sec. | 0.6% | 0.6% | 0.2% | 0.2%
| np | 0.2% | 0.2% | 75.1% | 75.0%
| cross term | 0.0% | 0.0% | 0.0% | 0.0%
$4p_{1/2}\rightarrow 4s;$ | total | 1.4$\times 10^{4}$ | 2.8$\times 10^{5}$ | 1.4$\times 10^{4}$ | 2.2$\times 10^{5}$
$3d_{3/2}\rightarrow 4s.$ | MS, sec. | 99.7% | 99.7% | 99.5% | 99.5%
| FS, h.o. | 0.1% | 0.1% | 0.1% | 0.1%
| FS, sec. | 0.0% | 0.0% | 0.1% | 0.1%
| np | 0.2% | 0.2% | 0.3% | 0.3%
| cross term | 0.0% | 0.0% | 0.0% | 0.0%
Table 7: The three-dimensional King plot NLs, see Eq. (53). Contributions to the three-dimensional NL are non-additive; we consecutively exclude each term from the total sum and record how the NL value is impacted. The second column names the one term excluded from the sum, whereas the third column shows the change of the NL in percentage of the ‘total’ value with all terms included. The abbreviations for the individual terms are the same as in Table 5. transitions |  contribution | NL[kHz3]
---|---|---
$4p_{1/2}\rightarrow 4s$; |  total | 8.2$\times 10^{3}$
$3d_{3/2}\rightarrow 4s$; |  no MS, sec. | 99.9%
$3d_{5/2}\rightarrow 4s$. |  no FS, h.o. | 0.0%
|  no FS, sec. | 0.0%
|  no np | 99.9%
|  no cross term | 0.0%
Table 8: NLs and their uncertainties for the two-dimensional King plots. For each entry, the upper line shows a NL calculated with method a) and its uncertainties, while the lower line gives the maximal experimental uncertainties of (both) transition frequencies that allow to detect the NL values listed in the upper line. $\epsilon_{\mathrm{const}}$ is the theoretical error of the NL value stemming from the numerical uncertainty of isotope-shift constants, whereas $\epsilon_{\mathrm{radii}}$ is the NL uncertainty originating from the experimental values of nuclear radii. The units are kHz.  Transitions |  | Value | $\epsilon_{\mathrm{const}}$ | $\epsilon_{\mathrm{radii}}$
---|---|---|---|---
$3d_{3/2}\rightarrow 4s;$ | 3rd point | 0.044 | 1.5 | 0.013
$3d_{5/2}\rightarrow 4s.$ |  | (0.005) | (0.20) | (0.002)
| 4th point | 0.181 | 1.5 | 0.023
|  | (0.012) | (0.12) | (0.002)
$4p_{1/2}\rightarrow 4s;$ | 3rd point | 14 | 2.4 | 4.0
$3d_{3/2}\rightarrow 4s.$ |  | (1.3) | (0.3) | (0.45)
| 4th point | 14 | 2.5 | 7.0
|  | (0.8) | (0.17) | (0.4)

##  V Conclusion

We performed a detailed study of isotope shifts in Ca+ transition energies.
The calculations of the first-order mass-shift and field-shift isotope-shift
constants were carried out within the relativistic many-body perturbation
theory for the $4s$, $4p_{1/2,\,3/2}$, and $4d_{3/2,\,5/2}$ states of Ca+. The
results are in good agreement with the previous calculations which were
performed using many-body perturbation theory [22] and coupled-cluster method
[23].

For the first time, the higher-order isotope-shift effects responsible for the
standard-model (SM) nonlinearities in King plots were calculated: second-order
mass shift, second- and higher-order field shifts, and nuclear polarization.
We analyzed the resulting Kings-plot nonlinearities in Ca+ with three
different methods and demonstrated that the dominant contributions originate
from the second-order mass shift and nuclear polarization.

Two pairs of transitions were examined in this work. The first pair,
($3d_{3/2}\rightarrow 4s$; $3d_{5/2}\rightarrow 4s$), consists of two narrow
transitions that can be measured very accurately, which makes it suitable for
new-physics searches with King plots. However, for this pair the present
theory can only predict an upper limit of the nonlinearities caused by the SM
effects, owing to large cancellations between the isotope shits. Based on our
results, we conclude that the SM nonlinearities in this transition pair can
become visible when the experimental accuracy is below 200 Hz, an accuracy
already achieved in modern experiments with Ca+. So, as long as experiments do
not detect any King-plot nonlinearity, this fact can be used for placing
bounds on possible new-physics effects, see Ref. [2]. However, as soon as the
King plot is measured to be nonlinear, it would be hard to discern whether
this nonlinearity is due to a SM effect or a new-physics interaction.

The second transition pair investigated in this work, ($4p_{1/2}\rightarrow
4s$; $3d_{3/2}\rightarrow 4s$), involves transitions with different principal
quantum numbers. The $4p$-$4s$ transition is rather wide and cannot be
measured as accurately as the $3d$-$4s$ one. However, the cancellation between
the isotope shifts is relatively small in this case and the theory can provide
a quantitative prediction of the SM nonlinearity. As follows from Table 8, the
King-plot nonlinearity for this transition pair can be observed already at the
1 kHz level of experimental uncertainty, which is feasible in modern
experiments. We conclude that the ($4p_{1/2}\rightarrow 4s$;
$3d_{3/2}\rightarrow 4s$) transition pair is a promising candidate for the
experimental identification of the SM King-plot nonlinearity and for obtaining
information about nuclear charge radii.

## Acknowledgements

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC-2123
QuantumFrontiers – 390837967 and under the project No. SU 658/4-2.

## References

  * Safronova et al. [2018] M. Safronova, D. Budker, D. DeMille, D. F. J. Kimball, A. Derevianko, and C. W. Clark, Reviews of Modern Physics 90, 025008 (2018).
  * Berengut et al. [2018] J. C. Berengut, D. Budker, C. Delaunay, V. V. Flambaum, C. Frugiuele, E. Fuchs, C. Grojean, R. Harnik, R. Ozeri, G. Perez, et al., Phys. Rev. Lett. 120, 091801 (2018).
  * Counts et al. [2020] I. Counts, J. Hur, D. P. L. Aude Craik, H. Jeon, C. Leung, J. C. Berengut, A. Geddes, A. Kawasaki, W. Jhe, and V. Vuletić, Phys. Rev. Lett. 125, 123002 (2020).
  * Schmidt et al. [2005] P. O. Schmidt, T. Rosenband, C. Langer, W. M. Itano, J. C. Bergquist, and D. J. Wineland, Science 309, 749 (2005).
  * Wineland et al. [2002] D. J. Wineland, J. C. Bergquist, J. J. Bollinger, R. E. Drullinger, and W. M. Itano, in _Frequency Standards and Metrology_ (World Scientific, 2002), pp. 361–368.
  * King et al. [2022] S. A. King, L. J. Spieß, P. Micke, A. Wilzewski, T. Leopold, E. Benkler, R. Lange, N. Huntemann, A. Surzhykov, and V. A. Yerokhin, Nature 611, 43 (2022).
  * Delaunay  _et al._ [2017] author C. Delaunay, R. Ozeri, G. Perez, and Y. Soreq, journal Phys. Rev. D 96, 093001 (2017).
  * Berengut et al. [2020] J. C. Berengut, C. Delaunay, A. Geddes, and Y. Soreq, Phys. Rev. Research 2, 043444 (2020).
  * Drake et al. [2021] G. W. F. Drake, H. S. Dhindsa, and V. J. Marton, Phys. Rev. A 104, L060801 (2021).
  * Hur et al. [2022] J. Hur, D. P. Aude Craik, I. Counts, E. Knyazev, L. Caldwell, C. Leung, S. Pandey, J. C. Berengut, A. Geddes, W. Nazarewicz, et al., Phys. Rev. Lett. 128, 163201 (2022).
  * Figueroa et al. [2022] N. Figueroa, J. Berengut, V. Dzuba, V. Flambaum, D. Budker, and D. Antypas, Phys. Rev. Lett. 128, 073001 (2022).
  * Allehabi et al. [2021] S. O. Allehabi, V. A. Dzuba, V. V. Flambaum, and A. V. Afanasjev, Phys. Rev. A 103, L030801 (2021).
  * Knollmann et al. [2019] F. W. Knollmann, A. N. Patel, and S. C. Doret, Phys. Rev. A 100, 022514 (2019).
  * Solaro et al. [2020] C. Solaro, S. Meyer, K. Fisher, J. C. Berengut, E. Fuchs, and M. Drewsen, Phys. Rev. Lett. 125, 123003 (2020).
  * Gebert et al. [2015] F. Gebert, Y. Wan, F. Wolf, C. N. Angstmann, J. C. Berengut, and P. O. Schmidt, Phys. Rev. Lett. 115, 053003 (2015).
  * Plunien et al. [1991] G. Plunien, B. Müller, W. Greiner, and G. Soff, Phys. Rev. A 43, 5853 (1991).
  * Nefiodov et al. [1996] A. V. Nefiodov, L. N. Labzowsky, G. Plunien, and G. Soff, Phys. Lett. A 222, 227 (1996).
  * Yerokhin et al. [2019] V. A. Yerokhin, K. Pachucki, and V. Patkóš, Annalen der Physik 531, 1800324 (2019).
  * Yerokhin et al. [2020] V. A. Yerokhin, R. A. Müller, A. Surzhykov, P. Micke, and P. O. Schmidt, Phys. Rev. A 101, 012502 (2020).
  * Blundell et al. [1987] S. A. Blundell, D. S. Guo, W. R. Johnson, and J. Sapirstein, Atomic Data and Nuclear Data Tables 37, 103 (1987).
  * Johnson et al. [1988] W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 2764 (1988).
  * Safronova and Johnson [2001] M. S. Safronova and W. R. Johnson, Phys. Rev. A 64, 052501 (2001).
  * Dorne et al. [2021] A. Dorne, B. K. Sahoo, and A. Kastberg, Atoms 9, 26 (2021).
  * Berengut et al. [2003] J. C. Berengut, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A 68, 022502 (2003).
  * Raman et al. [2001] S. Raman, C. W. Nestor, and P. Tikkanen, Atomic Data and Nuclear Data Tables 78, 1 (2001).
  * Angeli and Marinova [2013] I. Angeli and K. P. Marinova, Atomic Data and Nuclear Data Tables 99, 69 (2013).
  * Müller et al. [2021] R. A. Müller, V. A. Yerokhin, A. N. Artemyev, and A. Surzhykov, Physical Review A 104, L020802 (2021).
  * Shabaev et al. [2004] V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Phys. Rev. Lett. 93, 130405 (2004).
  * Flambaum et al. [2021] V. V. Flambaum, I. B. Samsonov, H. B. T. Tan, and A. V. Viatkina, Phys. Rev. A 103, 032811 (2021).
  * Puchalski et al. [2006] M. Puchalski, A. M. Moro, and K. Pachucki, Phys. Rev. Lett. 97, 133001 (2006).
  * King [2013] W. H. King, _Isotope Shifts in Atomic Spectra_ (Springer Science & Business Media, 2013).