text:
\section{Introduction}
\label{section1}

The Mg$^+$ 3s-3p fine structure doublet near 280nm (see fig.~\ref{fig:Mg_levels}) is a prominent feature in many astronomical spectra. It can be observed by ground based telescopes, once the cosmological redshift translates the strong UV lines into the transparency range of the atmosphere. This makes the transitions good probes of physics at early epochs. One example for their astrophysical importance is the ongoing quest for a cosmological space-time variation of the fine structure constant~$\alpha$. An investigation of narrow absorption lines from gas clouds which are back-illuminated by distant quasars (quasi-stellar objects, QSO), allows to compare the value of~$\alpha$ in distant regions of the universe with the value measured on earth~\cite{Bahcall1967}. The detection of a variation of the electromagnetic coupling constant~$\alpha$ would provide compelling evidence of physics beyond the standard model.

The Mg$^+$ doublet was first used in 1976~\cite{Wolfe} to constrain variations of~$\alpha$ via the alkali doublet method~\cite{Savedoff1956}. An analysis of fine structure doublets allows to distinguish between a variation in $\alpha$ and a frequency change due to the cosmological redshift, since the ratio between the fine structure and gross structure of atomic energy levels scales as~$\alpha^2$ in leading order. Dzuba et al.~introduced the many multiplet (MM) method in 1999~\cite{Dzuba99}, a refined approach which utilizes the fact that the transition frequencies itself contain relativistic corrections $\propto\alpha^2$. The relativistic corrections are most pronounced for the ground states, so the comparison of transitions in different atomic species from the same object can significantly improve the sensitivity.

\begin{figure}[b!]
\includegraphics[width=0.7\columnwidth]{batteiger_fig1.eps}
\caption{The Mg$^+$ resonance doublet near 280 nm is isoelectronic to the 3s-3p doublet in Na historically labeled as D$_1$ and D$_2$ Fraunhofer lines. From the D$_1$ and D$_2$ transition frequencies measured in $^{24}$Mg$^+$ and $^{26}$Mg$^+$ we deduce isotope shifts and fine structure splittings (FS). The measurement accuracy allowed to resolve the isotope shift of the fine structure splitting for the first time.}
\label{fig:Mg_levels}
\end{figure}

Webb and coworkers applied the MM method to Fe$^+$ transitions measured against the Mg$^+$ doublet and found indications that~$\alpha$ was smaller than the present day value for redshifts~$z\cong1$ (about $8\times10^9$~years ago,~$z=(\nu_{emit}-\nu_{obs})/\nu_{obs}$)~\cite{Webb99}. The same group extended their studies to larger samples and further transition lines~\cite{Murphy2001,Murphy2003}, which finally resulted in a statistically significant signature of a smaller fine structure constant $\Delta\alpha/\alpha=(-0.57\pm0.11)\times10^{-5}$. Hereby, 143 QSO absorption samples where evaluated (77 of which were low-$z$ Mg/Fe systems~\cite{Murphy2004}), all observed by the Keck observatory using the HIRES spectrometer~\cite{keckvlt}. Attempts to reproduce the results with 23 Mg-based samples observed by a different telescope and spectrometer (VLT/UVES~\cite{keckvlt}) suggest a drift consistent with zero~$\Delta\alpha/\alpha=(-0.06\pm0.06)\times10^{-5}$~\cite{Srianand2004}. This work and further investigations consistent with zero drift are summarized and disputed in~\cite{Murphy2008}. Recent ''Single Ion Differential $\alpha$ Measurement'' (SIDAM) studies on Fe$^+$ lines also gave a result compatible with zero drift~\cite{Molaro2008}. Next generation telescopes~\cite{Molaro2007}, possibly in combination with absolute frequency calibration~\cite{Murphy2007,Steinmetz2008,Braje2008,Li2008} may help to clarify the discrepancies on the reported uncertainty levels.

Comparison with QSO absorption spectra requires laboratory reference data with an accuracy of better than 10$^{-4}$\AA~\cite{Berengut2006}, which corresponds to about 40~MHz at 280~nm. Absolute frequency calibrated measurements~\cite{comb} in the favorable environment of an ion trap, as reported in this article, provides reference data~\cite{Wolf2008,Wolf2009,Herrmann2009} that may be regarded as exact for comparison with astrophysical spectra. Besides high accuracy it is particularly important to provide isotopically resolved data. Magnesium has three stable isotopes with mass numbers 24, 25 and 26 (natural abundance 79:10:11). Different isotope abundances in the investigated astronomical objects were identified as an important source of systematic errors in the MM analysis, moreover, an investigation of the evolution of isotope abundances itself is of interest~\cite{Murphy2001,Kozlov2004,Ashenfelter2004,Fenner2005}. For these reasons an isotope shift measurement on the D$_1$ transition was requested in~\cite{Berengut2006}. Also comparison with isotope shift theory demands refined reference data since the accuracy of theoretical predictions exceeded the measurement uncertainty \cite{berengut2003,Drullinger1980}.

So far the only isotopically resolved measurement was provided by Drullinger et al.~\cite{Drullinger1980}, where the $^{26,25,24}$Mg$^+$ isotope shifts of the D$_{2}$ line were determined with an uncertainty of 100~MHz. The most accurate literature values for the D$_{1}$ and D$_{2}$ transition frequencies originate from Fourier transform (FT) spectrometry using a hollow-cathode discharge lamp here an uncertainty of 60~MHz was reported for unresolved lines of natural isotopic abundance~\cite{Pickering1998,Pickering2000,Aldenius2006}. MM~analysis used laboratory reference data for the individual isotopic (and hyperfine structure) components~\cite{Murphy2003} that are derived from the unresolved FT measurements~\cite{Pickering1998,Pickering2000,Aldenius2006} in combination with the D$_{2}$ isotope shift results~\cite{Drullinger1980}. Our measurement of the D$_{1}$ and D$_{2}$ component of the $^{24}$Mg$^+$ and $^{26}$Mg$^+$ isotopes, improves the transition frequencies by more than two orders of magnitude. From these four absolute frequency measurements we derive $^{26-24}$Mg$^+$ isotope shifts. The isotope shift of the D$_{1}$~transition has not been measured before and the D$_{2}$ value is improved by three orders of magnitude. The accuracy of our measurement allows to resolve the $^{26-24}$Mg$^+$ isotope shift of the fine structure splitting. From the $^{24}$Mg$^+$ and $^{26}$Mg$^+$ measurements we additionally deduce $^{25}$Mg$^+$ transition frequencies, including the hyperfine structure, with an accuracy that is better than requested in~\cite{Berengut2006}.

\begin{figure}[tb]
\includegraphics[width=1\columnwidth]{batteiger_fig2.eps}
\caption{Basic idea. A tightly focused laser beam cools ions at one side only. A weaker spectroscopy beam is applied on axis to be insensitive to possible radial micro motion. An imaging photo detector collects the essentially background free signal from single, sympathetically cooled ions.}
\label{Chain}
\end{figure}

Finally, we note that strong constraints on \emph{present-day} variations of~$\alpha$ are set terrestrially by a comparison of optical clocks based on different atomic transitions. Though optical clocks may only be compared on time scales of years, the tremendous accuracy allows to constrain variations of~$\alpha$ about an order of magnitude more accurately than MM studies, if a linear drift over 10$^{10}$~years look back time is assumed. Recent measurements of an Al$^+$ clock against a Hg$^+$ clock at NIST suggested ${\dot{\alpha}}/{\alpha}=(-1.6\pm2.3)\times10^{-17}$/year~\cite{Rosenband2008}.

\section{Method}
\label{section2}

\begin{figure}[tb]
\includegraphics[width=1\columnwidth]{batteiger_fig3.eps}
\caption{Wide scan across the $^{24}$Mg$^+$ D$_1$ resonance. Fluorescence photons were collected for 3s per data point, the scan was randomized in order. For our measurements we restricted the tuning range to 180 MHz around the line center. Inset: For comparison we show an asymmetric $^{24}$Mg$^+$ D$_2$ line profile scanned from the low frequency side, overlayed is a Voigt fit to the rising edge, which is centered around the measured transition frequency. Here the ion was essentially uncooled, only a weak cooling beam was superimposed to prevent ion loss during the scan time.}
\label{fig:widescan2}
\end{figure}

 \begin{figure*}[tb]
\includegraphics[width=1.6\columnwidth]{batteiger_fig4.eps}
\caption{Schematic experimental setup. All involved lasers are phase locked (PLL) to a GPS disciplined hydrogen maser. The phase coherent link between the RF reference and the UV regime via a fiber laser frequency comb and a transfer diode laser is shown for the spectroscopy dye laser. The spectroscopy beam is intensity stabilized with an AOM and spatially filtered with a 15 $\mu$m pinhole. A quarterwave plate (QWP) rotates the polarization state of the spectroscopy beam. D$_2$ cooling light is provided by a frequency quadrupled fiber laser~\cite{Friedenauer2006}. An alternative setup generating both beams from a single laser output is shown in~\cite{Herrmann2009}.}
\label{fig:Setup}
\end{figure*}

Lines observed in QSO absorption spectra are usually dipole-allowed transitions from the ground state~\cite{Berengut2006}. Spectroscopy on such transitions in trapped ions requires different spectroscopic approaches than probing narrow clock transitions. Spectroscopy in the favorable sideband resolved regime~\cite{Bergquist1987}, that allows to avoid first order Doppler and recoil shifts, is hard to establish, since typical transition bandwidths $\Gamma_{trans}$ are broad compared to practical secular trap frequencies $\omega_{sec}$. In the weak binding limit ($\omega_{sec}\ll\Gamma_{trans}$) the absorption spectrum is strongly asymmetric if no further measures are taken and the transition is probed with a single laser beam. In this regime the motional sidebands overlap and for constant temperature a Voigt profile would be observed, just as for an unbound atom~\cite{Wineland1979}. However, scattering photons changes the temperature of the ion; on the low-frequency (``red'') side of the resonance the ion is cooled to a detuning dependent temperature, on the high-frequency (``blue'') side it is heated and expelled from the trap center. Therefore the Voigt spectrum is heavily distorted: On the red side due to a detuning dependent equilibrium temperature and thus varying Gaussian contribution to the linewidth. Near resonance the temperature diverges and the ion is lost, so the blue side is not observed~\cite{Nagourney1983} (see also fig.~\ref{fig:widescan2} inset). Therefore, in view of precision spectroscopy, the crux is to suppress adverse temperature variations induced by the spectroscopy laser. Since the statistical uncertainty determining the line center is given by the linewidth divided by the signal-to-noise ratio of the data, a suitable approach should ensure low temperature (i.e.~narrow linewidth) and a high scattering rate. Previously used approaches include a two laser method~\cite{Drullinger1980}, a double resonance scheme~\cite{Wineland1981}, buffer gas cooling~\cite{Nakamura2006} and chopped detection~\cite{Wolf2008}. Our novel spectroscopy technique allows to record essentially unperturbed Voigt profiles for the first time. The method is discussed in detail in~\cite{Herrmann2009}, here we will only sketch the basic idea. We prepare a crystallized chain of ions which is continuously Doppler cooled at one side only (see fig.~\ref{Chain}). A weaker spectroscopy laser probes sympathetically cooled ions at the other end of the chain, where we collect a basically background free fluorescence signal from single ions using an imaging photon detection
system. This avoids both the background and ac Stark shift from the cooling laser. Counterintuitive to a thermal conductivity picture, the symmetric eigenmode structure efficiently couples cooling and spectroscopy ions at opposite ends of the chain, i.e.~the sympathetically cooled ions are cooled just as efficiently as directly cooled ions. A detailed analysis of the considered transitions shows: Residual temperature variations lead to systematic shifts of the line center always below $10^{-4}$ of the linewidth, as long as the cooling laser damping exceeds the heating induced by the spectroscopy beam for all detunings. Therefore, the method enables to observe virtually unperturbed resonances from single, cold ions with high signal-to-noise ratio (see fig.~\ref{fig:widescan2}).

\section{Setup, Measurement and Results}
\label{section3}

\begin{figure*}[htb]
  \centering
   \subfigure[~$^{26}$Mg$^+$~D$_2$]{\includegraphics[width=0.65\columnwidth]{batteiger_fig5a.eps}}\qquad
   \subfigure[~$^{24}$Mg$^+$~D$_1$]{\includegraphics[width=0.62\columnwidth]{batteiger_fig5b.eps}}\qquad
   \subfigure[~$^{26}$Mg$^+$~D$_1$]{\includegraphics[width=0.62\columnwidth]{batteiger_fig5c.eps}}
  \caption{Measurement of two isotopic components of the D$_1$ and D$_2$ transition. Outer dashed lines represent the error bar on the absolute transition frequencies due to correlated amplitude and phase modulation, inner dotted lines the reduced error bar for the relative determination of isotope shifts and fine structure splittings. Squares indicate measurements with two laser systems according to fig.\ref{fig:Setup}, circles measurements with the double pass AOM scheme described together with the measurement of the $^{24}$Mg$^+$~D$_2$ component in ~\cite{Herrmann2009}.}
  \label{fig:Data}
\end{figure*}

The experimental setup is described in~\cite{Herrmann2009}, here we summarize basic features and illuminate further aspects in greater detail. The Mg$^+$ ions were stored in a linear Paul trap operated at 15.8~MHz. The radial secular frequencies were about 1~MHz, while the axial secular frequency was chosen to be about 60~kHz. Ions were loaded into the trap by two-step photo ionization via the 3s$_0$-3p$_1$ transition in neutral Mg with a frequency doubled dye laser at 285~nm~(Coherent 699). References~\cite{Kjaergaard2000,Gulde2001,Lucas2004,Steele2007,Tanaka2007} showed that isotope selective photo ionization is possible when atom and laser beams are carefully arranged. Our trap was loaded from neutral Mg vapor of natural abundance (Doppler broadening: about 3~GHz, $^{26-24}$Mg~3s$_0$-3p$_1$ isotope shift: 1.4~GHz~\cite{Salumbides2006}). Yet we observed purified chains of up to 30 $^{26}$Mg$^+$ ions, when we applied additional radio frequencies resonant to the $^{24}$Mg$^+$ and $^{25}$Mg$^+$ radial secular frequencies to an electrode parallel to the trap axis during the loading period. Cooling and spectroscopy light was either provided by another frequency doubled dye laser (Coherent 699) or a frequency quadrupled Yb fiber laser system~\cite{Friedenauer2006} at 1118~nm (Koheras Boostik), which could only interact with the D$_2$ transition (90~GHz tuning range). A schematic setup using both lasers is shown in fig.~\ref{fig:Setup}. After a failure of the fiber laser the cooling and spectroscopy beams were produced by the frequency doubled dye laser. In this case the UV output was split and passed through two double pass acousto-optic modulators  (AOM) setups~\cite{DPAOM}, which allowed to independently control the frequency and amplitude of both beams. The intensities were set to $s_{cool}\sim0.5$ and $s_{probe}\sim7\times10^{-4}$ in units of the saturation intensity ($I_{sat}\sim$~2.47~mW/mm$^2$). The cooling beam was detuned 64~MHz below resonance, the practical tuning range of the spectroscopy laser was 180~MHz, limited by the efficiency of the scanning AOM in double pass configuration. The cooling beam was tightly focused to a beam radius of $w_0\sim20\,\mu$m and enclosed an angle of about 15$^\circ$ with the trap axis predefined by the optical access of our vacuum vessel. Given that geometry we found optimal experimental conditions when we cooled 2-3 outer ions of a chain containing 8-10 ions. This allowed us to perform low background spectroscopy on 2-3 ions at the other end of the chain. Working with shorter chains decreased the number of low background spectroscopy ions, longer chains increased the probability that a MgH$^+$ ion was generated by a chemical reaction with the background gas. Dark MgH$^+$ ions embedded in the chain changed their site on time scales of few seconds. Once they enter the cooling or spectroscopy region it turns dark which distorts the line profile. We therefore reloaded our trap once an impurity ion was observed. The spectroscopy beam collinear with the ion chain was stabilized in intensity, spatially filtered with a small pinhole and loosely focused to about $w_0\sim200\,\mu$m in the trapping region. Additional spatial filtering was necessary since tuning of the double pass AOM led to frequency dependent distortions of the beam profile, observable with a beam profiling camera. We stabilized the \emph{total} intensity, but the ions sampled frequency dependent \emph{local} intensity variations, which resulted in an apparent alignment sensitivity of the line center. Filtering with a~15~$\mu$m size pinhole removed the beam profile distortions and the alignment sensitivity of the line center.

The RF reference for all involved frequency combs, phase locked loops and the spectroscopy AOM was provided by a hydrogen maser. A comparison between the free-running maser and the global positioning system showed a deviation smaller than $1.7\times10^{-13}$ over the whole measurement period, which is sufficient for our purpose. The phase coherent link between the RF reference and the optical regime was achieved with a 100~MHz repetition rate erbium doped fiber laser frequency comb. The frequency comb reaches well into the fourth sub-harmonic of the Mg$^+$ doublet at about 1120~nm. To bridge the gap to the green dye laser output we phase locked a diode laser at 1120~nm to a mode of the frequency comb and recorded a beat note between its frequency doubled output and the dye laser (see fig.~\ref{fig:Setup}). A synthesizer mixed the beat frequency down to a constant offset lock frequency, which allowed us to arbitrarily set the frequency of the phase locked dye laser with respect to the IR comb modes. This technique was also used to tune the spectroscopy laser, when we cooled the ion chain with the independent fiber laser system. The mode number of the IR comb, used for the measurement of the $^{24}$Mg$^+$~D$_2$ line was determined with a second frequency comb with repetition rate of 250 MHz (Menlo Systems)~\cite{Herrmann2009}.
We deduced the mode number for the remaining transitions relative to the known $^{24}$Mg$^+$~D$_2$ component using two independent wavemeters with an relative accuracy considerably better than the repetition rate (Burleigh WA 1500, High Finesse WS7).

Every single measurement consisted of 30 - 40 frequency settings within a 180~MHz ($>$~4~$\Gamma$) wide range around the line center, fluorescence photon counts were recorded for 3~s each. The chronological order of probe laser frequencies was randomized for each line scan to be insensitive to a great variety of possible drifts. Photons were collected with a f/2 imaging system attached to a single photon camera (Quantar Mepsicron II). The spatially resolved count rate was digitized in 512~$\times$~512 pixel images. To read out the fluorescence of a single ion we set a circular region of interest around the ion using software (see fig.~\ref{Chain}).

\begin{table*}
\caption{Absolute transition frequencies for individual isotopic components. Pickering et al.~do not explicitly specify an uncertainty for the $^{24}$Mg$^+$ component (see table~1 in \cite{Pickering1998}), for comparison we adopt the composite line uncertainty of 60~MHz.}
\begin{tabular}{|l|c|c|c|}
\hline
Measured transitions~~~~&  This work &~~~Drullinger et al.~\cite{Drullinger1980}~~~&~~~Pickering et al.~\cite{Pickering1998}~~~\\
\hline
\hline
~~$^{24}$Mg$^+$~D$_1$ (3s$_{1/2}$-3p$_{1/2}$) &~~~~1 069 338 342.56 (16) MHz~~~~& - & 1 069 338 293 (60) MHz \\
\hline
~~$^{26}$Mg$^+$~D$_1$ (3s$_{1/2}$-3p$_{1/2}$) &~~~~1 069 341 427.47 (16) MHz~~~~& - & - \\
\hline
~~$^{24}$Mg$^+$~D$_2$ (3s$_{1/2}$-3p$_{3/2}$)~\cite{Herrmann2009}~~~&~~~~1 072 082 934.33 (16) MHz~~~~&~~~1 072 082 833 (120) MHz~~~& 1 072 082 862 (60) MHz   \\
\hline
~~$^{26}$Mg$^+$~D$_2$ (3s$_{1/2}$-3p$_{3/2}$) &~~~~1 072 086 021.89 (16) MHz~~~~&~~~1 072 085 883 (156) MHz~~~& -  \\
\hline
\end{tabular}
\label{TableResults}
\end{table*}

The D$_2$ transition shows two magnetic sublevels for the s$_{1/2}$ ground state and four magnetic sublevels for the p$_{3/2}$ excited state. Magnetic fields in direction of the spectroscopy laser therefore shift the line center if the polarization is not perfectly linear. To compensate such shifts we introduced a quarter wave plate (QWP) into the beam path and measured the line center as function of the polarization, by rotating the QWP in 45$^\circ$ steps. We observed the expected sinusoidal modulation of the line center with an amplitude smaller than 950~kHz, a sinusoidal fit to the data compensates linear Zeeman shifts. Each measurement day in fig.~\ref{fig:Data}a consisted of at least one QWP rotation scan. Binning the data into measurement days serves as an important consistency check against possible systematics induced by operators and the lab environment. Since the s$_{1/2}$ and p$_{1/2}$ levels show two magnetic sublevels each, only line broadening is expected for the D$_1$ line. Nevertheless we continued with the QWP evaluation scheme to rule out unexpected polarization effects on the line center. No significant polarization dependence of the D$_1$ line was observed.

We assume that the limiting systematic uncertainty is introduced by possible asymmetries in the rather broad spectrum of the frequency doubled dye laser, which was about 10~MHz. The phase locked loop to the frequency doubled diode laser stabilizes the carrier phases~(see fig.~\ref{fig:Setup}). However, correlated amplitude- and phase modulation (AM/PM) can lead to asymmetries in the frequency spectrum with respect to the carrier frequency. Therefore the center-of-gravity of the spectroscopy laser spectrum can shift against its absolute frequency. Since the fluorescence profile is a convolution of the atomic response and the laser spectrum, the counted frequency may differ from the true transition frequency. To quantify the impact of this effect, we measured the in-loop heterodyne beat note between the diode laser and the dye laser with a spectrum analyzer referenced to the same RF reference as the involved offset lock frequency. We fitted a Gaussian profile to the beat note and found an average deviation from the offset lock frequency of 80(60)~kHz at 560~nm, therefore we assume an uncertainty of 160 kHz in the UV, an order of magnitude larger than e.g.~the 14~kHz statistical uncertainty of the $^{26}$Mg$^+$~D$_2$ measurement. A consistency check with the frequency quadrupled fiber laser directly locked to the IR frequency comb gave good agreement within the statistical error bars~\cite{Herrmann2009}. All other considered systematic effects, including method inherent shifts due to residual detuning dependent temperature variations, are much smaller: ac Stark shift due to residual background from the cooling laser (30 kHz), dc Stark shift from the trapping fields (0.1~Hz), line shape model (270~Hz), RF reference (180~Hz) and 2nd order Doppler shift (-0.3~Hz). We corrected the recoil shift of -106~kHz for $^{24}$Mg$^+$ and -98~kHz for $^{26}$Mg$^+$ transition frequencies.

We measured the $^{26}$Mg$^+$~D$_2$ line on six measurement days (152 line scans, statistical uncertainty~14~kHz) and the $^{24}$Mg$^+$~D$_1$ line on five measurement days (119 line scans, statistical uncertainty~20~kHz). Due to technical problems we were only able to measure a total of 67 lines in five
measurement days of the $^{26}$Mg$^+$~D$_1$ component. Still, the statistical uncertainty of the measurement reached 32~kHz and all binned measurement days agree within their error bars~(see fig.~\ref{fig:Data}c). The numerical values of the transition frequencies are given in table~\ref{TableResults}.

The reproducibility of our measurement over time periods of days and months, also indicates that the systematic uncertainty due to correlated AM/PM enters rather as a constant offset on the absolute frequency, than as a time-varying scatter. We therefore assume that the uncertainty of measured isotope shift frequencies and fine structure splittings is smaller than 160 kHz due to cancellation of correlated AM/PM effects. Possible shifts due to residual detuning dependent distortions of the intensity profile induced by the scanning double pass AOM should not influence measurements where the phase locked spectroscopy laser was tuned by changing the offset lock frequency with a synthesizer. Thus, the reproducibility over time and against changes in the setup indicates that the statistical uncertainty of a typical measurement day is an appropriate upper bound on residual systematic deviations due to drifting correlated AM/PM and residual distortions of the intensity profile~(see fig.~\ref{fig:Data}). We therefore quadratically add a systematic deviation of 60 kHz to the statistical uncertainty of the involved transition frequencies, when we determine the relative quantities summarized in table~\ref{TableResultsII}.

\begin{table} [tb]
\caption{Isotope shifts~($\delta\nu^{A',A}$), fine structure splittings (FS) and isotope shift of the fine structure splitting deduced from the D$_1$ and D$_2$ absolute frequency measurements. The error bars are reduced due to common mode rejection of systematic effects.}
\begin{tabular}{|l|c|}
\hline
~$\delta\nu^{26,24}$~D$_1$~(3s$_{1/2}$-3p$_{1/2}$)~~~& 3084.905 (93) MHz  \\
\hline
~$\delta\nu^{26,24}$~D$_2$~(3s$_{1/2}$-3p$_{3/2}$)~& 3087.560 (87) MHz \\
\hline
~$^{24}$Mg$^+$~FS~(3p$_{1/2}$-3p$_{3/2}$)~&~~~2 744 591.767 (88) MHz~~~\\
\hline
~$^{26}$Mg$^+$~FS~(3p$_{1/2}$-3p$_{3/2}$)~&~~~2 744 594.422 (92) MHz~~~\\
\hline
~$\delta\nu^{26,24}$~FS~(3p$_{1/2}$-3p$_{3/2}$)~& 2.66 (13) MHz \\
\hline
\end{tabular}
\label{TableResultsII}
\end{table}

\section{$^{25}$Mg$^+$ isotope shift and hyperfine structure}
\label{section4}

$^{25}$Mg$^+$ has a nuclear spin {\it I}~=~5/2. Its resulting hyperfine structure makes additional beams necessary to establish a cycling transition and repump from dark hyperfine levels, which were not available at the time of the measurement. Yet, having measured the $^{26-24}$Mg$^+$ isotope shifts, we may deduce $^{25}$Mg$^+$ transition frequencies with an uncertainty of about 20~MHz. The frequency shift for a transition between isotopes with atomic masses {\it A'} and~{\it A} respectively can be written as~\cite{berengut2003,Drake}
\begin{equation}
\label{shiftformula}
\delta\nu^{A',A}=(k_{\rm{NMS}}+k_{\rm{SMS}})(\frac{1}{A'}-\frac{1}{A})+F_{FS}~\delta\left\langle r^2\right\rangle^{A',A},
\end{equation}
where $\delta\nu^{A',A}~:=\nu^{A'}-\nu^{A}$, $F_{FS}$ is the field shift constant and $\delta\left\langle r^2\right\rangle^{A',A}$ is the difference between the root-mean-square (rms) nuclear charge radii of the involved isotopes. The mass shift arising from the nuclear recoil is usually treated as the sum of the normal mass shift NMS, that can easily be determined and the specific mass shift SMS, which is difficult to evaluate accurately and has been subject to a variety of theoretical investigations \cite{berengut2003,safronova2001,tupitsyn2003,Korol2007}. The field shift arising from different nuclear charge distributions for different isotopes, turns out to be a small contribution (about 20~MHz) compared to the total mass shift of the $^{26-24}$Mg$^{+}$ D$_1$ and D$_2$ line (about 3~GHz). We express the sum of the mass shift constants in eqn.~(\ref{shiftformula}) in terms of the measured $^{26-24}$Mg$^+$ isotope shift and arrive at an expression
\begin{equation}
\label{shiftformulaII}
\delta\nu^{25,24}=a~\delta\nu^{26,24}+F_{FS}~(\delta\left\langle r^2\right\rangle^{25,24}-a~\delta\left\langle r^2\right\rangle^{26,24}),
\end{equation}
where $a=\frac{m_{26}(m_{25}-m_{24})}{m_{25}(m_{26}-m_{24})}$, that allows to deduce the $^{25-24}$Mg$^+$ isotope shifts using field shift constant determinations $F_{FS}$  from~\cite{berengut2003,safronova2001,tupitsyn2003} and tabulated rms nuclear charge radii, which were measured by electron scattering and muonic atoms x-rays~\cite{angeli2004}. The uncertainty introduced by the measurement uncertainty of the rms nuclear charge radii dominates over the influence of different theoretical field shift constants derived in~\cite{berengut2003,safronova2001,tupitsyn2003,Korol2007}. The results are: $\delta\nu^{25,24}$Mg$^+$~D$_1$~=~1620(19)~MHz and $\delta\nu^{25,24}$Mg$^+$~D$_2$~=~1621(19)~MHz in agreement with~\cite{Drullinger1980}.

Knowing the $^{25-24}$Mg$^+$ isotope shift (i.e.~the hyperfine centroid) one can evaluate transition frequencies between $^{25}$Mg$^+$ hyperfine structure levels. The ground state hyperfine constant was accurately measured in~\cite{Itano1981}. Magnetic dipole and electric quadrupole hyperfine constants for 3p states in $^{25}$Mg$^+$ are evaluated in~\cite{Sur2005,Safronova1998}, uncertainty estimates are not given, but the computation reproduces the experimental value of the ground state hyperfine constant within 0.6\%. The largest shift to a 3p state amounts to 180 MHz (3p$_{1/2}$,{\it F}=2), the uncertainties of the 3p hyperfine constants (assuming that they are smaller than 10\%) may therefore be neglected against the 19~MHz uncertainty in the determination of the $^{25}$Mg$^+$ hyperfine centroid, which is smaller than the transition bandwidth $\Gamma_{trans}$~=~41~MHz and smaller than the 60~MHz requested for the modeling of isotopically unresolved Mg$^+$ lines in astronomical spectra~\cite{Berengut2006}. Using hyperfine constants tabulated in~\cite{Sur2005} one can deduce the whole hyperfine structure of 3p states in $^{25}$Mg$^+$. The absolute frequency of the 3s$_{1/2}$({\it F}=3)-3p$_{3/2}$({\it F}=4) cycling transition, which is widely used in trapped ion experiments, reads for instance: 1~072~085~225~(19)~MHz. This value may be updated when improved hyperfine constants or rms nuclear charge radii become available.

\section{Comparison with literature}
\label{section5}
\subsection{Comparison with Previous Experimental Data \label{section2a}}

\begin{table} [tb]
\caption{Comparison of transition frequencies for a composition of natural isotopic abundance. All frequencies are given in MHz.}
\begin{tabular}{|l|l|l|}
\hline
Reference &~~~~~D$_1$~[MHz] &~~~~~D$_2$~[MHz]~\\
\hline
\hline
This work&~1 069 338 844.2 (1.9)~&~1 072 083 436.4 (1.9)~\\
\hline
Drullinger~\cite{Drullinger1980}&~~~~~~~~~~~~-&~1 072 083 327 (121)\\
\hline
Pickering~\cite{Pickering1998}&~1 069 338 652 (60)&~1 072 083 252 (60)\\
\hline
Aldenius~\cite{Aldenius2006}&~1 069 338 802 (60)&~1 072 083 342 (60)\\
\hline
\end{tabular}
\label{TableComposite}
\end{table}

The only isotopically resolved reference data stood alone for nearly three decades, when Drullinger et al.~\cite{Drullinger1980} measured the D$_{2}$ line in a Penning trap on a cloud of Mg$^+$ ions. The transition frequency of the $^{24}$Mg$^+$ component was measured against an iodine reference line with an uncertainty of 120~MHz and the $^{26-24}$Mg$^+$ isotope shift was measured to be 3050~(100)~MHz in agreement with our results (see table~\ref{TableResults},\ref{TableShift}). So far, there was no isotopically resolved reference data for the D$_1$ line.

The most accurate literature values for the D$_{1}$ and D$_{2}$ transition frequencies originate from isotopically unresolved Fourier transform (FT) spectroscopy, here an uncertainty of 60~MHz was reported for composite lines of natural isotope abundance ~\cite{Pickering1998,Pickering2000,Aldenius2006}. In order to compare isotopically resolved results to unresolved measurements we determine the center-of-gravity for a composition of natural isotope abundance ($^{24}$Mg:~78.99\%, $^{25}$Mg:~10\%, $^{26}$Mg:~11.01\%). According to Land\'{e}'s interval rule $\delta\nu^{25,24}$ represents the isotope shift between the center-of-gravity of the $^{25}$Mg$^+$ hyperfine structure multiplet and the $^{24}$Mg$^+$ component~\cite{Sobelman}, thus the center-of-gravity of the composite line is shifted by $0.1\times\delta\nu^{25,24}+0.1101\times\delta\nu^{26,24}$ with respect to the $^{24}$Mg$^+$ transition frequency. The low natural abundance of $^{25}$Mg reduces the uncertainty of the center-of-gravity frequency to about~2~MHz.

Pickering~et~al.~\cite{Pickering1998,Pickering2000} measured the D$_1$ and D$_2$ transitions in the spectrum of a hollow-cathode discharge lamp. Isotope shifts and hyperfine constants reported in~\cite{Drullinger1980} were used for their line shape model fitted to the isotopically unresolved lines. The isotope shift measurement of the D$_2$ line was taken as an approximation for the D$_1$ line. While the transition frequencies deduced for the $^{24}$Mg$^+$ components are in reasonable agreement with our measurements (see table~\ref{TableResults}), the transition frequencies given for the composite lines significantly deviate (3~$\sigma$) from the results derived from our isotopically resolved measurement (see table~\ref{TableComposite}). We attribute this to discrepancies in the definitions of the composite line center. In~\cite{Pickering1998} the spacing between the $^{24}$Mg components and the center of the composite lines were reported to be 378~MHz~(D$_1$) and 360~MHz~(D$_2$), analyzing the center-of-gravity of isotopically resolved data we find 501~MHz~(D$_1$) and 502~MHz~(D$_2$, Drullinger~\cite{Drullinger1980}:~495~MHz).

Transition frequencies reported in~\cite{Aldenius2006} were acquired by fitting a single Voigt profile to the composite line-shape, the composite line transition frequencies are in reasonable agreement with our center-of-gravity results~(see table~\ref{TableComposite}).

\subsection{Comparison with Isotope Shift Theory \label{section2b}}

\begin{table}[tb]
\caption{Comparison of measured and predicted $^{26-24}$Mg$^+$ isotope shifts. The field shift contributions are included (each +18~(20)~MHz limited by experimental nuclear charge radii~\cite{angeli2004}) and set the uncertainty level for the comparison of our measurement with theoretical investigations. Uncertainties for the SMS contribution are only given in~\cite{berengut2003}.}
\centering
\begin{tabular}{|l|l|l|}
\hline
Reference &~~~~D$_1$ [MHz] &~~~~D$_2$ [MHz] \\
\hline
\hline
Exp. this work~~~~~~~&~~3084.905 (93)~~&~~3087.560 (87)~~\\
\hline
Exp. Drullinger~\cite{Drullinger1980}~&~~~~~~~~~-  &~~3050 (100) \\
\hline
Theory Safronova~\cite{safronova2001}~~&~~3059  &~~3060   \\
\hline
Theory Tupitsyn~\cite{tupitsyn2003}~~&~~3200  &~~~~~~~~~-  \\
\hline
Theory Berengut~\cite{berengut2003}~~&~~3094 (43)  &~~3099 (28)  \\
\hline
Theory Korol~\cite{Korol2007}~~&~~3107  &~~3109  \\
\hline
\end{tabular}
\label{TableShift}
\end{table}

The accuracy of theoretical isotope shifts used to exceed experimental data~\cite{berengut2003,Drullinger1980}, our measurement reverses this situation. Theoretical isotope shifts from recent investigations~\cite{berengut2003,safronova2001,tupitsyn2003,Korol2007} are shown in table~\ref{TableShift}. Different results arise mainly from the determination of the SMS, since differences in the NMS and differences in the field shift due to uncommon field shift constants $F_{\rm{FS}}$ are marginal, compared to the uncertainty introduced by experimental rms nuclear charge radii~\cite{angeli2004}. Only the authors of~\cite{berengut2003} give error bars for their evaluation of the SMS constant {\it k}$_{\rm{SMS}}$, the uncertainty levels introduced by {\it k}$_{\rm{SMS}}$ and the nuclear charge radii are both on the order of 20~MHz. Therefore the considerably more precise value of our measurement can neither set an accurate calibration point for SMS investigations nor give rigorous insight on isotopic differences in the nuclear charge distribution. Yet the results of~\cite{berengut2003} and our measurement are in good agreement. Combining our results with the values given for {\it k}$_{\rm{SMS}}$ and $F_{\rm{FS}}$ from~\cite{berengut2003} allows to deduce the difference in the rms nuclear charge radii. Here we find $\delta\left\langle r^2\right\rangle^{24-26}$=~0.06~(15)~fm$^2$ from the D$_2$ line and $\delta\left\langle r^2\right\rangle^{24-26}$=~0.07~(30)~fm$^2$ from the D$_1$ result, consistent with the tabulated value of 0.14~(16)~fm$^2$~\cite{angeli2004}.

The accuracy of our measurement allows to resolve the $^{26-24}$Mg$^+$ isotope shift of the forbidden 3p$_{3/2}$-3p$_{1/2}$ transition. Since the fine structure splitting is of purely relativistic origin our measurement may be an interesting test for relativistic effects in isotope shift investigations. One result of~\cite{Korol2007} was, that the final relativistic mass shifts for levels of a multiplet appear to be closer to each other than the mass shifts obtained with non relativistic operators, giving rise to a difference in the determined mass shifts of 2.2~MHz and -14.4~MHz for the relativistic and non relativistic treatment respectively (see table~V in~\cite{Korol2007}). Our measured value of 2.66~(13)~MHz for the total isotope shift supports this result.
\\

\section{Summary and Outlook}
\label{section6}

We have accurately measured the $^{24}$Mg and $^{26}$Mg components of the 3s-3p fine structure doublet in Mg$^+$ using a new spectroscopy method for trapped ions in the weak binding limit~\cite{Herrmann2009}. Absolute frequency calibration (relative to the definition of the SI second) allows us to directly deduce the corresponding isotope shifts and fine structure splittings. The measurement accuracy exceeds previous results by more than two orders of magnitude and allows to resolve the isotope shift of the relativistic fine structure splitting for the first time. We deduce transition frequencies for the $^{25}$Mg hyperfine structure and the center-of-mass of an isotope mixture of natural abundance. A further measurement of the $^{25}$Mg component would improve the uncertainty of composite transition frequencies, isotope shifts and hyperfine constants for the 3p$_{1/2}$ and 3p$_{3/2}$ level~\cite{Safronova1998,Sur2005}. In turn it would allow for a more stringent consistency check of theoretical hyperfine structure constants, isotope shift theories and nuclear charge distributions.

Though our measurement may be regarded as exact reference for present and near future astrophysical spectra and theoretical isotope shifts, future applications may require even more accurate calibration. Narrowing down the linewidth of the spectroscopy laser e.g.~by stabilization to an ultra-stable reference cavity would essentially remove the limiting systematic effect of our measurement. The statistical uncertainty of the $^{26}$D$_2$ line was only 14~kHz, indicating that our spectroscopy scheme may allow for further improvement of at least one order of magnitude.

This research was supported by the DFG cluster of excellence
``Munich Centre for Advanced Photonics''.


text:
\section*{INTRODUCTION}

\section*{Predicting folding cores}

In Ref.\ \cite{HeaRBF14} we determined the experimental folding core for the complex between the multifunctional protein cyclophilin A (CypA) and the immunosuppressant drug cyclosporin A (CsA) using hydrogen-deuterium exchange NMR (HDX). We showed that the established theoretical method of rigidity analysis implemented using the {\sc First} software provides a reasonable prediction of the experimental folding core  \citep{HesRTK02,RadB04,RadAIK04}. In addition to this, we generated the {\sc First} folding core (FIR) for unbound CypA, which broadly agreed with the previously published data from HDX experiments. Whilst impressive, the {\sc First} method is not perfect and does not capture the changes in the HDX folding core observed upon ligand binding.
Indeed, FIR is larger in the absence of the ligand, a result which contrasts with what is observed experimentally. That there exists some discrepancy between theory and experiment could be expected since {\sc First} is used to make inferences about protein flexibility based upon a {\it static} crystal input structure only. In addition to flexibility, HDX is dependent upon the surface exposure of the exchanging protons as well as protein motion \citep{Woo93,LiW99}.
With this in mind, we now seek to improve upon the predictive power of the rigidity analysis method by incorporating information on the surface exposure and dynamics of the protein.
We change the hydrogen bond (HB) and hydrophobic tether (HP) networks so that these interactions are absent for surface atoms. With this modification we enhance the flexibility of the protein surface so that the remaining rigid residues may better correlate with the experimental folding core.
Building upon this, we model the prospective dynamics of the protein using coarse-grained simulations. With these techniques we are able to probe the propensity for large-amplitude motions which may only be accessed over long timescales, such as the timescales of our HDX experiments.

\section*{Coarse-grained simulations of protein motion}

The study of protein motion is of fundamental importance in structural biology due to its intimate link with protein function \citep{HenK07}.
The most prevalent method for modeling protein motion computationally is through molecular dynamics (MD) simulations \citep{GohT06}, yielding detailed atomic trajectories \citep{KarP90}.
MD simulations are known to be computationally intensive; yet some of the most biologically relevant motions take place on the ms--s timescale which remains inaccessible to MD at present \citep{WelMHT05}.
Methods of coarse graining may be applied to reduce computational demand and to enable simulation of  larger amplitude protein motion \citep{Cle08}.
The constrained geometric simulation software {\sc Froda} (framework rigidity optimized dynamic algorithm) \label{froda} is a module within {\sc First} which rapidly simulates possible protein motion \citep{WelMHT05}.  It has been used to investigate protein-protein docking problems \citep{JolWHT06} as well as substrate recognition and conformational changes \citep{MacNBC07}. Large-amplitude motions in proteins with hundreds or thousands of residues can be rapidly explored using a desktop computer on a timescale of minutes \citep{JimFRW12} --- an improvement on the order of $10^4$ over MD in terms of computational speed for motions over comparable distances.
With {\sc Froda} we direct motion along a small number of low-frequency normal mode vectors \citep{Dia90,Hin98,SuhS04,PetP06,DobLS08}, the superposition of which well describes large-amplitude protein mobility \citep{KreAWE02,AleLEM05,JimFRW12}.  The resulting potential for {\sc Froda} to access the rarely sampled, high energy protein conformers  to which HDX is sensitive  \citep{SljW13} ensures that it is a natural choice of simulation method for improving our folding core predictions.

\section*{MATERIALS AND METHODS}

In addition to the FIR folding core,  we now introduce four other possible theoretical estimates for the folding core. These include information on the surface exposure and dynamics of the protein in addition to rigidity. The PDB structure 1CWA \citep{MikKPW93} was used for all simulations of protein rigidity and motion, before and after removal of the ligand CsA.
Initially we focus upon the CypA-CsA complex, before applying our analysis to the unbound CypA and comparing our results.

\section*{A modified bond network for FIR$_B$}
\label{sec:HP_settings}

The presence and strength of HBs within the constraint networks were determined based upon donor-acceptor geometry using {\sc First}. Bonds weaker than an energy cutoff value $E_{\mathrm{cut}}$ are excluded from the network. During a rigidity dilution  $E_{\mathrm{cut}}$ is systematically lowered and HBs are consequently removed.\footnote{We emphasize that $E_{\mathrm{cut}}$ should be thought of as an arbitrary parameter and not as directly related to an energy scale for motion --- such as $k_\mathrm{B}T$ --- in the structure \cite{JimFRW12}.}
In addition to the standard bond network defined in this way by {\sc First}, we will also consider a modified network based on the
surface exposure of atoms with HBs and HPs. By demanding that both interacting atoms in a HB or HP constraint are buried within the protein, i.e.\ not exposed on the surface, we enhance the flexibility of the protein surface. The folding core generated in this way consists of residues which are both rigid and buried within the protein.
HPs are indirect, entropy-driven interactions  \citep{FolRD08} thought to contribute significantly to protein folding \citep{Dil90}.
HPs have not always been included in {\sc First} simulations, and their effect on modeling the flexibility and mobility of proteins is still not fully understood. In early papers using {\sc First}, HPs were not included in the bond network \citep{ThoHYK00,JacRKT01}, whereas more recently it is usual practice to include HPs and maintain them throughout rigidity dilution  \citep{ZavLTD04,RadG08}  or to increase their number as $E_{\mathrm{cut}}$ is lowered \citep{RatRG12,PflRKR13}.
In {\sc First}, HPs are modelled as flexible constraints, restricting the separation distance between atoms involved in the interaction, but not the angle between them. In this way, the interacting atoms are permitted to slip relative to each other \citep{HesRTK02,GohKC04}.
HPs between carbon or sulfur atoms are typically included if the distance between these atoms is less than the sum of their van der Waals radii, $r_v$, plus a distance cutoff, $D_{\mathrm{HP}}$. For carbon and sulfur, $r_v = 1.7$~$\mathrm{\AA}$ and $1.8$~$\mathrm{\AA}$ respectively, and $D_{\mathrm{HP}}$ is typically set to 0.25~$\mathrm{\AA}$ \citep{HesRTK02,GohKC04,RadG08}. These distances will then allow to define the burial distance of a given atom (see below).

As a first improvement to FIR, we therefore modified the bond network in {\sc First}, requiring that each HP and HB is formed between two interacting atoms, both of which are buried within the protein. We implemented this bond network and performed a rigidity dilution on the structure 1CWA. As in Ref.\ \cite{HeaRBF14}, the folding core is determined from the dilution plot. In this case, the folding core consists of residues which are rigid \emph{and} buried, and we refer to this measure as FIR$_B$.

\section*{Burial distances}
\label{sec:bur_method}

Let us now define in detail how we distinguish buried and exposed residues.
For each atom in the protein, a sphere of radius $r_v + r_w$ is drawn with the atom in the centre, where $r_v$ is the van der Waals radius of the atom and $r_w = 1.4$~\AA is that of a water molecule; for nitrogen, $r_{v}(N) = 1.5$~\AA\.
The points forming the sphere are then individually checked for contact with the neighboring atoms. If any point on the sphere is not in contact with a neighbor then that is a potential solvent position and the atom is labeled as being exposed. If an atom is not exposed, it is buried and its burial distance, $R$, is the shortest distance to an exposed atom. An exposed atom has $R =0$~\AA.
For a given conformer, we calculated $R$ for each nitrogen atom in the amide backbone and used this  as the value of $R$ for the residue to which it belongs.

\section*{Coarse-grained mobility simulation for FRO}

\label{sec:Fr_method}
\label{sec:Pseudo_method}

The coarse-grained elastic network model implemented using the {\sc ElNemo} software \citep{SuhS04} allows for a rapid and accurate calculation of the normal modes of motion \citep{Tir96,Hin98}.  In Ref.\ \citep{WelJR09} we reported results for directed motion along these normal modes to generate new conformations satisfying the constraints of the {\sc First} bond network.
Following this work, we bias motion along both (parallel and antiparallel) directions of the ten lowest frequency non-trivial normal mode vectors, $m_7$ -- $m_{16}$ with directed and random step sizes both set to 0.01~\AA. These low frequency modes well capture the large scale motion of a protein \citep{JimFRW12}. A total of $2000$ conformers were generated for each mode and direction. Details of the algorithm can be found elsewhere \citep{WelMHT05}. In each of our {\sc Froda} simulations, we use  $E_{\mathrm{cut}} = -2.0$ kcal/mol to determine the HB network. These simulations have been repeated at five different $E_{\mathrm{cut}}$ values ($-0.5, -1.0, -1.5, -2.0$ and $-3.0$~kcal/mol), but the nature of the results is not particularly sensitive to this (data not shown, see \citep{Hea13}).

We carried out simulations of motion using {\sc Froda} and calculated burial distances for all residues in the final (2000$^{th}$) conformer  from each of the {\sc Froda} simulations biased along a low frequency normal mode vector. The residues were classed as being buried ($R > r_v(N)\equiv 1.5$~\AA) or exposed ($R \leq 1.5$~\AA) in each of the conformers. Those which remained buried in more than half of the final conformations formed the {\sc Froda} folding core, which we refer to as FRO.

\section*{Secondary structure estimate SeS and the FIR+FRO combination F+F}

As well as FIR, FIR$_B$ and FRO introduced above, we consider two further theoretical folding cores, SeS and F+F. These are formed of the secondary structural units of helices and sheets (SeS) and the intersection of results between FIR and FRO (F+F). We choose to include SeS because of its simplicity and lack of need for any sophisticated computational analysis. The estimate F+F was included in order to test whether there is a distinct improvement upon combining data from static (rigidity, FIR) and dynamic (mobility, FRO) measures.
In addition, we checked the intersection between FIR$_B$ and FRO. This is missing just one residue Val12 from F+F, but is otherwise identical. Since Val12 is part of the HDX folding core, the  FIR$_B$+FIR folding core makes a slightly weaker prediction than F+F in the present case.

\section*{Quantitative measures for comparing folding cores}

Let  $N_\mathrm{Ex}$ and $N_\mathrm{Th}$ be the number of residues contained in an experimentally and theoretically determined folding core, respectively. Clearly, $N_\mathrm{Ex}=N_\mathrm{Th}$ is a necessary condition for  agreement of experimentally and theoretically estimated folding cores. However, not only the number of residues, but of course the agreement of the specific set of residues in both experimental and theoretical folding cores is most important. In order to capture this, let us define ${\cal T}$ as the number of residues correctly identified by a theoretical prediction of the experimental folding core, i.e.\ ${\cal T}$ is the number of true positives. For a perfect agreement, we have ${\cal T}=N_\mathrm{Th}$ while the expected ${\cal T}$ attained randomly is $N_\mathrm{Th}N_\mathrm{Ex}/N$. Here $N$ denotes the total number of amino acids in the protein \citep{RadB04}.
We can define the \emph{specificity}, $\alpha$, and  \emph{sensitivity}, $\gamma$, of a theoretical folding core prediction as
\begin{equation}
\alpha = \frac{{\cal T}}{N_\mathrm{Th}}
\label{eqn:alpha}
\end{equation}
and
\begin{equation}
\gamma = \frac{{\cal T}}{N_\mathrm{Ex}}.
\label{eqn:gamma}
\end{equation}
The specificity $\alpha$ measures the fraction of residues identified by the theoretical method which are also part of the experimental folding core. The sensitivity $\gamma$ shows the proportion of residues in the experimental core which have been correctly predicted by the theoretical method. A perfect correspondence between theory and experiment, ${\cal T}=N_\mathrm{Th}=N_\mathrm{Ex}$, would yield $\alpha = \gamma = 1$ while for a completely wrong identification, ${\cal T}=0$, we have $\alpha = \gamma = 0$.

Another, previously defined \citep{RadB04}, quantitative measure is the so-called folding core identification \emph{enhancement factor}, $\epsilon$.
This is the ratio of ${\cal T}$ and the number of residues expected to be identified by a random selection,
\begin{equation}
\epsilon = \frac{{\cal T} N}{{N_\mathrm{Th}  N_\mathrm{Ex}}}.
\label{eqn:epsilon}
\end{equation}
\label{eqn:s(AB)}
A theoretical method with random probability of success has $\epsilon = 1$, and when $\epsilon > 1$, the match is better than random. We shall also compare our $\alpha$ and $\gamma$ measures to $\epsilon$.

\section*{RESULTS AND DISCUSSION}

\section*{Effects of the modified bond network and the inclusion of mobility}

We measured the surface exposure of each amide nitrogen in the structure 1CWA, cp.\ Figure \ref{fig:stat_bur}. For the small protein CypA, we find that the majority of the residues are somewhat exposed to the protein surface, with $R<r_v(N)=1.5$~\AA. Since we measure the burial of the amide nitrogen and the length of the N-H bond is $1.5$~\AA, if the H atom is exposed to the surface, then the amide nitrogen has $R \leq 1.5$~\AA.
\begin{figure}[tbp]
\begin{center}
\includegraphics[width=\columnwidth]{Fig-p2-1CWA_stat_burial_line.pdf}
\caption[]{The burial distance $R$ for each residue of CypA in the static structure of the CypA-CsA complex. The blue dashed line at $R = 1.5$~\AA\ is used to separate exposed residues (below the line) from buried residues (above the line). Below the horizontal axis, we schematically show the extent of the FIR$_B$ folding core along the residues via coloured regions. Those residues which are part of FIR$_B$ are coloured red; other rigid residues are coloured otherwise. The thin horizontal black line denotes residues that are flexible.}
\label{fig:stat_bur}
\end{center}
\end{figure}
In Figure~\ref{fig:stat_bur} we also show a superposition of FIR$_B$ along the protein backbone as taken from the rigidity dilution plot. We note that most of the residues that are highly buried within the protein structure ($R\gtrsim 2.0$~\AA) correspond to regions of the protein that are part of the folding core. This shows that the majority of those residues that are buried in the CypA-CsA complex are in fact rigid \emph{and} buried.

Let us now see how our improved FIR$_B$ folding core estimate is different from the FIR folding core of Ref.\ \cite{HeaRBF14}. Figures \ref{fig:HB0_FC} (a) and (b) show a cartoon representation of the residues identified as belonging to the FIR$_B$ folding core superimposed on the 3D structure of the protein. We also indicate in Fig.\ \ref{fig:HB0_FC} (c) the FIR$_B$ folding core along the protein chain as in Figure~\ref{fig:stat_bur}. This is similar to the cartoon representation of the comparison between the HDX folding core and the FIR estimate given in Ref.\ \cite{HeaRBF14}.
\begin{figure}[tbp]
\begin{center}
\includegraphics[width=\columnwidth]{Fig-p2-1CWA_FIR_FIRB.pdf}
\caption[]{(a) {\sc PyMOL} cartoon of 1CWA with FIR$_B$ in blue and the ligand, CsA, in black. Those residues which are part of FIR but not FIR$_B$ are colored red, with the other residues of CypA in green (the colour scheme is similar to the one used in Fig.\ \ref{fig:stat_bur}). (b) The surface of the protein is illustrated as a mesh, with CsA represented as spheres. The image has been rotated from its orientation in (a) by $180^{\circ}$ about the $45^\circ$ diagonal. (c) A linear representation  by residue number of the rigid cluster decomposition at $E_{\mathrm{cut}} = -1.061$~kcal/mol shown as thick colored blocks (rigid) or thin black lines (flexible). The mutually rigid red blocks show the amino acids  which constitute FIR$_B$.}

\label{fig:HB0_FC}
\end{center}
\end{figure}
As expected, FIR$_B$ is smaller than FIR, since there are fewer constraints present throughout the rigidity dilution. The amino acids $12$--$18$ and $41$--$47$ are part of FIR but not part of FIR$_B$. These residues are close to the protein surface on the opposite side of the protein to the ligand binding site as can be seen in Figure \ref{fig:HB0_FC}(b).

Figure \ref{fig:froda_FC}(a) shows a qualitative comparison between FIR$_B$ and the HDX folding core, depicted on the 3D structure of the CypA-CsA complex.
\begin{figure*}[bht]
\begin{center}
\includegraphics[width=\textwidth]{Fig-p2-1CWA_FIRB_FRO.pdf}
\caption[]{Theoretical folding cores (a) FIR$_B$ and (b) FRO for the CypA-CsA complex compared with the HDX folding core. In both cases, residues are colored purple (part of both folding cores),
      blue (HDX folding core only), red (theoretical folding core
      only) or green (not part of either folding core). CsA is
      colored black.}
\label{fig:froda_FC}
\end{center}
\end{figure*}
We find that for most of the residues in the protein, the theoretical and experimental folding cores are in agreement. However, residues 138 -- 143, forming the $\alpha$-helix at the forefront of the image, are part of the HDX folding core, but are \emph{not} part of FIR$_B$. With reference to the representation of FIR$_B$ along the protein chain, as shown in Figure \ref{fig:HB0_FC}(c), we see that these residues are indeed rigid at $E_\mathrm{cut}$ corresponding to FIR$_B$, but they are not part of the largest rigid cluster. In this case, loosening the definition of the {\sc First} folding core to include residues that are rigid independently of the largest rigid cluster would capture this slowly exchanging helix as part of the rigid folding core.

In Figure \ref{fig:froda_FC}(b) we show the {\sc Froda} folding core compared qualitatively with the HDX folding core.
Notably, residues $138$ -- $143$ in the helix are now part of the FRO folding core (but not of FIR$_B$). We emphasise that the results of this mobility-based analysis depend indirectly on the results of the rigidity analysis; {\sc First} identifies a constraint network and {\sc Froda} then explores the motion that is possible within those constraints. Therefore it is possible for a residue to be (i) rigid but exposed --- thus lying in FIR or FIR$_B$ but not in FRO --- or (ii) non-rigid but well protected by burial within the protein --- thus lying in FRO but not in FIR or FIR$_B$.
This suggests that in order to predict the results of the HDX folding core, which depend upon surface exposure as well as dynamics (inherited from flexibility) a combination of the approaches employed by {\sc First} and by {\sc Froda} may be necessary.

In Figure \ref{fig:FCs}, each of the five theoretical folding cores are compared with HDX folding cores for the CypA-CsA complex. In each case, folding core residues along the protein backbone are shown as bold blocks of colour.
\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=\textwidth]{Fig-p2-5methods.pdf}
\caption[Theoretical and experimental folding cores]{Comparison of the experimentally determined HDX folding core for CypA-CsA (blue and shown above and below the other folding cores) with the five folding cores computed using the theoretical approaches outlined in the text and labeled accordingly (from top to bottom: SeS (yellow), FIR (red), FIR$_\mathrm{B}$ (orange), FRO (purple), and F+F (pink)).
In each case the folding core residues in the primary structure of CypA are represented as colored blocks while non-folding core residues are kept white. The residue
      numbers along the protein backbone are indicated horizontally, with thin black vertical lines
      added every ten residues for clarity.}
\label{fig:FCs}
\end{center}
\end{figure*}
FIR, FIR$_B$ and FRO have been discussed above, so we here briefly discuss the remaining folding cores definition SeS and F+F.
We find that SeS consists almost entirely of slowly exchanging amino acids, although there are also many slowly exchanging residues which are not part of a secondary structure unit.
The residues in FRO are also for the most part slowly exchanging.
We note that some parts of the HDX folding core, for example residues $86$--$90$, are poorly predicted by all of the methods. This region corresponds to a flexible, surface-exposed, unstructured region in the crystal structure, indicated by an arrow in Figure \ref{fig:F+F_FC}(a).
\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=\textwidth]{Fig-p2-F+F.pdf}
\caption{Cartoon representation of a comparison of HDX folding cores with F+F for (a) the CypA-CsA complex and (b) unbound CypA. In both cases, residues are colored purple (part of both folding cores), blue (HDX folding core only), red (F+F folding core
      only) or green (not part of either folding core). CsA is
      colored black in (a). The arrow in (a) indicates residues $86$--$90$.}
\label{fig:F+F_FC}
\end{center}
\end{figure*}
Our methods, based around predicting rigid, buried and immobile regions, fall short of identifying this area as slowly exchanging.
Figure \ref{fig:F+F_FC}(b) presents a graphical comparison between the HDX folding core and F+F on the structure of CypA.  We observe that in the absence of the ligand, residues  $86$--$90$ are mostly absent from both theoretical and experimental folding cores. The relatively small number of false positives (red) to false negatives (blue) demonstrate the high ratio of specificity to sensitivity for this method.

\section*{Quantitative analysis of theoretical folding cores}

Table \ref{tab:sAB} shows the number of residues in each theoretical folding core along with the number of true positives ${\cal T}$.  The comparison pictured in Figure \ref{fig:FCs} is made quantitative with the specificity $\alpha$, the sensitivity $\gamma$, and the enhancement factor $\epsilon$, as defined in Equations (\ref{eqn:alpha}), (\ref{eqn:gamma}) and (\ref{eqn:epsilon}), respectively.
\begin{table*}[bt]
\caption[Quantitative measures of agreement between theoretical methods and HDX folding core for CypA-CsA and unbound CypA]{Quantitative measures of agreement between theoretical folding cores and the experimental HDX folding core for (top) CypA-CsA  and (bottom) unbound CypA with $N=165$. The bold numbers indicate the largest selectivity $\alpha$, sensitivity $\gamma$ and enhancement $\epsilon$ for CypA-CsA and unbound CypA. The error estimates for $\alpha$ and $\gamma$ reflect an assumed $\pm 5\%$ variation in $N_\mathrm{Ex}$, $N_\mathrm{Th}$ and ${\cal T}$. The two estimates for the enhancement $\epsilon$ show the variation when using the full set of residues with $N=165$ and when excluding the $6$ proline residues of CypA such that $N=159$.
}
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|c|c|c|}
\hline
& Folding core & $N_\mathrm{Ex}$ & $N_\mathrm{Th}$ & ${\cal T}$  & $\alpha$ & $\gamma$ & $\epsilon_{165}$  & $\epsilon_{159}$ \\
\hline
\hline
& SeS      & 80 &	75	&	59	 &	$0.79\pm 0.04$ 	&	{$\mathbf{0.74}\pm 0.04$} &	$1.62$	& $1.56$	\\
& FIR	& 80 &	86	&	57	 &	$0.66\pm 0.04$ 	&	$0.71\pm 0.04$ &	$1.37$	& $1.32$	\\
CypA-CsA
& FIR$_B$ & 80 &	72 	&	53	 &	$0.74\pm 0.04$ 	&	$0.66\pm 0.04$  &	$1.52$	& $1.46$ \\
& FRO	& 80 &	80	&	56	 &	$0.70\pm 0.04$ 	&	$0.70\pm 0.04$ &	$1.44$	& $1.39$	\\
& F+F	& 80 &	54	&	44	 &	$\mathbf{0.81}\pm 0.05$ & $0.55\pm 0.03$ &	$\mathbf{1.68}$	 & $\mathbf{1.62}$	\\
\hline
& SeS	& 73 & 75	& 58  & $0.77\pm 0.04$ & $0.79\pm 0.04$ & $1.75$ & $1.68$  \\
& FIR 	& 73 & 112	& 66  & $0.59\pm 0.03$ & $\mathbf{0.90}\pm 0.05$ & $1.33$ & $1.28$ \\
unbound CypA
& FIR$_B$ & 73 & 72 	& 52 & $0.72\pm 0.04$ & $0.71\pm 0.04$ & $1.63$ & $1.57$  \\
& FRO 	& 73 & 76 	& 54  & $0.71\pm 0.04$& $0.74\pm 0.04$  & $1.61$ & $1.54$\\
& F+F 	& 73 & 59 	& 51 & $\mathbf{0.86}\pm 0.05$ & $0.70\pm 0.04$ & $\mathbf{1.95}$ & $\mathbf{1.88}$  \\
\hline
\end{tabular}
\end{center}
\label{tab:sAB}
\end{table*}%
Each of the five theoretical folding cores are significant improvements upon random selections, as demonstrated by the values of $\epsilon$, the lowest of which is $1.33$. The four theoretical methods introduced in this paper all outscore FIR with respect to this enhancement measure, with F+F scoring highest for both the CypA-CsA complex and the unbound CypA, with SeS second. SeS is highly selective and sensitive in both cases, with specificity $\alpha \geq 0.77$ and sensitivity $\gamma \geq 0.74$. Unfortunately, the simple SeS method of folding core selection is clearly inappropriate for predicting changes upon ligand binding, since by definition it does not change.
As discussed in Ref.\ \cite{HeaRBF14}, FIR is sensitive to the bond network and its increase in size upon ligand removal is in disagreement with the experimental results, which show a larger folding core in the presence of CsA. The most accurate method for CypA-CsA and unbound CypA is F+F, with $\alpha = 0.81$ and $0.86$.  This high specificity comes at the cost of lower sensitivity, meaning that it predicts with high accuracy but does not capture all of the HDX folding core. A larger proportion of the HDX folding core needs to be accurately determined by this method in order to successfully capture the impact of ligand binding on HDX.
We have also estimated the variation in $\alpha$, $\gamma$ and $\epsilon$ in Table \ref{tab:sAB}, assuming a $\pm 5\%$ error in $N_\mathrm{Ex}$, $N_\mathrm{Th}$ and ${\cal T}$. This leads to some variation in the results for $\alpha$ and $\gamma$ and allows us to estimate the accuracy of the numbers quoted here. In addition, due to their lack of an amide proton, the six proline residues in CypA do not appear in the $^1$H-$^{15}$N HSQC spectrum and so cannot be part of the HDX folding core. For this reason, we also take $N=159$ when calculating $\epsilon$. We find that in most cases the largest values of $\epsilon$, $\alpha$ and $\gamma$ remain significantly larger. This shows that our findings are indeed robust.

\section*{Effect of ligand binding on folding cores}

As demonstrated in Ref.\ \cite{HeaRBF14}, for CypA there are only small changes to the HDX folding core upon ligand binding, and these are not captured using rigidity analysis alone. Indeed, only seven residues are part of the HDX folding core for the complex but not for the unbound protein.
Upon ligand removal, FIR$_B$ changes by a single residue, which is not one of these seven. As discussed in Ref.\ \cite{HeaRBF14}, the FIR changes unexpectedly upon ligand binding. These two folding cores solely derived using {\sc First} do not reflect the effect observed experimentally. We see a larger effect of ligand binding on FRO.
Figure \ref{fig:burial_abs} shows how ligand binding affects the burial distance $R$. We calculated the average $R$ for each residue from the final conformers of the {\sc Froda} simulations, and plotted the absolute difference, $\Delta R$, between these values.
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\columnwidth]{Fig-p2-abs_diff_bur.pdf}
\caption[Ligand binding affects $R$ for binding site residues]{The difference in burial distance $\Delta R$ (black bars) between the CypA-CsA complex and the unbound protein plotted for each amino acid. The arrows and vertical dotted lines indicate the binding residues. As in Figures \ref{fig:HB0_FC} and \ref{fig:stat_bur}, FIR$_B$ is also shown.}
\label{fig:burial_abs}
\end{center}
\end{figure}
Of the eleven residues for which $\Delta R > 0.5$ $\mathrm{\AA}$, five are binding site residues and none are further than two residues along the backbone from a binding site residue.
We also calculated $\Delta R$ for each residue  from {\sc Froda} simulations carried  out at $E_{\mathrm{cut}}$ = $-0.5$, $-1.0$, $-1.5$ and $-3.0$ kcal/mol.
In each case, the residues which are part of, or close to, the binding site have the largest $\Delta R$ (data not shown).
 The HDX folding core changes in sections of the protein more distant from the binding site, and so tracking the protein surface in {\sc Froda} does not capture this effect.

\section*{CONCLUSION}
We have introduced and discussed rapid computational methods for adding information on surface exposure and protein dynamics into a rigidity-based folding core definition. These improve upon the {\sc First} folding core for predicting the HDX folding cores of both CypA-CsA and unliganded CypA. In addition, we show that they are less sensitive to the erroneous increase in folding core size upon ligand removal observed using {\sc First} with the standard bond network. On balance, F+F, which combines rigidity and motion, appears as the best choice for a computational determination of the folding core from a protein's structural information alone.
Nevertheless, no method achieves a perfect score of $\alpha=\gamma=1$, neither for the uncomplexed CypA nor for the CypA-CsA complex. This insensitivity to ligand binding is disappointing but unsurprising due to the scale of the challenge --- for CypA-CsA and in general.

We have sought to obtain a practical balance between the computational demands of our predictions and their utility. The requirement for large-scale motions make more computationally sophisticated methods such as MD inappropriate. Methods using simpler input data than ours, making sequence-based predictions of surface exposure and then in turn predictions of HDX folding cores, are currently unsatisfactory. We feel that although balance between sophistication of input data and method and the accuracy of the prediction is yet to be struck, it may be achieved through coarse-grained approaches. The {\sc froda} alternative {\sc frodan} has more flexibility of hydrogen bond and hydrophobic contact geometry built in to its bond networks. A comparison between the predictions made using {\sc frodan} instead of {\sc froda} would be a valuable extension to this work.

Our overall strategy is based on using the results from HDX experiments as the benchmark against which theoretical  results are compared. The measures $\alpha$, $\gamma$ and $\epsilon$  score most highly when the match of a theoretical folding core with the HDX result is perfect. However, the HDX data may contain errors themselves and this additional source of variation did not play a role in our judgements of the success of the computational approaches.
Furthermore, for the CypA-CsA complex, the protein is large relative to the ligand and the ligands effect on HDX is rather small.
We expect that an application of our methods to a larger protein that exhibits a more significant conformational change upon ligand binding will also give a clearer change in the values of $\alpha$, $\gamma$ and $\epsilon$.

\section*{ACKNOWLEDGMENTS}

\ack{We thank C.\ Blindauer for help with the experimental NMR aspects of this work. The CypA plasmid was kindly provided by G.\ Fischer from the Max Planck Research Unit for Enzymology of Protein Folding in Halle, Germany. We gratefully acknowledge funding from the EPSRC Life Sciences Interface programme (MOAC DTC EP/F500378/1). JWH thanks the Institute of Advanced Study in Warwick for an Early Career Fellowship; SAW is grateful to the Leverhulme Trust for an Early Career Fellowship.}\vspace*{6pt}

\bibliographystyle{biophysj}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

The observation of natural phenomena and of
laboratory experiments provides a wide spectrum of experimental results
showing distributions of data that deviate from exponential decay as
predicted for Boltzmann behavior~\cite{davis}. It was the goal of
non-extensive statistical mechanics developed originally by Tsallis~\cite%
{tsallis88} to offer a new approach to explain the non-Boltzmann behavior of
non-equilibrium systems \footnote{We note that other approaches exist such
as super-statistics \cite{beck-cohen} and the fractional and nonlinear
Fokker-Planck equations  (see e.g. \cite{lutsko-boon} and references therein)}.
More precisely the primary motivation for
non-extensive thermodynamics is as a way to understand deformed exponential
distributions (such as $q$-exponentials exhibiting long tails when $q>1$)
found empirically in many areas of physics and other scientific disciplines~%
\cite{swinney-tsallis}. The interest raised by this new approach has grown
over the years and has produced an abundant literature~\cite{tsallis09}
reflecting new theoretical developments and a considerable number of
applications to subjects as diverse as defect turbulence, energy
distribution in cosmic rays, earthquake magnitude value distributions and
velocity distributions in micro-organism populations or distributions of
financial market data. Parallel to these developments, critical analyses
were presented questioning the merits of the non-extensive theory~\cite%
{opponents} and thereby of its applications as well. These criticisms raised
questions that often gave rise to ontological conflicts~\cite{conflict}.
Most of the criticisms are based on phenomenological analyses and
thermodynamic arguments questioning the compatibility of the theory with
classical statistical thermodynamics. More recently, the non-extensive
theory was also re-examined on the basis of analyses demonstrating the
necessary {{discreteness}} of systems to which the theory applies~\cite{abe_10}
and the limits of validity of the non-extensive formalism for a Hamiltonian
system, the $q$-ideal gas, a model system of independent quasi-particles
where the interactions are incorporated in the particles properties~{{\cite%
{abe_99, boon_lutsko_11}}}. Here we merely adopt the viewpoint of analytical rigor to
establish the validity limits of the non-extensive formalism for the general
class of classical Hamiltonian systems with continuous variables and
consequently of the class of physical systems to which the non-extensive
interpretation applies.

In the formalism, the $q$-exponential distributions arise as the result of a
variational calculation maximizing the generalized entropy, the $q$-entropy,
under the constraints of normalizability of the distribution function and of
a prescribed average energy~\cite{tsallis09}. The goal of the present work
is to investigate the exact form of the distribution so derived for
classical Hamiltonian systems from both the usual Tsallis entropy~\cite%
{tsallis88} and the homogeneous entropy~\cite{lutsko_boon_grosfils}. The
conclusion is found to be the same in both cases: the theory is only
consistent in the thermodynamic limit for $0\leq q\leq1$. For finite systems
of $N$ particles, the upper limit is $1+a/N$ for $a\sim{\mathcal{O}}(1)$.
This means (i)~that the non-extensive thermodynamics formalism cannot be
used, at least in any straightforward way, to explain phenomena for which one
observes that $q$ takes a value $q>1$ {{(corresponding to asymptotically power-law
distributions)}} i.e. when the distributions exhibit extended exponential
forms with long tails, and (ii)~that, within the validity domain, the
distribution has finite support, thus implying that configurations with e.g.
energy above some limit have zero probability, a strange situation for
systems with typical molecular potentials which are steeply repulsive at
small distances. We treat successively the case according to the development
based on the Tsallis entropy~$S_{q}$ and the case based on the homogeneous
entropy~$S^{H}_{q}$. {{The reason for examining the two cases is the criterion of
stability against perturbations of the probability distribution function, or
{\it {Lesche stability}}~\cite{lesche_abe}}}  : the
homogeneous entropy was shown to be Lesche stable while the Tsallis entropy
is not~\cite{lutsko_boon_grosfils}. We also present the results for the generalized
($q$-dependent) thermodynamic properties of Hamiltonian systems in both cases.

\section{Tsallis entropy}

 Non-extensive statistical mechanics is developed
on the basis of three axioms: (i)~the $q$-entropy for systems with
continuous variables is given by~\cite{tsallis09}
\begin{equation}
S_{q}\,=\,k_{B}\,\frac{1-{K\int \rho ^{q}\left( \Gamma \right) d\Gamma }}{q-1%
}\,,  \label{STq}
\end{equation}%
where $\Gamma$ is the phase space variable and
$K$ must be a quantity with the dimensions of $\left[ \Gamma \right]
^{q-1}$, i.e. $K=\hbar ^{DN\left( q-1\right) }$
($D$ denotes the dimension of the system and $N$ the number of particles)
and the classical Boltzmann-Gibbs entropy is retrieved in the limit $q\rightarrow 1$;
(ii)~the distribution function $\rho \left( \Gamma \right) $ is slaved to the
normalization condition
\begin{equation}
1\,=\,\int \rho \left( \Gamma \right) d\Gamma \,;  \label{Norm}
\end{equation}%
(iii)~the internal energy is measured as
\begin{equation}
U\,=\,{\int P_{q}\left( \Gamma \right) H\,d\Gamma }\,=\,\frac{\int \rho
^{q}\left( \Gamma \right) H\,d\Gamma }{\int \rho ^{q}\left( \Gamma \right)
d\Gamma }\,,  \label{Uesc}
\end{equation}%
where $P_{q}\left( \Gamma \right) $ is the escort probability distribution~%
\cite{beck_schlogl} which is the actual probability measure. We consider
Hamiltonian systems with $H=T+V$ where $T$ is the classical kinetic energy.
The distribution function $\rho \left( \Gamma \right) $ is obtained by
maximizing the $q$-entropy subject to the normalization (\ref{Norm}) and
average energy (\ref{Uesc}).
Introducing the Lagrange multipliers $\bar{\alpha}$ and $\bar{\gamma}$, the variational method  leads to
\begin{equation}
0 = \frac{k_{B}Kq}{1-q}\rho^{q-1}(\Gamma)-\bar{\alpha} - q\bar{\gamma}\frac{(H-U)}{\int \rho ^{q}\left( \Gamma \right)d\Gamma }\rho^{q-1}(\Gamma) \,,
\label{vareq}
\end{equation}
which is solved to yield
\begin{equation}
\rho \left( \Gamma \right) =\frac{\left( 1-\left( 1-q\right) \gamma \left(
H-U\right) \right) ^{\frac{1}{1-q}}}{\displaystyle\int \left( 1-\left(
1-q\right) \gamma \left( H-U\right) \right) ^{\frac{1}{1-q}}d\Gamma }\,,
\label{rhoT1}
\end{equation}%
with  $\gamma = \frac{\bar{\gamma}}{K\int \rho ^{q}\left( \Gamma \right)d\Gamma }$ and
where the normalization condition (\ref{Norm}) has been used to eliminate the
multiplier $\bar{\alpha}$. The other, $\gamma $, is determined from the energy
constraint (\ref{Uesc}); using \Ex{rhoT1} in (\ref{Uesc}), we have
\begin{equation}
0=\frac{\displaystyle \int \left( 1- \left( 1-q\right)\,{\gamma}\, %
\left( H-U\right)\right) ^{\frac{q}{1-q}}\left(
H-U\right) d\Gamma }{\displaystyle \int \left( 1- \left( 1-q\right)\,{\gamma}\, %
\left( H-U\right)\right)^{\frac{q}{1-q}}d\Gamma } \,,
\end{equation}%
which by multiplying the numerator by $(1-q)\gamma$ and adding and subtracting one gives
\begin{eqnarray}
{\int \left( 1-\left( 1-q\right) \gamma \,\left( H-U\right) \right) ^{\frac{q%
}{1-q}}d\Gamma } =
{\int \left( 1-\left( 1-q\right) \gamma \,\left( H-U\right)
\right) ^{\frac{1}{1-q}}d\Gamma }\,.  \label{UT}
\end{eqnarray}%
Note that the sign of the factor $f(\Gamma )\equiv 1-\left( 1-q\right)
\gamma \left( H-U\right) $ is, so far, arbitrary so that we allow for the
cancellation of a (possible imaginary) factor throughout these equations.
However, this is only possible in the expression for the distribution
\Ex{rhoT1} if the sign of $f(\Gamma )$ is independent of $\Gamma $. Given
this fact, the sign can be fixed by making the substitution $f(\Gamma
)=s\left|f(\Gamma )\right|$ in the energy constraint which requires that $s^{\frac{q}{%
1-q}}=s^{\frac{1}{1-q}}$ so that $s=1$. Since the kinetic energy is in
principle unbounded, the constant sign of the argument means that the
distribution may have to be restricted to some subset of phase space so that
it should be written as
\begin{equation}
\rho \left( \Gamma \right) =\frac{\left( 1- \left( 1-q\right) \gamma %
\left( H-U\right) \right)^{\frac{1}{1-q}} \Theta\left( 1- \left( 1-q\right) \gamma %
\left( H-U\right) \right)}%
{\displaystyle \int \left( 1-\left( 1-q\right) \gamma%
\left( H-U\right)\right)^{\frac{1}{1-q}} \Theta\left( 1- \left( 1-q\right) \gamma %
\left( H-U\right) \right)d\Gamma}\,,
\label{rhoT11}
\end{equation}%
where $\Theta (x)$ is the step function which is one for $x>0$
and zero otherwise. We note that the introduction of the step function is in fact a redefinition of the variational problem in the sense that we have replaced $\rho(\Gamma)$ with $\bar{\rho}(\Gamma)\Theta(f(\Gamma))$ in Eq.(\ref{STq}-\ref{Uesc}) and then maximized with respect to  $\bar{\rho}(\Gamma)$. If we do not do this, then there is simply no solution to the variational problem which is real and non-negative for all $\Gamma$.
To see this, we note from the variational equation (\ref{vareq})
that if $\rho(\Gamma)$ vanishes for some value of $\Gamma$ then it necessarily follows that $\bar{\alpha} = 0$, which is untenable since $\bar{\alpha}$ will generally assume a non-zero value due to the normalization condition. Therefore, we can only restrict the support of the distribution by redefining the variation problem.
We now turn to the determination of the possible values for~$q$.

{\it{Case:} $q<1$}. The exponent occuring in the distribution is
$\frac{1}{1-q}>0$ so that if $\gamma>0$, then the distribution must have finite
support since for some sufficiently large value of $T$, the argument, %
$f(\Gamma)$, becomes zero, and from (\ref{rhoT1}) we have
\begin{equation}
\rho\left( \Gamma\right) =\frac{\exp_{q}\left( -\gamma\left( H-U\right)
\right) }{\int\exp_{q}\left( -\gamma\left( H-U\right) \right) d\Gamma }\,,
\label{rhoTq}
\end{equation}
where {{$\exp_q(x) \equiv (1+(1-q)x)^{\frac{1}{1-q}}\Theta(1+(1-q)x)$ }}is the
$q$-exponential function.

If $\gamma<0$, then $f(\Gamma)$ is always positive, and the distribution
function never goes to zero so that there can be no finite support. This
however leads to another problem since $f(\Gamma)$ is unbounded as the
kinetic energy increases which, in turn, means that the integral of $%
f(\Gamma)$ over momenta will diverge so that the distribution cannot be
normalized. We conclude that $\gamma<0$ is not allowed. We remark that it
might be thought that one could introduce a limit on the kinetic energy, but
this is not in keeping with the proposal that the nonextensive distribution
be a generalization of the canonical distribution which describes {\it{open
systems}} in contact with a reservoir.

{\it{Case:} $q>1$}. The exponent being then $\frac{1}{1-q}<0$,
we write the distribution function as%
\begin{equation}
\rho \left( \Gamma \right) =\frac{\left( 1+ \left\vert 1-q\right\vert \gamma%
 \left( H-U\right) \right) ^{-\left\vert \frac{1}{1-q}\right\vert } %
\Theta \left(1+ \left\vert 1-q\right\vert \gamma%
\left( H-U\right) \right)}%
{\displaystyle \int \left( 1+ \left\vert 1-q\right\vert \gamma \left( H-U\right) \right) ^{-\left\vert \frac{1}{1-q}\right\vert } %
\Theta \left(1+ \left\vert 1-q\right\vert \gamma%
\left( H-U\right) \right) d\Gamma } \,.
\end{equation}%
If $\gamma >0$, the distribution will be normalizable provided
that $\left( \sum p_{i}^{2}\right) ^{-\left\vert \frac{1}{1-q}\right\vert }$
is integrable over $d^{ND}p$ (for large $p$) which is to say that $\left(
P^{2}\right) ^{-\left\vert \frac{1}{1-q}\right\vert }P^{ND-1}dP$ be
integrable for $P\rightarrow \infty $; this means we must have $%
-1>ND-1-2\left\vert \frac{1}{1-q}\right\vert $,
or $1<q<1+\frac{2}{ND}$
which condition reduces to the classical Boltzmann result ($q=1$)
in the thermodynamic limit~\footnote{Note that $\rho^{q}\left( \Gamma\right) $ must also be
integrable which imposes that $-1 > ND-1-\frac{2q}{q-1}$ or $q < 1+\frac {2}{%
ND-2}$, but this condition is weaker than the constraint
$1<q<1+\frac{2}{ND}$.}.

The other possibility is $\gamma <0$. We then write the distribution as%
\begin{equation}
\rho \left( \Gamma \right) =\frac{\left( 1-\left\vert \gamma \right\vert
\left\vert 1-q\right\vert \left( H-U\right) \right)%
^{-\left\vert \frac{1}{1-q}\right\vert} %
\Theta \left( 1-\left\vert \gamma \right\vert \left\vert 1-q\right\vert
\left( H-U\right) \right) }%
{\displaystyle \int \left( 1-\left\vert \gamma \right\vert \left\vert
1-q\right\vert \left( H-U\right) %
\right) ^{-\left\vert \frac{1}{1-q}\right\vert }\Theta
\left( 1-\left\vert \gamma \right\vert \left\vert 1-q\right\vert %
\left( H-U\right) \right)
d\Gamma }
\end{equation}%
When integrated over momenta, this expression would show a
singularity at
$T=\frac{1}{\left\vert \gamma \right\vert \left\vert 1-q\right\vert }+U-V$,
unless \\
$\int {\left( X-\sum_{i}p_{i}^{2}\right) ^{-\left\vert \frac{1}{1-q}\right\vert }}dp^{3N}
 \sim \int_{0}^{\sqrt{X}}{\left( X-P^{2}\right)^{\left\vert \frac{-1}{1-q}\right\vert }}P^{3ND-1}dP
\sim \int_{0}^{X}{\left(X-Y\right) ^{\left\vert \frac{-1}{1-q}\right\vert }}Y^{\frac{3ND-2}{2}}dY$
is finite; this requires that $1-\left\vert \frac{1}{1-q}\right\vert =1-%
\frac{1}{q-1}>0$, i.e. $q>2$ . But $\rho ^{q}\left( \Gamma \right) $ must
also be integrable and in order that the singularity be integrable imposes $%
1-\frac{q}{q-1}>0$ which is incompatible  with the  condition $q>2$.

Furthermore the maximization condition demands that the second derivative of the
$q$-entropy \Ex{STq} be $- k_B\,K\, q \,\rho ^{q - 2} < 0$, which is satisfied if $q > 0$.
So in summary, the distribution function exists  in the thermodynamic limit for $\gamma> 0$
and  $0\leq q \leq 1$, and  $\rho\left( \Gamma\right) $ has the form of a $q$-exponential
with finite support.

\section{Homogeneous entropy}

The homogeneous entropy was proposed as an alternative to the Tsallis entropy because it has various desirable properties that the Tsallis entropy does not share such as being Lesche stable and giving positive-definite heat capacities \cite{lutsko_boon_grosfils}. It is therefore interesting to ask whether it is subject to the same limitations as found for the Tsallis entropy. The analysis follows essentially the same lines as above except that in this case $\rho\left(\Gamma\right) $ is the physical probability~\cite{lutsko_boon_grosfils} and
so the energy is computed with the normal average. The homogeneous entropy
reads~\cite{lutsko_boon_grosfils}
\begin{equation}
S_{q}^{H}=k_{B}\,\frac{1-\left( K\int\rho^{\frac{1}{q}}\left( \Gamma\right)
d\Gamma\right) ^{q}}{1-q}\,,
\label{SHq}
\end{equation}
where $K$ is a quantity with the dimensions $\left[ \Gamma\right] ^{\frac{1-q%
}{q}}$, i.e. $K=\hbar^{DN\,\frac{1-q}{q}}$, and the normalization and energy
constraints are
\begin{equation*}
1=\int\rho\left( \Gamma\right) d\Gamma\;\;\;\;;\;\;\;\;U={\int\rho\left(
\Gamma\right) H\,d\Gamma}\,.
\end{equation*}
Following a similar analysis as given in the previous section,
the conclusions are that the condition for normalizability is $0\leq q\leq1$ and
$\gamma>0$ in which case the distribution function reads

\begin{equation}
\rho\left( \Gamma\right) \,=\,\frac{\left( \exp_{q}\left( -\,\gamma \,\left(
H-U\right) \right) \right) ^{q}}{\int\,d\Gamma\,\left( \exp _{q}\left(
\,-\gamma\,\left( H-U\right) \right) \right) ^{q}}\;.
 \label{rhoHq}
\end{equation}

Here
$\gamma = \frac{\bar{\gamma}}{\left(K \int \rho ^{\frac{1}{q}} \left( \Gamma \right) d\Gamma \right)^q}$
where $\bar{\gamma}$ is the Lagrange multiplier used to fix the average energy
and $\rho\left( \Gamma\right)$ has finite support.

\section{Thermodynamic properties}

Given that the formalism is constructed
solely on the basis of the three axioms \Ex {STq}, \Ex{Norm} and \Ex{Uesc},
the consistent way to define the thermodynamic temperature is through the
thermodynamic definition ${\partial S}/{\partial U}= 1/T$.
Considering $q<1$ and  $\gamma >0$,
we obtain
from the Tsallis entropy~\Ex{STq}
\begin{eqnarray}
{T^T_q}=\frac{1}{ k_B\,\gamma}%
\left( K^{\frac{1}{1-q}} \int   \exp_q\left(-\gamma%
\left(H-U\right)\right) d\Gamma \right)^{q-1} \,,
\label{TqT}
\end{eqnarray}%
and from the homogeneous entropy~\Ex{SHq}
\begin{eqnarray}
T^H_q \,=\,\frac{q}{k_B\,\gamma}
\left( K^{\frac{q}{q-1}} \int \exp_q \left( - \frac{\gamma}{q}\left(H-U\,\right)\right)%
d\Gamma \right)^{1-q}  \label{THq}
\end{eqnarray}
With the explicit expressions of ${\gamma}$ (see \Ex{rhoTq} and \Ex{rhoHq}), \Ex{THq}
gives ${\bar{\gamma}}\,=\,{1}/{(k_B\,T^H_q)}$, the analog of the classical expression, whereas
the equivalent relation for the Tsallis temperature is only obtained in the limit $q \rightarrow 1$.

Correspondingly, the expressions for the specific heat  $C_{V}=\left( \frac{\partial U}{\partial T_q}\right)$
are given by
\begin{eqnarray}
C^T_{V}=
\frac{{\beta^T}/{\gamma}}{\left[ \frac{1}{q}\left(\frac{\beta^T}{\gamma}\right)^{4}\frac{1}{K^{2}%
\int\rho^{2q-1}\left(  \Gamma\right)  \left(  \beta^T\left(  H-U\right)
\right)^{2}d\Gamma}-2\left(  1-q\right)  \right]}
\end{eqnarray}
with the classical notation $\beta^T = 1/(k_BT^T_q)$, and by
\begin{equation}
C^H_{V} =  q^{\frac{1}{q}} \left( {\frac{\beta^H}{\gamma}} \right)^{\frac{1-q}{q}} \, {\frac{k_{B}}{K}}%
\int\rho^{\frac{2q-1}{q}}\left(  \Gamma\right)  \left(  \beta^H\left(H-U\right)  \right)  ^{2}d\Gamma
\label{CvH}
\end{equation}
where $\beta^H = 1/(k_BT^H_q)$, or
\begin{equation}
C^H_{V} =  {\frac{1}{k_{B}\,(T_q^H)^2}}%
\int d\Gamma \rho \left(  \Gamma\right) \left(H-U\right)^{2}\,\mathsf{C}_q(\Gamma)
\label{CvH1}
\end{equation}
with $\mathsf{C}_q(\Gamma) = q^{\frac{1}{q}}\left(\frac{\beta^H}{\gamma}\right)^{\frac{1-q}{q}}%
\frac{\rho^{\frac{q-1}{q}}\left(  \Gamma\right)}{K}$. \Ex{CvH1} is the generalization of the classical
expression of the specific heat given in terms of the energy fluctuations:
$C_{V} =  {\langle (\Delta E)^2\rangle}/({k_{B}\,T^2})$.
The thermodynamic temperatures are both positive while the generalized specific heat is only
positive-definite in the case of the homogeneous entropy.\footnote{{It was indeed shown that in the
Tsallis formulation the specific heat can be negative~\cite{abe_99}.}}

\section{Concluding comments}

Non-exponential distributions are widely observed in nature. Non-extensive
thermodynamics was motivated, in part, as a means of explaining the origin
of such distributions which arise naturally as a result of maximizing the
generalized entropy with the usual constraints of the normalization of the
distribution and of fixed average energy. We have shown that when this
procedure is applied to Hamiltonian systems, the resulting distributions
only exist for the restricted range of $0<q<1+{\mathcal O}\left( \frac{1}{N}\right) $.
Since the so-called "fat-tailed" distributions correspond to $q>1$, this
means that generalized thermodynamics cannot be seen as an explanation of
their occurrence for these systems. The problem with larger values of $q$ has
to do with the existence of the integrals over momenta due to the
unboundedness of the kinetic energy. One way around this would be to
redefine the formalism so as to restrict the momenta a priori by making the
ansatz ${\rho\left( \Gamma\right) =\Theta}\left( T_{\ast}-T\right) {\rho}%
_{\ast}{\left( \Gamma\right) }$, for some fixed positive number $T_{\ast}$,
throughout the variational problem and maximizing to determine ${\rho}_{\ast}%
{\left( \Gamma\right) }$. However, this is obviously quite artificial and
{\it ad hoc} since, for example, one could replace the step function by any
function of $T$ that goes to zero sufficiently quickly as $T$ grows.
This suggests the more straightforward conclusion that in the case of classical
Hamiltonian systems, nonextensive thermodynamics does not provide a simple,
natural explanation of distributions with fat tails.

\acknowledgments
This work was {partly} supported by the European Space Agency under contract
number ESA AO-2004-070.

\bigskip