\section{Introduction}

Astronomical observations performed using ground-based telescopes are affected by wavelength-dependent atmospheric dispersion when taken at
a zenithal angle different from zero. The atmospheric dispersion is due to the variation with wavelength of the refractive index of the atmosphere.
By definition, an atmospheric dispersion model computes the refractive index of dry air at sea level as function of wavelength, n($\lambda$).
Since observatories are usually located at higher altitudes, the refractive index should be corrected for temperature, pressure, and relative humidity (RH).
From \cite{Smart1931}, we have that the atmospheric dispersion, in milli-arcseconds (mas) in the sky, is given by:
\begin{align}
   \Delta R (\lambda) &= R(\lambda) - R(\lambda_{\rm ref}) \nonumber \\
   \Delta R (\lambda) &\approx 206265 \left[n(\lambda) - n(\lambda_{\rm ref}) \right] \times \tan Z,
   \label{eq:dispersion}
\end{align}
where  $\lambda_{\rm ref}$ is the reference wavelength, and $Z$ is the zenithal angle of observation. \\
The problem of atmospheric dispersion is particularly serious for high-resolution spectrographs aiming at sub-pixel RV precision. In order to detect Earth-like planets orbiting Sun-like stars, a RV precision of around 10 cm s$^{-1}$ must be achieved. With the objective of developing high-resolution spectrographs that can reach such a RV precision \citep{Fischer2016}, several instrumentation challenges need to be solved including a correction of atmospheric dispersion variation down to this level of positioning accuracy. In order to reach the highest reproductibility in collected spectra and RV properties, an atmospheric dispersion corrector \citep[ADC,][]{Avila1997} is thus mandatory.
\\
Even assuming a perfect telescope centering and guiding at a given reference wavelength, atmospheric dispersion will shift the barycenter of the image as a function of wavelength \citep{Cohen1988}. This will strongly affect the light collected by a spectrograph, especially in high-resolution spectrographs aiming at high-fidelity spectra \citep{Blackman2019}. The loss of light, due to the elongated PSF that can be bigger than the size of the fiber, as a function of wavelength introduces a wavelength dependent photon-noise that becomes particularly stronger at shorter wavelengths where the dispersion is stronger, reducing the ability to measure efficiently the line's position and its properties. The importance of increasing the collected light at shorter wavelength comes from two aspects:
\begin{enumerate}[label=(\roman*)]
\item the first from the fact that the efficiency of many optical components is lower at shorter wavelengths;
\item while the second is from the fact that around 380 nm, there's an abundance of CaII lines that are used as stellar activity indicators \citep{Santos2008}.
\end{enumerate}
In addition, imperfect atmospheric dispersion will change the continuum slope \citep{Pepe2008} which will alter the weight of spectral lines and influence the mean radial velocity (RV) value. \\

According to \cite{Blackman2019}, the ADC coupled with the injection on the fiber (controlled by a guiding camera) might affect the observations and introduce chromatic changes in the recorded spectra. To our knowledge, an exhaustive analysis of the problem in the visible range of the spectrum has never been performed. \cite{Bechter2018} estimated a RV error of 6 cm s$^{-1}$ for the combined adaptive optics (AO) and ADC correction in the near-infrared without mentioning the level of residuals. Also, \cite{Halverson2016} associated an error of 6.9 cm s$^{-1}$ from an ADC with 100 mas of residuals.
\\
While designing an ADC and setting the requirements on the residuals, one should pay attention to not only the effect of atmospheric dispersion residuals on RV, but also on the flux loss associated. The influence of optical glass melt data is also important. In fact, different batches from the manufacturer will end up with different Sellmeier coefficients \footnote{These are the coefficients of the Sellmeier equation that establish the relation between the refractive index and wavelength for a particular transparent medium; it is used to determine the dispersion of light in the medium.} due to melt procedures \citep{Langenbach1999}. The design of an ADC is based on catalog data that is not optimized to the melt data. From an instrumentation point of view, setting the requirements on the ADC to be so demanding (< 50 mas peak-to-valley (PV), maximum amount of residuals within the wavelength range of interest) is becoming more and more difficult. Not only the fact that finding glasses that can return such residuals, and the fact that atmopsheric models are not as accurate (to the level of the requirements), but also the manufacturing process is not as easy as it looks like. A manufacturing error of 1 arcsecond on the prisms angles will create a variation of 12\% on the residuals.\\
 In this paper, in order to fully understand the effects of ADC residuals, we will simulate its effect on synthetic spectra in section \ref{sec:simulations}, show the results and the analysis in sections \ref{sec:results} \& \ref{sec:discussion}.

 \section{Simulations}
 \label{sec:simulations}
 To test the effect of atmospheric dispersion on astronomical observations, and specially the impact on RV, we generated synthetic spectra affected by dispersion. These were created with a tool developed for the work of \cite{Martins2018} which create spectra that mimic reduced HARPS observations. HARPS, \citep[High-Accuracy Radial velocity Planet Searcher,][]{Mayor2003}, is one of the most precise spectrographs on Earth, with a design that was reproduced on HARPS-N and that greatly inspired ESPRESSO. It is considered the prototypical fiber-fed spectrograph. Since we had access to the HARPS pipeline, to expedite and simplify the reduction process with no loss of generality, we decided to use the HARPS spectrograph as our spectral template. It is important to state that no ADC is assumed in all the simulations during this work except for section \ref{subsec:adcresiduals}.

\subsection{Star and fiber simulations}
\label{subsec:star_simulations}
In order to simulate the point spread function (PSF) of the scientific object on the guiding camera, we start with a 2D Gaussian function: \\
$ I(x,y) = e^{-\frac{(x-x_0)^2}{2\sigma_x^2}-\frac{(y-y_0)^2}{2\sigma_y^2}}$ \\
where I is the intensity, (x,y) represent the displacement with respect to the origin ($x_0$,$y_0$);
$\sigma_x$ and $\sigma_y$ represent the sigma of the distribution in arcseconds. Since we are interested in studying the effect of atmospheric dispersion, it is important to understand how it will affect the position of the barycenter of the PSF. The atmospheric dispersion is a function of wavelength which will shift the barycenter as function of wavelength as well. For a star at zenithal angle Z, the angular atmospheric dispersion is proportional to Z (see Eq. \ref{eq:dispersion}). At the focus of a telescope, this will produce a linear dispersion with a well defined direction \citep{Wynne1986}.
For simplification purposes but with no loss of generality, we will assume that this displacement is along the $x$-axis of the guiding camera, while the displacement along the $y$-axis will remain constant. We also assumed that the seeing is constant as a function of wavelength.
As for $\sigma_x$ and $\sigma_y$ they are functions of the seeing and will be considered equal to 0.8", a typical value of La Silla observatory where HARPS is installed.
We will use the dispersion model presented by \cite{Filippenko1982}, hereafter the Filippenko model, to estimate the dispersion for the wavelength range of interest between 380 nm and 700 nm, the wavelength coverage of HARPS,
(see Figure \ref{Fig:filippenko}). The atmospheric dispersion computed is always relative to a reference wavelength. In this case, we assumed that the reference wavelength is 380 nm, centered at the center of the guiding camera. The barycenter of the PSF corresponding to the reference wavelength is considered the origin of displacement ($x_0$,$y_0$). The fiber size is assumed to be 1", a value equal to the HARPS fiber aperture.
\begin{figure}
\centering
\includegraphics[width=\hsize]{filippenko.png}
\caption{Differential atmospheric dispersion for a zenith angle of 60$^{\circ}$ using the Filippenko model for La Silla observatory.}
\label{Fig:filippenko}
\end{figure}
\\
In addition, we take into consideration in the simulations the quantum efficiency (QE) of the guiding camera by multiplying the intensity of each wavelength by the corresponding QE. The variation of QE with wavelength will affect the position of the barycenter of the PSF.
For that, we assumed two different cameras, one more efficient in the blue and the other being more efficient in the red, the Bigeye G-132 and G-283, respectively.
Figure \ref{Fig:qe} shows the variation of the QE of both cameras of interest where it is clear that the Bigeye G-132 has higher efficiency in the blue part of the spectra.

\begin{figure}
\centering
\includegraphics[width=\hsize]{qe_both.png}
\caption{Absolute quantum efficiency of Bigeye G-132 and G-283 scaled to wavelength range of interest.}
\label{Fig:qe}
\end{figure}

The atmospheric dispersion is also a function of the zenithal angle of observation (Z; see Eq. \ref{eq:dispersion}). Therefore we simulated a PSF, taking into consideration the amount
of dispersion (0 to 2 arcseconds) with Z taking the values between 0$^{\circ}$ and 60$^{\circ}$ presented in Table \ref{Table:disp}.

\begin{table}
\caption{Zenithal angles (and airmass) of interest and their corresponding levels of atmospheric dispersion for La Silla observatory.}
\label{Table:disp}
\centering
\begin{tabular}{c c c}
\hline
Z ($^{\circ}$) & Airmass & Dispersion (mas) \\
\hline
0 & 1 &0 \\
5 & 1.004 & 115 \\
10 & 1.015 & 230 \\
15 & 1.035 & 350 \\
30 & 1.154 & 750 \\
45 & 1.414 & 1300 \\
60 & 2 & 2250 \\
\hline
\end{tabular}
\end{table}

The case of no dispersion (Z = 0$^{\circ}$) will be used as reference in the remaining of this study. In Figure \ref{Fig:dispersion}, we show the effect of atmospheric dispersion on the shape of the PSF for 2 cases:
no dispersion (Z = 0$^{\circ}$), and high dispersion (Z = 60$^{\circ}$). The amount of dispersion will shift the position of the barycenter of the PSF; while the variation of QE will affect the intensity levels of each image. In the case of no dispersion, all the spots of different wavelengths will fall on the same location on the guiding camera and hence explaining the circular shape of the final spot. As for the case of high dispersion, each spot will be displaced from the reference spot (the spot corresponding to the reference wavelength) by an amount equal to the corresponding atmospheric dispersion which explain the elliptical shape of the final spot.
\begin{figure}
\centering
\includegraphics[width=\hsize]{dispersion.png}
\caption{Effect of atmospheric dispersion on the simulated spots taking into consideration the QE of the camera. Left: the case where there is no dispersion (Z = 0$^{\circ}$); right: the case where the dispersion is the highest (Z = 60$^{\circ}$). The dispersion is shown relative to the extreme blue wavelength. The dashed lines show the position of the origin (x$_0$,y$_0$).
}
\label{Fig:dispersion}
\end{figure}
\\
In order to compute the fraction of light at each wavelength that will enter the fiber, the fiber need to be centered at the barycenter of each spot for all the cases of interest, which will maximize the amount of light entering the fiber.
The barycenter is computed as follows:
\begin{equation}
\bar{x} = \dfrac{\sum\limits_{i=0}^{n}x_{i} \times m_{i}}{\sum\limits_{i=0}^{n}m_{i}} \\
\bar{y} = \dfrac{\sum\limits_{i=0}^{n}y_{i} \times m_{i}}{\sum\limits_{i=0}^{n}m_{i}} ,
\label{eq:barycenter}
\end{equation}
where $m_i$ is the intensity of each point of coordinates $x_i,y_i$; and $n$ corresponds to the number of images used (equivalent to the wavelength step between 380 and 780 nm; 84 images in our case). Due to atmospheric dispersion, the elongation of the PSF will increase compared to the no dispersion case.
This will vary the fraction of each wavelength entering the fiber for different values of dispersion.
From Figure \ref{Fig:fiber_eff}, we notice that with increasing Z, less light will enter the fiber and specially at the two limits of the wavelength range since the guiding cameras are usually wavelength centered (around 550 nm). We also notice that the position of the peak is not at the same wavelength with increasing Z. This variation indicates that the position of the barycenter is not constant. We can also notice that the amount of light entering, in the blue, is higher for the blue camera (Bigeye G-132), which is expected. In fact the guiding camera has a strong impact on which wavelengths are lost, and on which proportion.
\begin{figure}
\centering
\includegraphics[width=\hsize]{fiber_efficiency_espresso.png}
\includegraphics[width=\hsize]{fiber_efficiency.png}
\caption{Percentage of light of each wavelength entering the fiber for all values of Z. Bigeye G-132 (top); Bigeye G-283 (bottom).}
\label{Fig:fiber_eff}
\end{figure}

\subsection{Synthetic spectra simulations}
\label{subsec:spectra_simulations}
In order to simulate our synthetic spectra, we followed the method described by \cite{Martins2018}. It is important to state that in order to take into consideration the atmospheric dispersion effect on the spectra, the original code was updated to our requirements. During the design of an ADC, one of the important factors is the wavelength range where the ADC will be installed. Therefore, we were not interested in this work, in testing the effect of different spectral types and we limited our tests to the solar spectrum \citep{Kurucz2006} of G spectral type. Here we summarize the main steps and explain how the atmospheric dispersion effect
was taken into consideration:
\begin{itemize}
\item We started from the very high resolution (500000) solar spectrum of \cite{Kurucz2006}.
\item The spectrum  was degraded to the intended resolution of HARPS using the algorithm presented in \cite{Neal2017}.
\item The continuum was scaled in flux to match the expected flux of the target as computed from HARPS exposure time calculator (Table \ref{Table:ETC}).
\begin{table}
\caption{HARPS ETC settings. The settings not in this table were left to their default values.}
\label{Table:ETC}
\centering
\begin{tabular}{l r}
\hline
\textbf{Target input flux} &  \\
\textbf{distribution} & \\
Template spectrum & G2V (Kurucz) \\
\textbf{Sky conditions} & \\
Airmass & Table \ref{Table:disp} \\
Seeing & 0.8" \\
\textbf{Instrument set-up} & \\
Exposure time & 300s, 100s, 50s \\
Signal to Noise & 220, 130, 90 \\
\hline
\end{tabular}
\end{table}
\item Each spectrum is then multiplied by the corresponding fiber efficiency ( Figure \ref{Fig:fiber_eff}) based on the assumed zenithal angle (amount of dispersion).
\item The noise is assumed to be Gaussian photon noise.
\end{itemize}
We simulated 10$^4$ spectra with different photon noise component for each dispersion case.
These spectra will be then used to test and evaluate the effect of atmospheric dispersion.

 \subsection{From synthetic spectra to radial velocities}
 \label{subsec:spectrarv}
 In order to reduce these spectra and compute the RV, HARPS pipeline (version 3.4) was used.
 We produced two sets of reduced data by turning on and off the flux correction function, which is a pipeline variable.
 The flux correction function will compare the flux from each spectra to a template based on the spectral type of the target.
 These templates are built from very high signal-to-noise ratio (SNR) observations \citep{espr2019} of standard stars at low airmass.
 Then, it will apply a correction in order to scale the flux to the same level of the template by multiplying it by the ratio of flux between the observed and the template.
 In that case, all our spectra will have the same flux distribution as the template, which ensures that variable atmospheric conditions will not
 induce any systematic effects in the cross-correlation function (CCF) computation. In case the correction was too small or too big due to mismatch with the flux template,
 the flux correction function is automatically switched off. In fact, according to the pipeline, if the ratio between the observed flux and the template is lower than 0.25 or higher than 3, the function will be off. We will not focus here on these limits as it is behind the scope of this paper. \\ In our case, the spectra were simulated assuming a G2 solar spectral type.
 Based on HARPS pipeline, no flux template exist for this spectral type. Hence, we built our own G2 flux template using the same technique described above. We simulated a very high SNR spectrum at an airmass of 1, where the dispersion is null. In this study, all the spectra were corrected using our template. As mentioned above, HARPS pipeline was used to compute the RV. The RV in HARPS pipeline is computed following the CCF technique where high-resolution spectra are cross-correlated with numerical masks over a range around the host's RV \citep{Baranne1996}. The CCF function from the HARPS pipeline is used with the same inputs for all the spectra tested of all the cases. In that case we guarantee that the variation of RV, if any, is only due to the atmospheric dispersion effect.

\section{Results}
\label{sec:results}
As mentioned before, several tests were done in order to understand and evaluate the effect of the atmospheric dispersion to be able
to justify the requirements on the residuals of an ADC.

As well known, and discussed, e.g. in \cite{Pepe2008} and \cite{Fischer2016}, atmospheric dispersion will introduce a slope variation in the spectral continuum.
To verify this effect, we divided spectra, simulated with atmospheric dispersion not corrected by an ADC, of different zenithal angles by the reference (Z = 0$^{\circ}$) to see if there is any slope variation.
This effect can be clearly seen in Figure \ref{Fig:slope}. We only show four cases for visualization purposes. It is clear that starting at Z = 15$^{\circ}$, the slope variation starts to be more severe until it reaches the worst case for Z = 60$^{\circ}$ where a big part of the blue light is lost due to atmospheric dispersion. This slope variation should be corrected with the flux correction function at the level of the CCF.
\begin{figure}
\centering
\includegraphics[width=\hsize]{slope_variation.png}
\caption{Continuum slope variation for different Z cases compared to the reference Z = 0$^{\circ}$. The slopes were extracted from the synthetic spectra affected by atmospheric dispersion not corrected by an ADC.}
\label{Fig:slope}
\end{figure}

\subsection{RV precision degradation}
In order to test the effect of atmospheric dispersion on RV precision, if any, we simulated several spectra for each zenithal angle in Table \ref{Table:ETC}. To be able to characterize the effect of both random noise and systematic effects, 10$^4$ spectra were generated for each case assuming no ADC to correct the atmospheric dispersion. From Figure \ref{Fig:added}, we notice that after 7000 spectra, the results converge to the average value, allowing us to obtain both random and systematic error information. We also tested that the spectra have a gaussian distribution, well characterized by its mean and standard deviation.

\begin{figure}
   \centering
   \includegraphics[width=\hsize]{added.png}
   \caption{Systematic RV variation as function of number of simulations. It is clear that after 7000 spectra the results converge to the average value plotted as red dashed line.}
   \label{Fig:added}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\hsize]{RV_variation_off.png}
\caption{Computed RV values of each set of simulated spectra when the flux correction function was \textbf{off}. The error bars correspond to the 1 $\sigma$ random error.}
\label{Fig:rv_off}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\hsize]{RV_variation_on.png}
\caption{Computed RV values of each set of simulated spectra when the flux correction function was \textbf{on}. The error bars correspond to the 1 $\sigma$ random error.}
\label{Fig:rv_on}
\end{figure}

As it was shown before (Figure \ref{Fig:slope}), the atmospheric dispersion will introduce a slope variation of the continuum. This variation will hence give different weights
to the spectral lines that will be used in the computation of the CCF. This will introduce a variation in the final RV. In figure \ref{Fig:rv_off}, we plot
this variation in RV as function of zenithal angles. It is clear that this variation increases with Z, and hence with the amount of increasing dispersion.
By definition, the flux correction function is supposed to correct for these atmospheric conditions variation which is also clear in Figure \ref{Fig:rv_on}. For a zenithal angle above Z = 40$^{\circ}$, the flux correction function is automatically switched off in the pipeline, explaining the sudden increase in the systematic error. The error bars shown are the random errors due to photon noise. \\
Since we are interested only in the systematic variation of the RV, which corresponds to the RV variation of all the zenithal cases compared to the one where no atmospheric dispersion was considered (zero-point RV at Z = 0$^{\circ}$), we decided to turn off the noise component in the simulated spectra in the rest of the paper, as we believe that the photon noise might affect our analysis. In this case, we will be only testing the effect of atmospheric dispersion.

\subsection{Flux losses due to different guiding cameras}
Due to the fact that atmospheric dispersion will increase the elongation of the PSF of the target on the focal plane where light is injected in the fiber (and imaged by the guiding camera), flux losses might be introduced. In fact, the atmospheric dispersion, if not corrected by an ADC, can be higher than 2", a value that is twice the typical size of the fibers used to inject the light into a spectrograph. Another problem, is the issue of guiding. Different guiding cameras, can lead to different fractions of each wavelength entering the fiber due to different QE and guiding algorithms. In this section, we will only focus on the effect due to different QE. We tested two different cameras in order to understand the effect of atmospheric dispersion on the variation of the barycenter position, the amount of flux that will enter the fiber, and the effect on the RV, if any.
\\
In order to test the effect of atmospheric dispersion on the position of the barycenter, we simulated the spot of the target as described in section 2.
These simulations were done for 2 different CCDs: Bigeye G-132 and Bigeye G-283, with high QE in the blue part of the spectra and the red one, respectively. \\

We were able to measure the position of the barycenter on each case of Z using Eq. \ref{eq:barycenter}.
These positions were plotted as function of atmospheric dispersion (Figure \ref{Fig:bary_disp}). As expected, the position of the barycenter is proportional to the amount of dispersion: the higher it is, the further the barycenter is from the origin.

\begin{figure}
\centering
\includegraphics[width=\hsize]{barycenter_shift_mas.png}
\caption{Variation of the position of the barycenter with the variation of zenithal angle, for the two cameras of interest. These positions were measured relative to the position of no atmospheric dispersion (Z = 0$^{\circ}$) as reference.}
\label{Fig:bary_disp}
\end{figure}

We also computed the maximum flux loss at 380 nm for the two cameras tested. It is clear from Figure \ref{Fig:flux_loss_max} that the choice of the guiding camera can increase or decrease the flux losses. In Figure \ref{Fig:flux_loss_max}, for Z = 60$^{\circ}$, the 2 cameras of interest suffers from total loss of flux at 380 nm due to the huge amount of atmospheric dispersion.

\begin{figure}
\centering
\includegraphics[width=\hsize]{flux_losses_380.png}
\caption{Maximum flux loss at 380 nm for the two cameras of interest.}
\label{Fig:flux_loss_max}
\end{figure}

\subsection{Flux losses due to different melt data}
\label{subsec:melt}
The melt data, glass characteristics of each batch produced, can cause a serious problem when designing an ADC \citep{Wehbe2019}. In fact, different melt data result with different Sellmeier coefficients
than the ones in the catalog normally used during the design phase of an ADC. \cite{Wehbe2019} provides a detailed analysis of the different cases that we will not
mention here. We will use one of the cases, to test the effect, on the flux losses. In fact, different cases of melt data will result different shapes of
ADC residuals and different PV that might alter the flux losses specially in the blue wavelengths. The variation between the melt cases will result in a variation in the Sellmeier coefficients that will affect the computed index of refraction. The latter will introduce a variation in the dispersion of the medium and hence the amount of residuals expected from the ADC. In Figure \ref{Fig:melt_data_residuals}, we show the different shapes and amounts of residuals
of different melt data (Table \ref{Table:melt}) for an ADC combination of K7 + S-FPL51Y \citep[the glass combination used for the ADC of ESPRESSO,][]{Cabral2012}. It is also important to state that the only variable between these curves was the melt data. The
prism angles and the reference wavelength were maintained as computed during the design phase using the catalog Sellmeier coefficients. It is clear that not only the shape and the amount is changing in Figure \ref{Fig:melt_data_residuals}, but also the value of the residuals for the wavelength of 380 nm. In the catalog case, the ADC was designed to have a minimum of residuals at 380 nm. Different melt data, return different residuals at 380 nm. This variation is severe when it comes to flux losses.
\begin{figure}
   \centering
   \includegraphics[width=\hsize]{melt_data_residuals.png}
   \caption{Residuals for all the melt data cases of  K7 + S-FPL51Y.}
   \label{Fig:melt_data_residuals}
\end{figure}

In order to test the effect of melt data on flux losses, each residuals curve from Figure \ref{Fig:melt_data_residuals} was used to create a fiber efficiency injection curve as described in section \ref{subsec:star_simulations}.
These new fiber curves were then used to compute the flux floss at 380 nm for each case. The results are shown in Figure \ref{Fig:melt_data_flux_loss}. It is clear that different melt data cases, with different PV, similar to sky dispersion values, return different flux losses. This variation between the sky dispersion and the ADC residuals is due to the variation of the shape of the dispersion between on-sky and ADC.

\begin{table}
   \caption{Melt data designation}
   \label{Table:melt}
   \centering
   \begin{tabular}{l r}
   \hline
   A & K7 \\
   B & S-FPL51Y \\
   No number & Catalogue data \\
   With number & Melt data  \\
   \hline
   \end{tabular}
   \end{table}

\begin{figure}
   \centering
   \includegraphics[width=\hsize]{flux_losses_sky_adc_380_annotate.png}
   \caption{Maximum flux loss at 380 nm for all the melt cases.}
   \label{Fig:melt_data_flux_loss}
\end{figure}

\section{Discussion}
\label{sec:discussion}
A key role when designing an ADC, is the amount of residuals set by the requirements (50 mas in the case of ESPRESSO). The tougher the requirements are, the harder it is to design, as there are less options in terms of glass pairs and the higher will be the cost. Building an ADC capable of returning the required residuals is not only a manufacturing challenge, but also a design challenge from chosing the correct glasses, to the correct angles, to the choice of the reference wavelength where the ADC will be optimized and the choice of the guiding camera and its QE. In fact, QE effect can be corrected with the use of a spectral filter in the guiding camera, nevertheless this is an option that should be analysed carefully as it would reduce the number of photons on the guiding camera and could jeopardize the efficiency of this task that has a huge impact on the total observation flux efficiency. All these factors, play a role when setting the residuals amount that, to our knowledge, is not a straightforward and well investigated area. In order to be able to set the requirements, we need to start by evaluating the effects of atmospheric dispersion on astronomical observations to be able to define the residuals of an ADC. In this section we will do so.

\subsection{Impact of atmospheric dispersion}
Without an ADC, the atmospheric dispersion can be easily modelled using the Filippenko's model taking into consideration the different atmospheric parameters (temperature, pressure, RH and Z). By varying the zenithal angle, one can vary the amount of atmospheric dispersion as the shape will remain constant and similar to Figure \ref{Fig:filippenko}. It is clear from Figure \ref{Fig:rv_off} that the systematic error on RV start to be larger than 1 m/s for Z = 5$^{\circ}$ when the flux correction function is off. \\
When the flux correction function is on, these systematic errors are reduced to below 0.2 cm/s for Z = 5$^{\circ}$ (see Figure \ref{Fig:systematic_nonoise}) which implies the need of such a function in the pipelines. It also implies that an ADC is certainly needed to reduce the systematic errors to the level of 0.1 m/s. \\
To go further into investigation, we simulated spectra using the same technique described in \ref{subsec:spectra_simulations}, but this time by setting the noise to zero. In that case we will be able to fully characterize the systematic errors due to atmopsheric dispersion. It is clear from Figure \ref{Fig:systematic_nonoise}, that the effects of atmospheric dispersion up to Z = 5$^{\circ}$ (equivalent to $\sim$ 100 mas), on RV, are almost negligible (for a 0.1 m/s objective). With increasing zenithal angle, hence increasing the amount of dispersion, the amount of blue light of the spectra that enters the fiber decreases (see Figure \ref{Fig:fiber_eff}). The loss of light, lines in the spectra, will increase the error on RV as seen in Figure \ref{Fig:systematic_nonoise} with increasing Z.

\begin{figure}
   \centering
   \includegraphics[width=\hsize]{systematic_error_nonoise.png}
   \caption{Systematic error variation as function of zenithal angle. The systematic error is behaving as expected, increasing with increasing Z.}
   \label{Fig:systematic_nonoise}
\end{figure}

\subsection{Impact of ADC residuals}
\label{subsec:adcresiduals}
Since an ADC is crutial in order to reduce the effects of atmospheric dispersion on astronomical observations, it is important to understant the effect of the ADC residuals. This effect, in general, should be similar to the effect of atmospheric dispersion by decreasing the zenithal angle to 5$^{\circ}$ or below which is corresponding to a dispersion of 100 mas and below, a value close to the current ADC residuals. In fact, the shape of the ADC residuals is not similar to the shape of the atmospheric dispersion, this is why we decided to use the ADC residuals curve (we simulated the case of ESPRESSO's ADC \citep{Cabral2012}, Figure \ref{Fig:melt_data_residuals}) instead of the Filippenko's curve (Figure \ref{Fig:filippenko}) and do a similar test as it was done previously in this work. It is also important to state the fact that ADC residuals are also proportional to Z. In this subsection, we will consider only the residuals corresponding to Z = 60$^{\circ}$, which is the worst case scenario. The residuals curves shown in Figure \ref{Fig:melt_data_residuals}, are used to produce curves similar to Figure \ref{Fig:fiber_eff}, which are used then to simulate spectra assuming the existence of an ADC. These new spectra were used to compute the RV as described in section \ref{subsec:spectrarv}. In this subsection, we will evaluate the effect of atmospheric dispersion residuals in mas instead of zenithal angles. \\
As discussed in subsection \ref{subsec:melt}, different melt data will return different amounts and shapes of ADC residuals. We will compare the different cases of melt, with the atmospheric dispersion results of the same amounts.

\begin{figure}
   \centering
   \includegraphics[width=\hsize]{systematic_error_nonoise_sky_adc.png}
   \caption{Similar to figure \ref{Fig:systematic_nonoise}, but but adding the different melt data residuals results. The points from Figure \ref{Fig:systematic_nonoise} represent the case where an ADC is used.}
   \label{Fig:systematic_nonoise_melt}
\end{figure}

The results from Figure \ref{Fig:systematic_nonoise}, also represent the case where an ADC is used. Each zenithal angle in that case will represent an amount of dispersion residuals as indicated in Table \ref{Table:disp}. It is clear from Figures \ref{Fig:melt_data_flux_loss} and \ref{Fig:systematic_nonoise_melt}, that ADC residuals up to 100 mas will create an error of $\sim$ 1 cm/s (mainly due to melt data since they produce different residuals shapes) and $\sim$ 1\% on the RV precision and flux losses, respectively. It is also clear from the same figures, that different melt data can return different residuals, different RV variation and different flux losses. It is important, when designing an ADC, to give more attention to the melt data as they can cause a serious problem specially when it comes to flux losses that seems more critical than RV variation.

\section{Conclusions}
As we showed above, atmospheric dispersion will largely contribute to the errors if not properly corrected by an ADC. The residuals of an ADC will affect the RV precision as well as the flux losses when designing a high-resolution spectrograph. Focusing on the spectral range between 380 nm and 780 nm, the case used in this analysis, ADC residuals of the order of 100 mas are guaranteed not to introduce large errors that could impact a precision of the level of 10 cm/s. In fact, 100 mas PV, will have an impact of $\sim$ 1 cm/s in terms of RV precision. This same 100 mas, will have an impact of $\sim$ 1\% in flux losses. The latter could be the most important factor of the requirements on the ADC specifications. These results are summarized in Figures \ref{Fig:systematic_nonoise_melt}, and \ref{Fig:conclusion}.

\begin{figure}
   \centering
   \includegraphics[width=\hsize]{conclusion_plot.png}
   \caption{Conclusion plot summarizing the requirements on ADC residuals. Top: expected RV error due to ADC residuals (similar to Figure \ref{Fig:systematic_nonoise}); bottom: expected flux losses at 380 nm due to ADC residuals (similar to Figure \ref{Fig:flux_loss_max}).}
   \label{Fig:conclusion}
\end{figure}

We also noticed that the choice of the guiding camera can play a role and affect specially the flux losses. We presented a study that will allow to define the specifications of an ADC based on the requirements on RV precision and flux losses, taking into account the guiding system, responsible of the fiber-injection, and on the optical optimization based on melt data. This procedure will allow to define the specifications of an ADC in a more accurate way avoiding the (typical) conservative approach that increases (unnecessarily) the complexity and the cost of an ADC.

\section*{Acknowledgements}
The first author is supported by an FCT fellowship (PD/BD/135225/2017), under the FCT PD Program PhD::SPACE (PD/00040/2012). Jorge H. C. Martins is supported in the form of a work contract funded by national funds through Funda\c{c}\~{a}o para a Ci\^{e}ncia e Tecnologia (DL 57/2016/CP1364/CT0007). This work was supported by FCT/MCTES through national funds and by FEDER - Fundo Europeu de Desenvolvimento Regional through COMPETE2020 - Programa Operacional Competitividade e Internacionaliza\c{c}\~{a}o by these grants: UID/FIS/04434/2019; PTDC/FIS-AST/32113/2017 \& POCI-01-0145-FEDER-032113.

\bibliographystyle{mnras}


\section{\textbf{Introduction}}

The Dirac operator is the relativistic Schr\"{o}dinger operator in quantum
physics.The basic and comprehensive results about Dirac operators were given
in \cite{Lev}. Inverse problems for the Dirac operators have been
extensively well studied in various publications (see \cite{Alb}, \cite{Gsy1
, \cite{Gus}, \cite{Hor}, \cite{Ksk2}, \cite{Ksk} , \cite{wang2} and the
references therein).The subject of fractional calculus has acquired
significant popularity and major attention from several authors in various
science due mainly to its direct involvement in the problems of differential
equations in mathematics, physics, engineering and others. This topic is
initiated by (\cite{Abdel}, \cite{halil}). The fractional calculus has
gained an interesting area in mathematical research and it has attracted the
attention of several scholars. A\ variety of new results can be found (\cit
{abu}, \cite{abu1} and references therein). In recent years, scholars have
focussed on a fractional generalization of the well known Sturm-Liouville
and Dirac problems (\cite{al}, \cite{al1}, \cite{tuna}, \cite{tuna2}, \cit
{tuna3}, \cite{gul}, \cite{kh}, \cite{kl} and \cite{rv}).

Inverse nodal problem was started for the Sturm--Liouville operator by
McLaughlin \cite{mc1} in 1988. In 1989, Hald and McLaughlin proved that it
is sufficient to know the nodal points to determine the potential function
of the regular Sturm--Liouville problem\ with more general boundary
conditions uniquely and gave some numerical schemes for the reconstruction
of the potential from nodal points \cite{H}. Yang proposed an algorithm to
solve an inverse nodal problem for the Sturm--Liouville operator in 1997
\cite{yang}. Such problems have been considered by several researchers in
\cite{Br2}, \cite{ch},, \cite{guo1}, \cite{Law}, \cite{Ozkan}, \cite{wang},
\cite{wang3}, \cite{wei}, \cite{Yur}, \cite{yng1}, \cite{Yang3} and \cit
{Yang4} ) and other works.The inverse nodal problems for the Dirac operators
with various boundary conditions have been studied and shown that a dense
subset of the zeros of the first component of the eigenfunctions alone can
determine the coefficients of discussed problem \cite{Guo}, \cite{Yang2} and
\cite{Yang5}. In \cite{mor}, the authors have developed the spectral theory
for a conformable fractional Sturm-Liouville problem and have proved
uniqueness theorem with respect to the nodal points.

Nowadays, the studies concerning the perturbation of a differential operator
by a Volterra type integral operator, namely the integro-differential
operator has acquired significant popularity and major attention from
several authors and \ take significant place in the literature.\cite{But},
\cite{But 2}, \cite{G}, \cite{Kur}, \cite{Kur2} and \cite{B} ).
Integro-differential operators are nonlocal, and therefore they are more
difficult for investigation, than local ones. New methods for solution of
these problems are being developed. For Sturm-Liouville type
integro-differential operators, there exist some studies about inverse
problems but there is very little study for Dirac type integro-differential
operators. The inverse nodal problem for Dirac type integro-differential
operators was first studied by \cite{kesk3}. In their study, it is shown
that the coefficients of the differential part of the operator can be
determined by using nodal points and nodal points also gives the partial
information about integral part. In \cite{kesk4}, the authors considered
Dirac type integro-differential operators with boundary conditions depend on
the spectral parameter linearly.

\section{\textbf{Conformable\ Fractional Preliminaries}}

Firstly, we want to recall some basic definitions and properties of
conformable fractional calculus which can be found in [\cite{halil}, \cit
{Abdel}].

\begin{definition}
Let $f:[0,\infty )\rightarrow
\mathbb{R}
$ be a given function. Then the conformable fractional derivative of $f$ of
order $\alpha $ is defined by:

$D_{x}^{\alpha }f(x)=\underset{\in \rightarrow 0}{\lim }\dfrac{f(x+\epsilon
x^{1-\alpha })-f(t)}{\epsilon }\alpha ,$ $D_{x}^{\alpha }f(0)=\underset
t\rightarrow 0^{+}}{\lim }D^{\alpha }f(x),$\newline
for all $x>0,$ $\alpha \in (0,1].$ If this limit exist and finite at $x_{0},$
we say $f$ is $\alpha -$differentiable at $x_{0}.$ Note that if $f$ is
differentiable, then $D^{\alpha }f(x)=x^{1-\alpha }f^{\prime }(x).$
\end{definition}

\begin{definition}
The conformable fractional Integral starting from $0$ of order $\alpha $ is
defined by

$I_{\alpha }f(x)=\int_{0}^{x}f(t)d_{\alpha }t=\int_{0}^{x}t^{\alpha
-1}f(t)dt,$\newline
for all $x>0.$
\end{definition}

\begin{lemma}
Let $f:[a,\infty )\rightarrow
\mathbb{R}
$ be any continuous function. Then, for all $x>a$, we have $D_{x}^{\alpha
}I_{\alpha }f(x)=f(x).$
\end{lemma}

\begin{lemma}
Let $f:(a,b)\rightarrow
\mathbb{R}
$ be any differentiable function. Then, for all $x>a$, we have $I_{\alpha
}D_{x}^{\alpha }f(x)=f(x)-f(a).$
\end{lemma}

\begin{theorem}
\textbf{(}$\alpha -$\textbf{integration by parts): }Let $f,g:[a,b
\rightarrow
\mathbb{R}
$ be two conformable fractional differentiable functions. Then,

$\int_{a}^{b}f(x)D_{x}^{\alpha }g(x)d_{\alpha
}x=f(b)g(b)-f(a)g(a)-\int_{a}^{b}g(x)D_{x}^{\alpha }f(x)d_{\alpha }x$
\end{theorem}

\begin{theorem}
($\alpha -$\textbf{Leibnitz Rule}) Let $t^{\alpha -1}f(x,t)$ and $t^{\alpha
-1}f_{x}(x,t)$ be continuous in $x$ on some regions of the $(x,t)-$plane,
including $a(x)\leq t\leq b(x),$ $x_{0}\leq x\leq x_{1}.$ If $a(x)$ and
b(x) $ are both $\alpha -$differentiable for $x_{0}\leq x\leq x_{1},$ then

$D_{x}^{\alpha }\left[ \int_{a(x)}^{b(x)}f(x,t)d_{\alpha }t\right]
=f(x,b(x))b^{\alpha -1}(x)D_{x}^{\alpha }b(x)-f(x,a(x))a^{\alpha
-1}(x)D_{x}^{\alpha }a(x)+\int_{a(x)}^{b(x)}D_{x}^{\alpha }f(x,t)d_{\alpha
}t.$
\end{theorem}

\begin{definition}
The space $C_{\alpha }^{n}[a,b]$ consists of all functions defined on the
interval $[a,b]$ which are continuously $\alpha -$differentiable up to order
$n.$
\end{definition}

\section{\textbf{Conformable Fractional Dirac Systems}}

In this work, we consider the following one-dimentional conformable
fractional Dirac type integro-differential syste
\begin{equation}
BY+\Omega (x)Y+\int\limits_{0}^{x}M(x,t)Yd_{\alpha }t=\lambda Y,\text{ \
x\in (0,\pi ),
\end{equation
with the boundary condition
\begin{eqnarray}
\sin \theta y_{1}(0)+\cos \theta y_{2}(0) &=&0\medskip  \\
\sin \beta y_{1}(\pi )+\cos \beta y_{2}(\pi ) &=&0
\end{eqnarray
where, $\theta $ is a real constant, $\lambda $ is the spectral parameter,
\newline
$B=\left(
\begin{array}{cc}
0 & D_{x}^{\alpha } \\
-D_{x}^{\alpha } &
\end{array
\right) ,$ \ $\Omega (x)=\left(
\begin{array}{cc}
p(x) & 0 \\
0 & r(x
\end{array
\right) ,$ $M(x,t)=\left(
\begin{array}{cc}
M_{11}(x,t) & M_{12}(x,t) \\
M_{21}(x,t) & M_{22}(x,t
\end{array
\right) ,$ \ $Y=\left(
\begin{array}{c}
y_{1} \\
y_{2
\end{array
\right) ,$ $p(x)$, $r(x)$, $M\left( x,t\right) $ are real-valued conformable
fractional differentiable functions and $x^{\alpha -1}p(x)$ and $x^{\alpha
-1}r(x)$ are continuous on $(0,\pi )$.

Let $\varphi (x,\lambda )=\left( \varphi _{1}(x,\lambda ),\varphi
_{2}(x,\lambda )\right) ^{T}$ be the solution of (1) satisfying the initial
condition $\varphi (0,\lambda )=(\cos \theta ,-\sin \theta )^{T}$. It is
clear that $\varphi (x,\lambda )$ is an entire function of $\lambda $
satisfies the following conformable fractional Volterra integral equations
\begin{equation}
\begin{array}{l}
\varphi _{1}(x,\lambda )=\cos \theta \cos \left( \lambda \dfrac{x^{\alpha }}
\alpha }\right) +\sin \theta \sin \left( \lambda \dfrac{x^{\alpha }}{\alpha
\right) \medskip \\
+\int_{0}^{x}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha
\right) p(t)\varphi _{1}(t,\lambda )d_{\alpha }t+\int_{0}^{x}\cos \left(
\lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha }\right) r(t)\varphi
_{2}(t,\lambda )d_{\alpha }t\medskip \\
+\int_{0}^{x}\int_{0}^{t}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}
\alpha }\right) \left\{ M_{11}(t,\xi )\varphi _{1}(\lambda ,\xi
)+M_{12}(t,\xi )\varphi _{2}(\lambda ,\xi )\right\} d_{\alpha }\xi d_{\alpha
}t\medskip \\
+\int_{0}^{x}\int_{0}^{t}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}
\alpha }\right) \left\{ M_{21}(t,\xi )\varphi _{1}(\lambda ,\xi
)+M_{22}(t,\xi )\varphi _{2}(\lambda ,\xi )\right\} d_{\alpha }\xi d_{\alpha
}
\end{array
\end{equation
\begin{equation}
\begin{array}{l}
\varphi _{2}(x,\lambda )=\cos \theta \sin \left( \lambda \dfrac{x^{\alpha }}
\alpha }\right) -\sin \theta \cos \left( \lambda \dfrac{x^{\alpha }}{\alpha
\right) \medskip \\
-\int_{0}^{x}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha
\right) p(t)\varphi _{1}(t,\lambda )d_{\alpha }t+\int_{0}^{x}\sin \left(
\lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha }\right) r(t)\varphi
_{2}(t,\lambda )d_{\alpha }t\medskip \\
-\int_{0}^{x}\int_{0}^{t}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}
\alpha }\right) \left\{ M_{11}(t,\xi )\varphi _{1}(\lambda ,\xi
)+M_{12}(t,\xi )\varphi _{2}(\lambda ,\xi )\right\} d_{\alpha }\xi d_{\alpha
}t\medskip \\
+\int_{0}^{x}\int_{0}^{t}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}
\alpha }\right) \left\{ M_{21}(t,\xi )\varphi _{1}(\lambda ,\xi
)+M_{22}(t,\xi )\varphi _{2}(\lambda ,\xi )\right\} d_{\alpha }\xi d_{\alpha
}
\end{array
\end{equation}

\begin{lemma}
Let $\ f(x)$ be a function in $C_{\alpha }^{1}[0,\pi ],$ then

$\underset{\left\vert \lambda \right\vert \rightarrow \infty }{\lim }\exp
(-\left\vert \func{Im}\lambda \dfrac{\pi ^{\alpha }}{\alpha }\right\vert
)\int_{0}^{\pi }f(x)\cos \lambda \dfrac{x^{\alpha }}{\alpha }d_{\alpha }x=0$

and

$\underset{\left\vert \lambda \right\vert \rightarrow \infty }{\lim }\exp
(-\left\vert \func{Im}\lambda \dfrac{\pi ^{\alpha }}{\alpha }\right\vert
)\int_{0}^{\pi }f(x)\sin \lambda \dfrac{x^{\alpha }}{\alpha }d_{\alpha }x=0
\bigskip

\begin{theorem}
For $\left\vert \lambda \right\vert \rightarrow \infty ,$ the following
asymptotic formulae are valid
\begin{eqnarray}
&&\left. \varphi _{1}(x,\lambda )=\cos \left[ \lambda \dfrac{x^{\alpha }}
\alpha }-\mu (x)-\theta \right] +\frac{1}{2\lambda }\upsilon (x)\cos \left[
\lambda \dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \right. \medskip
\notag \\
&&-\frac{1}{2\lambda }\upsilon (0)\cos \left[ \lambda \dfrac{x^{\alpha }}
\alpha }-\mu (x)+\theta \right] +\frac{1}{2\lambda }\sin \left[ \lambda
\dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \int_{0}^{x}\upsilon
^{2}(t)d_{\alpha }t\medskip \\
&&-\frac{1}{2\lambda }K\left( x\right) \cos \left[ \lambda \dfrac{x^{\alpha
}{\alpha }-\mu (x)-\theta \right] -\frac{1}{2\lambda }L\left( x\right) \sin
\left[ \lambda \dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \medskip
\notag \\
&&+o\left( \frac{1}{\lambda }\exp (\left\vert \tau \right\vert \dfrac
x^{\alpha }}{\alpha })\right) ,  \notag
\end{eqnarray
\begin{eqnarray}
&&\left. \varphi _{2}(x,\lambda )=\sin \left[ \lambda \dfrac{x^{\alpha }}
\alpha }-\mu (x)-\theta \right] -\frac{1}{2\lambda }\upsilon (x)\sin \left[
\lambda \dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \right. \medskip
\notag \\
&&-\frac{1}{2\lambda }\upsilon (0)\sin \left[ \lambda \dfrac{x^{\alpha }}
\alpha }-\mu (x)+\theta \right] -\frac{1}{2\lambda }\cos \left[ \lambda
\dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \int_{0}^{x}\upsilon
^{2}(t)d_{\alpha }t\medskip \\
&&-\frac{1}{2\lambda }K\left( x\right) \sin \left[ \lambda \dfrac{x^{\alpha
}{\alpha }-\mu (x)-\theta \right] +\frac{1}{2\lambda }L\left( x\right) \cos
\left[ \lambda \dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right] \medskip
\notag \\
&&+o\left( \frac{1}{\lambda }\exp (\left\vert \tau \right\vert \dfrac
x^{\alpha }}{\alpha })\right) ,  \notag
\end{eqnarray
uniformly in $x\in \lbrack 0,\pi ],$ where, $\mu (x)=\dfrac{1}{2
\int_{0}^{x}(p(t)+r(t))d_{\alpha }t,$ $\upsilon (x)=\dfrac{1}{2}\left(
p(x)-r(x)\right) ,$ $K(x)=\int_{0}^{x}(M_{11}(t,t)-M_{22}(t,t))d_{\alpha }t,$
$L(x)=\int_{0}^{x}(M_{12}(t,t)-M_{21}(t,t))d_{\alpha }t$ \newline
and $\tau =\func{Im}\lambda .$
\end{theorem}
\end{lemma}

\begin{proof}
we denote
\begin{eqnarray*}
&&\left. \varphi _{1,0}(x,\lambda )=\cos \left( \lambda \dfrac{x^{\alpha }}
\alpha }-\theta \right) ,\medskip \right. \\
&&\left. \varphi _{1,n+1}(x,\lambda )=\int_{0}^{x}\sin \left( \lambda \dfrac
x^{\alpha }-t^{\alpha }}{\alpha }\right) p(t)\varphi _{1,n}(t,\lambda
)d_{\alpha }t\right. \\
&&+\int_{0}^{x}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha
\right) r(t)\varphi _{2,n}(t,\lambda )d_{\alpha }t\medskip \\
&&+\int_{0}^{x}\int_{0}^{t}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha
}{\alpha }\right) \left\{ M_{11}(t,\xi )\varphi _{1,n}(\lambda ,\xi
)+M_{12}(t,\xi )\varphi _{2,n}(\lambda ,\xi )\right\} d_{\alpha }\xi
d_{\alpha }t\medskip \\
&&+\int_{0}^{x}\int_{0}^{t}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha
}{\alpha }\right) \left\{ M_{21}(t,\xi )\varphi _{1,n}(\lambda ,\xi
)+M_{22}(t,\xi )\varphi _{2,n}(\lambda ,\xi )\right\} d_{\alpha }\xi
d_{\alpha }t,\medskip
\end{eqnarray*
\begin{eqnarray*}
&&\left. \varphi _{2,0}(x,\lambda )=\sin \left( \lambda \dfrac{x^{\alpha }}
\alpha }-\theta \right) ,\medskip \right. \\
&&\left. \varphi _{2,n+1}(x,\lambda )=-\int_{0}^{x}\cos \left( \lambda
\dfrac{x^{\alpha }-t^{\alpha }}{\alpha }\right) p(t)\varphi _{1,n}(t,\lambda
)d_{\alpha }t\right. \medskip \\
&&+\int_{0}^{x}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha }}{\alpha
\right) r(t)\varphi _{2,n}(t,\lambda )d_{\alpha }t\medskip \\
&&-\int_{0}^{x}\int_{0}^{t}\cos \left( \lambda \dfrac{x^{\alpha }-t^{\alpha
}{\alpha }\right) \left\{ M_{11}(t,\xi )\varphi _{1,n}(\lambda ,\xi
)+M_{12}(t,\xi )\varphi _{2,n}(\lambda ,\xi )\right\} d_{\alpha }\xi
d_{\alpha }t\medskip \\
&&+\int_{0}^{x}\int_{0}^{t}\sin \left( \lambda \dfrac{x^{\alpha }-t^{\alpha
}{\alpha }\right) \left\{ M_{21}(t,\xi )\varphi _{1,n}(\lambda ,\xi
)+M_{22}(t,\xi )\varphi _{2,n}(\lambda ,\xi )\right\} d_{\alpha }\xi
d_{\alpha }t,
\end{eqnarray*
applying successive approximations method to the equations (4) and (5) and
using Lemma3, we get the estimates (6) and (7)$,$
\end{proof}

\bigskip

The characteristic function $\Delta (\lambda )$ of the problem (1)-(3) is
defined by the relation
\begin{equation}
\Delta (\lambda )=\varphi _{1}(\pi ,\lambda )\sin \beta +\varphi _{2}(\pi
,\lambda )\cos \beta ,
\end{equation
It is obvious that $\Delta (\lambda )$ is an entire function and its zeros,
namely $\left\{ \lambda _{n}\right\} _{n\in
\mathbb{Z}
}$ ,\ coincide with the eigenvalues of the problem (1)-(3). Using the
asymptotic formulae (6) and (7), one can easily obtai
\begin{eqnarray}
&&\left. \Delta (\lambda )=\sin \left[ \lambda \dfrac{x^{\alpha }}{\alpha
-\mu (x)-\theta +\beta \right] -\frac{1}{2\lambda }\upsilon (x)\sin \left[
\lambda \dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta -\beta \right] \right.
\medskip  \notag \\
&&-\frac{1}{2\lambda }\upsilon (0)\sin \left[ \lambda \dfrac{x^{\alpha }}
\alpha }-\mu (x)+\theta +\beta \right] \medskip  \notag \\
&&-\frac{1}{2\lambda }\cos \left[ \lambda \dfrac{x^{\alpha }}{\alpha }-\mu
(x)-\theta +\beta \right] \int_{0}^{x}\upsilon ^{2}(t)d_{\alpha }t\medskip \\
&&-\frac{1}{2\lambda }K\left( x\right) \sin \left[ \lambda \dfrac{x^{\alpha
}{\alpha }-\mu (x)-\theta +\beta \right] \medskip  \notag \\
&&+\frac{1}{2\lambda }L\left( x\right) \cos \left[ \lambda \dfrac{x^{\alpha
}{\alpha }-\mu (x)-\theta +\beta \right] +o\left( \frac{1}{\lambda }\exp
(\left\vert \tau \right\vert \dfrac{x^{\alpha }}{\alpha })\right) ,  \notag
\end{eqnarray
for sufficiently large $\left\vert \lambda \right\vert .$ Since the
eigenvalues of the problem (1)-(3) are the roots of $\Delta (\lambda _{n})=0
, we can write the following equation for them
\begin{equation*}
\begin{array}{l}
\left( 1-\dfrac{1}{2\lambda _{n}}\upsilon (\pi )\cos 2\beta -\dfrac{1}
2\lambda _{n}}\upsilon (0)\cos 2\theta -\dfrac{1}{2\lambda _{n}}K(\pi
)\right) \tan (\lambda _{n}\dfrac{\pi ^{\alpha }}{\alpha }-\mu (\pi )-\theta
+\beta )=\medskip \\
-\dfrac{1}{2\lambda _{n}}\upsilon (\pi )\sin 2\beta -\dfrac{1}{2\lambda _{n}
\upsilon (0)\sin 2\theta +\dfrac{1}{2\lambda _{n}}\int_{0}^{\pi }\upsilon
^{2}(t)d_{\alpha }t-\dfrac{1}{2\lambda _{n}}L(\pi )+o\left( \dfrac{1}
\lambda _{n}}\right) \medski
\end{array
\end{equation*
which implies tha
\begin{equation*}
\begin{array}{l}
\tan (\lambda _{n}\dfrac{\pi ^{\alpha }}{\alpha }-\mu (\pi )-\theta +\beta
)=\medskip \\
\left( 1-\dfrac{1}{2\lambda _{n}}\upsilon (\pi )\cos 2\beta -\dfrac{1}
2\lambda _{n}}\upsilon (0)\cos 2\theta -\dfrac{1}{2\lambda _{n}}K(\pi
)\right) ^{-1}\times \medskip \\
\times \left( -\dfrac{1}{2\lambda _{n}}\upsilon (\pi )\sin 2\beta -\dfrac{1}
2\lambda _{n}}\upsilon (0)\sin 2\theta +\dfrac{1}{2\lambda _{n}
\int_{0}^{\pi }\upsilon ^{2}(t)d_{\alpha }t-\dfrac{1}{2\lambda _{n}}L(\pi
)+o\left( \dfrac{1}{\lambda _{n}}\right) \right
\end{array
\end{equation*
for sufficiently large $n.$

We obtain from the last equation
\begin{eqnarray}
&&\left. \lambda _{n}=\dfrac{\alpha }{\pi ^{\alpha -1}}n+\alpha \dfrac
\theta +\mu (\pi )-\beta }{\pi ^{\alpha }}\medskip \right.  \notag \\
&&+\dfrac{\alpha }{2n\pi ^{\alpha }}\left( \upsilon (\pi )\sin 2\beta
-\upsilon (0)\sin 2\theta +\int_{0}^{\pi }\upsilon ^{2}(t)d_{\alpha }t-L(\pi
)\right) \medskip \\
&&+o\left( \dfrac{1}{n}\right)  \notag
\end{eqnarray
for $\left\vert n\right\vert \rightarrow \infty .$

\section{\textbf{Main Results}}

In this section, we obtain the asymptotic formula for the nodal points of
considered problem and prove an inverse nodal problem for the
one-dimensional conformable fractional Dirac-type integro differential
system.

\begin{lemma}
For sufficiently large $n$, the first component $\varphi _{1}(x,\lambda
_{n}) $ of the eigenfunction $\varphi (x,\lambda _{n})$ has exactly $n$
nodes $\left\{ x_{n}^{j}:j=\overline{0,n-1}\right\} $ in the interval
\left( 0,\pi \right) $:\newline
$0<x_{n}^{0}<x_{n}^{1}<...<x_{n}^{n-1}<\pi $. The numbers $\left\{
x_{n}^{j}\right\} $ satisfy the following asymptotic formula
\begin{equation}
\begin{array}{l}
\left( x_{n}^{j}\right) ^{\alpha }=\dfrac{\left( j+1/2\right) \pi ^{\alpha
}{n}+\dfrac{\mu (x_{n}^{j})+\theta }{n\pi ^{1-\alpha }}\medskip \\
-\dfrac{\left( j+1/2\right) \pi ^{\alpha }}{n}\left( \dfrac{\theta +\mu (\pi
)-\beta }{n\pi }\right) -\dfrac{\theta +\mu (\pi )-\beta }{\pi ^{2-\alpha
}n^{2}}\left( \mu (x_{n}^{j})+\theta \right) \medskip \\
+\dfrac{\alpha }{2n^{2}}\left( \upsilon (0)\sin 2\theta
+\int_{0}^{x_{n}^{j}}\upsilon ^{2}(t)d_{\alpha }t-L\left( x_{n}^{j}\right)
\right) +o\left( \frac{1}{n^{2}}\right)
\end{array
\end{equation}
\end{lemma}

\begin{proof}
From (6), the following asymptotic formula can be written for sufficiently
large $n
\begin{eqnarray*}
&&\left. \varphi _{1}(x,\lambda _{n})=\cos \left[ \lambda _{n}\dfrac
x^{\alpha }}{\alpha }-\mu (x)-\theta \right] +\frac{1}{2\lambda _{n}
\upsilon (x)\cos \left[ \lambda _{n}\dfrac{x^{\alpha }}{\alpha }-\mu
(x)-\theta \right] \right. \medskip \\
&&-\frac{1}{2\lambda _{n}}\upsilon (0)\cos \left[ \lambda _{n}\dfrac
x^{\alpha }}{\alpha }-\mu (x)+\theta \right] +\frac{1}{2\lambda _{n}}\sin
\left[ \lambda _{n}\dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta \right]
\int_{0}^{x}\upsilon ^{2}(t)d_{\alpha }t\medskip \\
&&-\frac{1}{2\lambda _{n}}K\left( x\right) \cos \left[ \lambda _{n}\dfrac
x^{\alpha }}{\alpha }-\mu (x)-\theta \right] -\frac{1}{2\lambda _{n}}L\left(
x\right) \sin \left[ \lambda _{n}\dfrac{x^{\alpha }}{\alpha }-\mu (x)-\theta
\right] \medskip \\
&&+o\left( \frac{1}{\lambda _{n}}\exp (\left\vert \tau \right\vert \dfrac
x^{\alpha }}{\alpha })\right) ,
\end{eqnarray*
from $\varphi _{1}(\left( x_{n}^{j}\right) ^{\alpha },\lambda _{n})=0,$ we
ge
\begin{eqnarray*}
&&\left. \cos \left[ \lambda _{n}\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}
\alpha }-\mu (x_{n}^{j})-\theta \right] =-\frac{1}{2\lambda _{n}}\upsilon
(x_{n}^{j})\cos \left[ \lambda _{n}\dfrac{\left( x_{n}^{j}\right) ^{\alpha
}{\alpha }-\mu (x_{n}^{j})-\theta \right] \right. \medskip \\
&&+\frac{1}{2\lambda _{n}}\upsilon (0)\cos \left[ \lambda _{n}\dfrac{\left(
x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu (x_{n}^{j})-\theta \right] \cos
2\theta \medskip \\
&&-\frac{1}{2\lambda _{n}}\upsilon (0)\sin \left[ \lambda _{n}\dfrac{\left(
x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu (x_{n}^{j})-\theta \right] \sin
2\theta \medskip \\
&&-\frac{1}{2\lambda _{n}}\sin \left[ \lambda _{n}\dfrac{\left(
x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu (x_{n}^{j})-\theta \right]
\int_{0}^{x_{n}^{j}}\upsilon ^{2}(t)d_{\alpha }t\medskip \\
&&+\frac{1}{2\lambda _{n}}K\left( x_{n}^{j}\right) \cos \left[ \lambda _{n
\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu (x_{n}^{j})-\theta
\right] \medskip \\
&&+\frac{1}{2\lambda _{n}}L\left( x_{n}^{j}\right) \sin \left[ \lambda _{n
\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu (x_{n}^{j})-\theta
\right] +o\left( \frac{1}{\lambda _{n}}\right) ,
\end{eqnarray*
$\medskip
\begin{eqnarray*}
&&\left. \tan \left[ \lambda _{n}\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}
\alpha }-\mu (x_{n}^{j})-\theta -\frac{\pi }{2}\right] \left( 1+\frac{1}
2\lambda _{n}}\upsilon (x_{n}^{j})-\frac{1}{2\lambda _{n}}\upsilon (0)\cos
2\theta -\frac{1}{2\lambda _{n}}K\left( x_{n}^{j}\right) \right) =\right.
\medskip \\
&&\frac{1}{2\lambda _{n}}\upsilon (0)\sin 2\theta +\frac{1}{2\lambda _{n}
\int_{0}^{x_{n}^{j}}\upsilon ^{2}(t)d_{\alpha }t-\frac{1}{2\lambda _{n}
L\left( x_{n}^{j}\right) +o\left( \frac{1}{\lambda _{n}}\right) ,
\end{eqnarray*
\newline
Taking into account Taylor's expansion formula for the arctangent, we get
\medskip $
\begin{equation*}
\lambda _{n}\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}{\alpha }-\mu
(x_{n}^{j})-\theta -\frac{\pi }{2}=j\pi +\dfrac{1}{2\lambda _{n}}\left(
\upsilon (0)\sin 2\theta +\int_{0}^{x_{n}^{j}}\upsilon ^{2}(t)d_{\alpha
}t-L\left( x_{n}^{j}\right) \right) +o\left( \dfrac{1}{\lambda _{n}}\right)
.\medskip
\end{equation*
It follows from the last equality

\begin{equation*}
\dfrac{\left( x_{n}^{j}\right) ^{\alpha }}{\alpha }=\dfrac{\left( j+\frac{1}
2}\right) \pi +\mu (x_{n}^{j})+\theta }{\lambda _{n}}+\dfrac{1}{2\lambda
_{n}^{2}}\left( \upsilon (0)\sin 2\theta +\int_{0}^{x_{n}^{j}}\upsilon
^{2}(t)d_{\alpha }t-L\left( x_{n}^{j}\right) \right) +o\left( \dfrac{1}
\lambda _{n}^{2}}\right) .
\end{equation*
The relation (11) is proven by using the asymptotic formul
\begin{equation*}
\lambda _{n}^{-1}=\dfrac{\pi ^{\alpha -1}}{n\alpha }\left\{ 1-\dfrac{\mu
(\pi )+\theta -\beta }{n\pi }-\dfrac{\left( \upsilon (\pi )\sin 2\beta
-\upsilon (0)\sin 2\theta +\int_{0}^{\pi }\upsilon ^{2}(t)d_{\alpha }t-L(\pi
)\right) }{2n^{2}\pi }+o\left( \dfrac{1}{n^{2}}\right) \right\}
\end{equation*}
\end{proof}

Let $X$ be the set of nodal points and $\mu (\pi )=0.$ For each fixed $x\in
\left( 0,\pi \right) $ and $\alpha \in (0,1]$ $\ $we can choose a sequence
\left( x_{n}^{j}\right) \subset X$ so that $x_{n}^{j}$ converges to $x.$
Then the following limits are exist and finite
\begin{equation*}
\underset{\left\vert n\right\vert \rightarrow \infty }{\lim }n\left( \left(
x_{n}^{j(n)}\right) ^{\alpha }-\dfrac{\left( j+1/2\right) \pi ^{\alpha }}{n
\right) =f(x),
\end{equation*
wher
\begin{equation}
f(x)=\dfrac{\mu (x)+\theta }{\pi ^{1-\alpha }}-\dfrac{x^{\alpha }}{\pi
\left( \theta +\mu (\pi )-\beta \right)
\end{equation
an
\begin{equation*}
\underset{\left\vert n\right\vert \rightarrow \infty }{\lim }2n^{2}\left(
\left( x_{n}^{j(n)}\right) ^{\alpha }-\dfrac{\left( j+1/2\right) \pi
^{\alpha }}{n}-\dfrac{\mu (x_{n}^{j})+\theta }{n\pi ^{1-\alpha }}+\dfrac
\left( j+1/2\right) \pi ^{\alpha }}{n}\left( \dfrac{\theta +\mu (\pi )-\beta
}{n\pi }\right) \right) =g(x),
\end{equation*
wher
\begin{equation}
g(x)=\alpha \left( \upsilon (0)\sin 2\theta +\int_{0}^{x}\upsilon
^{2}(t)d_{\alpha }t-L\left( x\right) \right)
\end{equation
Therefore, proof of the following theorem is clear. Let $\mu (\pi )=0$

\begin{theorem}
The given dense subset of nodal points $X$ uniquely determines the
coefficients $\theta $ and $\beta $ of the boundary conditions and if $L(x)$
is known, $X$ also uniquely determines the potential $\Omega (x)$ a.e. on
\left( 0,\pi \right) $ . Moreover, $\Omega (x),$ $L(x),$ $\theta $ and
\beta $ can be reconstructed by the following formulae:

\textbf{Step-1:} For each fixed $x\in (0,\pi )$ and $\alpha \in (0,\pi ],$
choose a sequence $\left( x_{n}^{j(n)}\right) \subset X$ such that $\underse
{\left\vert n\right\vert \rightarrow \infty }{\lim }x_{n}^{j(n)}=x;$

\textbf{Step-2: }Find the function $f(x)$ from (12) and calculate
\begin{eqnarray}
\theta &=&f(0)  \notag \\
\beta &=&f(\pi )\theta ^{1-\alpha }  \notag \\
D_{x}^{\alpha }\mu (x) &=&D_{x}^{\alpha }f(x)+\dfrac{\alpha }{\pi }\left[
f(0)-f(\pi )f(0)^{1-\alpha }\right]
\end{eqnarray
\textbf{Step-3: }Find the function $g(x)$ from (13) and calculate
\begin{equation}
\medskip \upsilon (x)=\dfrac{1}{\sqrt{\alpha }}\sqrt{D_{x}^{\alpha
}(g(x)+\alpha L(x))}
\end{equation
\textbf{Step-4: } If $L(x)$ is known then from (14) and (15) calculat
\begin{eqnarray*}
p(x) &=&\upsilon (x)+D_{x}^{\alpha }\mu (x) \\
r(x) &=&D_{x}^{\alpha }\mu (x)-\upsilon (x)
\end{eqnarray*
If $p(x)$ and $r(x)$ are known then from (13) calculat
\begin{equation*}
L\left( x\right) =\upsilon (0)\sin 2\theta +\int_{0}^{x}\upsilon
^{2}(t)d_{\alpha }t-\frac{g(x)}{\alpha }
\end{equation*}
\end{theorem}


\section*{Electronic states in a multi-electron quantum dot}

\subsection*{Single-particle states in the presence of a perpendicular magnetic field}

In this section we investigate the spectrum of a multi-electron quantum dot in the presence of a magnetic field perpendicular to the plane of the dot, where we will largely follow the method used in Ref.~\onlinecite{SEto1997}. We assume strong confinement of the electrons along the $z$-direction, perpendicular to the plane, and a circularly symmetric parabolic in-plane  confinement. Under these assumptions we write a single-particle Hamiltonian in the $xy$-plane for electron $i$:
\begin{align}\label{s:ham0}
H_0^{(i)} = {}&{}  \frac{1}{2 m^*} [\mathbf p_i + e \mathbf A (\mathbf r_i)]^2 + \frac{1}{2} m^* \omega_0^2 \mathbf r_i^2 + \frac{1}{2}g\mu_{\rm B}B \sigma_i^{z},
\end{align}
where $\mathbf A (\mathbf r) = \frac{1}{2}B(x\hat y - y\hat x)$ is the vector potential describing the magnetic field ${\bf B} = B\hat z$, $m^*$ is the effective mass of the electrons, $\omega_0$ defines the strength of the in-plane confinement, such that $\sigma_0=\sqrt{\hbar/(m^* \omega_0)}$ gives the effective radius of the dot in the absence of a magnetic field, and $g$ is the effective $g$-factor of the host material.

The eigenstates of this Hamiltonian are the Fock-Darwin states,
\begin{align}
\psi_{n,l,\eta} (\mathbf r_i) = {}&{} \sqrt{\frac{n!}{\pi \sigma^2 (n+|l|)!}}
\, \rho_i^{|l|} e^{-\rho_i^2/2} L_n^{|l|}\left( \rho_i^2 \right) e^{-i l \theta_i},\label{s:fdstates}
\end{align}
written in terms of polar coordinates $\rho_i = r_i/\sigma$ and $\theta_i$.
Here, $\sigma = \sqrt{\hbar/m^* \Omega}$ is the magnetic-field-dependent effective dot radius, with $\Omega = \sqrt{\omega_0^2 + \omega_c^2/4}$ and $\omega_c = e B/m^*$, and
$L_a^{b}\left( x \right)$ is the associated Laguerre polynomial.
The quantum number $n = 0,1,2,\dots$ labels the radial orbital degree of freedom, the quantum number $l \in \mathbb{Z}$ the orbital angular momentum, and $\eta=\pm 1$ the spin of the electron.
The corresponding eigenenergies are
\begin{align}\label{s:FD}
E_{n,l,\eta}^{(i)} = {}&{} \hbar \Omega (2 n + |l| + 1) - \frac{1}{2} \hbar \omega_c l + \frac{1}{4}g\frac{m^*}{m_e}\hbar\omega_c\eta,
\end{align}
where $m_e$ is the bare electron rest mass.
The first term contributes the regular two-dimensional harmonic-oscillator energies, but with a magnetic-field-dependent oscillator frequency $\Omega$, the second term adds the direct coupling of the angular momentum $l$ to the perpendicular magnetic field, and the last term accounts for the Zeeman effect.
In Fig.~\ref{s:fd} we show this spectrum as a function of $\omega_c$, where the boxed labels indicate the quantum numbers $(n,l)$ of the lowest few states.
For clarity we omitted the Zeeman effect, i.e., we set $g=0$, meaning that all lines are still twofold degenerate.
The levels plotted in red are the ones we used to construct the many-particles states that formed the basis for our analytic and numerical calculations, see below.
\begin{figure}[t]
	\includegraphics[scale=0.75]{fock_darwin}
	\caption{The Fock-Darwin spectrum as a function of $\omega_c = eB/m^*$ for $g=0$, see Eq.~(\ref{s:FD}), where the boxed labels indicate the orbital and angular momentum quantum numbers $(n,l)$.
	The levels plotted in red are the ones we included in constructing the many-particle states that we used as a basis for our calculations.}
	\label{s:fd}
\end{figure}

\subsection*{Many-particle Hamiltonian including interactions}

In a quantum dot with more than one electron one has to account for electron-electron interactions as well, which we describe by the Hamiltonian
\begin{align}\label{s:coulomb}
V = {}&{} \sum_{i<j} \frac{e^2}{4 \pi \varepsilon |\mathbf{r}_i - \mathbf{r}_j|},
\end{align}
where $\varepsilon$ is the effective dielectric constant of the surroundings of the quantum dot.
For a system with two electrons analytical diagonalization of the Hamiltonian $\sum_i H_0^{(i)} + V$ is possible~\cite{SGolovach2008}, but for more than two electrons there is no obvious solution.
We thus treat this many-body problem by working in a restricted configuration space, where we construct a basis of many-particle states from products of single-particle states and impose a cutoff on the quantum numbers $n$ and $l$ as was done in Ref.~\cite{SEto1997}.
In this basis, we can write explicit expressions for the matrix elements of the interaction Hamiltonian (\ref{s:coulomb}).
We then proceed by (\textit{i}) numerical diagonalization of the resulting Hamiltonian matrix and (\textit{ii}) applying perturbation theory in the interaction Hamiltonian $V$, which works best in the weak-coupling limit, where $e^2/4 \pi \varepsilon \sigma_0 \ll \hbar \omega_0$~\cite{SKouwenhoven1997}.

We start by constructing a basis of many-particle states from antisymmetrized products of Fock-Darwin states.
For a system of $M$ electrons one such product state, which we denote $\ket{s}$, is characterized by a set of quantum numbers $s=\{ n_{s_1}, l_{s_1},\eta_{s_1}; n_{s_2}, l_{s_2},\eta_{s_2}; \dots ; n_{s_M}, l_{s_M},\eta_{s_M} \}$.
The antisymmetrized wave function in position space $\langle \mathbf r_1, \mathbf r_2,\dots \mathbf r_M    |   s \rangle = \phi_s(\mathbf r_1,\mathbf r_2,\dots \mathbf r_M)$ can then be written as
\begin{align}
\phi_s(\mathbf r_1,\mathbf r_2,\dots \mathbf r_M) = \mathcal A [\psi_{n_{s_1},l_{s_1},\eta_{s_1}} (\mathbf r_1) \psi_{n_{s_2},l_{s_2},\eta_{s_2}} (\mathbf r_2) \dots \psi_{n_{s_M},l_{s_M},\eta_{s_M}} (\mathbf r_M)],
\end{align}
where $\mathcal A$ is the antisymmetrization operator.

Since $\sum_i H_0^{(i)}$ is diagonal in the basis of these product states, we can write the full Hamiltonian as
\begin{align}\label{s:h1}
H = {}&{} \sum_s \left[\hbar \Omega (2 N_s + K_s + M) - \frac{1}{2} \hbar \omega_c L_s + \frac{1}{2}g\frac{m^*}{m_e}\hbar\omega_c S_s + V_{ss}  \right] \ket s \bra s + \sum_{s\neq r} V_{sr} \ket s \bra r,
\end{align}
where $N_s = \sum_i n_{s_i}$, $K_s = \sum_i |l_{s_i}|$, $L_s = \sum_i l_{s_i}$, $S_s = \frac{1}{2} \sum_i \eta_{s_i}$, and $V_{sr}= \bmk{s}{V}{r}$.
To write an explicit matrix form of this Hamiltonian we therefore need to evaluate the integrals
\begin{align}
V_{sr} = {}&{} \sum_{i<j}^M \int d\mathbf r_1 \cdots d\mathbf r_M\, \phi_s^*(\mathbf r_1,\dots \mathbf r_M) \frac{e^2}{4 \pi \varepsilon |\mathbf{r}_i - \mathbf{r}_j|} \phi_r(\mathbf r_1,\dots \mathbf r_M)
\end{align}
for all sets of quantum numbers $s$ and $r$. Since the Coulomb potential couples electrons pairwise, we only need to evaluate integrals of the form
\begin{align}
v_{n_1,l_1;n_2,l_2;n_3,l_3;n_4,l_4} \equiv \int d\mathbf r_1 d\mathbf r_2 \,
\psi_{n_1,l_1,\eta}^* (\mathbf r_1) \psi_{n_2,l_2,\eta'}^* (\mathbf r_2) \psi_{n_3,l_3,\eta'} (\mathbf r_2) \psi_{n_4,l_4,\eta} (\mathbf r_1)
\frac{e^2}{4 \pi \varepsilon |\mathbf{r}_1 - \mathbf{r}_2|}.
\end{align}
Using the Fock-Darwin states as given in Eq.~(\ref{s:fdstates}) the result of the integral can be written in a closed form~\cite{SRontani1999a},
\begin{align}
v_{n_1,l_1;n_2,l_2;n_3,l_3;n_4,l_4} = {} & {}
\frac{e^2}{4\sqrt 2 \pi\varepsilon\sigma} \delta_{l_1+l_2,l_3+l_4} \sqrt{ \prod_{i=1}^4 \frac{n_i!}{2^{|l_i|}(n_i+|l_i|)!}}
\sum_{j_1=0}^{n_1}\sum_{j_2=0}^{n_2}\sum_{j_3=0}^{n_3}\sum_{j_4=0}^{n_4}
\left[ \prod_{k=1}^4 \frac{(-\frac{1}{2})^{j_k}}{j_k!} \left( \begin{array}{c} n_k+|l_k| \\ n_k - j_k \end{array} \right) \right]
\nonumber\\ {} & {}
\times \sum_{\lambda_1 = 0}^{\alpha_{1}}\sum_{\lambda_2 = 0}^{\alpha_{2}}\sum_{\lambda_3 = 0}^{\alpha_{3}}\sum_{\lambda_4 = 0}^{\alpha_{4}} \delta_{\lambda_1+\lambda_2,\lambda_3+\lambda_4} \left[ \prod_{t=1}^4 \left(\begin{array}{c} \alpha_{t} \\ \lambda_t \end{array}\right)\right]
(-1)^{\alpha_2+\alpha_3-\lambda_2-\lambda_3}\,
\Gamma \left( \frac{\Lambda+2}{2} \right)\Gamma\left( \frac{A-\Lambda +1}{2} \right),
\end{align}
where $\alpha_{i} = j_i+j_{5-i}+\frac{1}{2}(|l_i|+l_i+|l_{5-i}|-l_{5-i})$, $\Lambda = \sum_{i=1}^4 \lambda_i$, and $A = \sum_{i=1}^4 \alpha_i$.
This allows us to find analytic expressions for all $V_{sr}$ in (\ref{s:h1}) and thus to write $H$ in a closed matrix form.

\subsection*{Numerical results}

We will investigate many-particle states with up to 5 electrons in a single quantum dot, assuming relatively small applied magnetic fields, $\omega_c/\omega_0 \leq 0.3$, and not too strong interactions, characterized by the dimensionless parameter $\kappa = e^2/4\pi\varepsilon\sigma_0 \hbar\omega_0 \lesssim 1$.
Since we are only interested in finding the lowest few levels for $M \leq 5$, the 28 single-particle states with $n\in \{ 0,1 \}$, $l\in \{ 0,\pm 1,\pm 2,\pm 3 \}$ and $\eta = \pm 1$ (the levels plotted in red in Fig.~\ref{s:fd}, which all still have a twofold spin degeneracy) form in this case a reasonable set to construct our basis of many-particle states from.

\begin{figure}[t]
	\includegraphics[scale=0.63]{spectrum_supp}
	\caption{The lowest few four-electron levels in a single planar quantum dot as a function of the applied perpendicular magnetic field $\omega_c = eB/m^*$ with $g=-0.4$, for three different interaction parameters: (a) $\kappa=0.2$, (b) $\kappa=0.5$, and (c) $\kappa=1.5$.
	The dots are the result of numerical diagonalization of the Hamiltonian (\ref{s:h1}) written in the basis of antisymmetrized product states of the 28 single-particle states with $n\in \{ 0,1 \}$, $l\in \{ 0,\pm 1,\pm 2,\pm 3 \}$, and $\eta = \pm 1$.
	States that have a singlet spin configuration are plotted in green and states with a triplet spin configuration are plotted in blue.
	The solid lines show the analytic expressions for the field-dependent energies (\ref{s:eigsT}--\ref{s:gammas2}) as found from second-order perturbation theory in $\kappa$.
	The thin red dashed lines in (c) show the result of fitting the numerical data to Eqs.~(\ref{s:fit1}) and (\ref{s:fit2}), which have the same algebraic structure as the perturbative results given in (\ref{s:eigsT}) and (\ref{s:eigsS1}).}
	\label{s:spec}
\end{figure}
We now focus on the case of four electrons and proceed by diagonalizing the Hamiltonian numerically to obtain the eigenstates and -energies of the lowest electronic states in the multi-electron quantum dot.
The results are shown as dots in Fig.~\ref{s:spec}, where we used (a) $\kappa = 0.2$, (b) $\kappa=0.5$, and (c) $\kappa = 1.5$, and we set $g=-0.4$.
Green dots correspond to eigenstates that have a many-particle spin-singlet structure and blue dots to states with a spin-triplet structure.
Under the simplest assumptions that the surroundings of the dot are made of pure GaAs and that there are no additional screening effects, we use $\varepsilon = 12.9\,\varepsilon_0$ and $m^* = 0.067\,m_e$ to estimate the dot size and orbital level splitting corresponding to these two values of $\kappa$, giving $E_{\rm orb} \approx 40~$meV for $\kappa=0.5$ and $E_{\rm orb} \approx 5~$meV for $\kappa=1.5$.
For the most common gate-defined dots $\kappa = 1.5$ seems thus to be more realistic, although we note that the actual value of $\varepsilon$ is probably hard to predict since it can be affected severely by structural inhomogeneities, screening effects due to nearby metallic gates, or the underlying three-dimensional nature of the electronic wave function~\cite{SZhang1986}.

The overall structure of the low-field part of the spectra shown in Fig.~\ref{s:spec} is, however, the same for all values of $\kappa$:
At zero field, exchange effects arising from the Coulomb interaction favor a spin triplet ground state, which for small $\kappa$ has the orbital configuration $(0,0)^2(0,1)^1(0,-1)^1$, where the superscript denotes the number of electrons in the state $(n,l)$~\cite{SEto1997}.
The first excited states are two spin singlets with the configurations $(0,0)^2(0,1)^2$ and $(0,0)^2(0,-1)^2$, and the next excited state is the singlet with the configuration $(0,0)^2(0,1)^1(0,-1)^1$.
When the field is increased, the most pronounced effects are: (\textit{i}) all orbital energies increase due to the magnetic compression of the wave functions, i.e., the dependence of $\Omega$ on $\omega_c$, (\textit{ii}) the first two excited singlets split in energy due to their total projected orbital angular momentum of $L_z = \pm 2\hbar$, and (\textit{iii}) the three triplet states split due to the Zeeman effect.

All together, this leads to a singlet-triplet crossing at $\omega_c/\omega_0 \sim 0.1$ after which a four-particle singlet becomes the ground state.
Close to this degeneracy the next excited state is typically $\sim 0.2\,\hbar\omega_0$ higher in energy, which is far enough to treat the lowest four levels to first approximation as a well-separated subsystem.
One could use levels that cross in this subsystem to construct a singlet-triplet qubit~\cite{SJacak2001a}, but in this case qubit control would still rely on modulation of the magnetic field.
If, instead, we add two more quantum dots with a single electron on each, then we can create a triple-dot exchange-only singlet-only qubit similar to the quadruple-dot qubit proposed in Ref.~\onlinecite{SSala2017}, where the tunability of the singlet-triplet splitting of the two central electrons, through the detuning and coupling between the two central dots, is now replaced by tunability of the splitting through the external magnetic field.

\vspace{6ex}
\subsection*{Analytic results}

Since we have closed-form expressions for all elements of the interaction Hamiltonian, we can do perturbation theory in small $\kappa$ to arrive at analytic expressions for the lowest few eigenenergies in a multi-electron quantum dot~\cite{SRontani1999}.
For the four-electrons case studied above we find for the three lowest triplet levels
\begin{align}\label{s:eigsT}
E_{T_\beta^{(S)}} = 6\hbar\Omega + \frac{S}{2}g\frac{m^*}{m_e}\hbar\omega_c + \gamma^{(1)}_T \kappa \hbar\sqrt{\Omega\omega_0} + \gamma_T^{(2)}\kappa^2\hbar\omega_0,
\end{align}
where $S \in \{-1,0,1\}$ labels the total spin projection of the triplet and we used the coefficients
\begin{align}
	\gamma_T^{(1)} = 2\sqrt{2 \pi}, \qquad
	\gamma_T^{(2)} = -\frac{195893509 \pi}{805306368} \approx -0.764,
\end{align}
that determine the prefactor of the first- and second-order correction, respectively.
The two lowest singlet levels have
\begin{align}\label{s:eigsS1}
	E_{S_{\alpha,\gamma}} =  6\hbar\Omega - \frac{L}{2}\hbar\omega_c + \gamma^{(1)}_{S1} \kappa \hbar\sqrt{\Omega\omega_0} + \gamma_{S1}^{(2)}\kappa^2\hbar\omega_0,
\end{align}
where $L$ labels the total orbital angular momentum projection of the state, i.e., $L = 2$ for the lowest singlet $\ket{S_\alpha}$ and $L = -2$ for the first excited singlet $\ket{S_\gamma}$.
Further,
\begin{align}
	\gamma_{S1}^{(1)} = \frac{67}{16} \sqrt{\frac{\pi}{2} }, \qquad
	\gamma_{S1}^{(2)} = -\frac{38109479 \pi}{134217728} \approx -0.892.
\end{align}
Finally, for the singlet $\ket{S_\beta}$ that lives in the same combination of orbital states as the lowest triplet we find
\begin{align}\label{s:eigsS2}
	E_{S_\beta} =  6\hbar\Omega  + \gamma^{(1)}_{S2} \kappa \hbar\sqrt{\Omega\omega_0} + \gamma_{S2}^{(2)}\kappa^2\hbar\omega_0,
\end{align}
with
\begin{align}
	\gamma_{S2}^{(1)} = \frac{35}{8} \sqrt{\frac{\pi}{2}}, \qquad
	\gamma_{S2}^{(2)} = -\frac{1391260025 \pi}{4294967296} \approx -1.02.\label{s:gammas2}
\end{align}
These results are shown as solid lines in Fig.~\ref{s:spec}(a,b) and match the numerical data reasonably well---better for smaller values of $\kappa$---but for $\kappa = 1.5$ the perturbative results are off by $\sim$~5\%.
Using these results, we can also easily write down an expression for the singlet-triplet splitting denoted $E_{ST}$ in the main text,
\begin{align}
	 E_{T_\beta^0} - E_{S_\alpha} = \hbar\omega_c - \frac{3}{16} \sqrt{\frac{\pi}{2}} \kappa \hbar\sqrt{\Omega\omega_0} + \frac{32763365 \pi}{805306368} \kappa^2\hbar\omega_0,
\end{align}
which for small $\omega_c/\omega_0$ can be very well approximated by
\begin{align}
	E_{T_\beta^0} - E_{S_\alpha} \approx \hbar\omega_c  + \left( - \frac{3}{16} \sqrt\frac{\pi}{2} \kappa + \frac{32763365 \pi}{805306368} \kappa^2\right) \hbar\omega_0
	\approx \hbar\omega_c  + \left( -0.235\,\kappa + 0.128\,\kappa^2\right) \hbar\omega_0.
\end{align}

For larger $\kappa$, such as $\kappa = 1.5$, the level structure of the low-energy part of the spectrum still looks qualitatively very similar to the small-$\kappa$ case, which suggests that we can use the same type of algebraic expression to describe the energies.
If we write for the lowest four levels
\begin{align}
	E_{T_\beta^{(S)}} = {} & {} 6\hbar\Omega + \frac{S}{2}g\frac{m^*}{m_e}\hbar\omega_c + a_{1} \hbar\sqrt{\Omega\omega_0} + b_1 \hbar\omega_0,\label{s:fit1}\\
	E_{S_\alpha} =  {} & {} 6\hbar\Omega - \hbar\omega_c + a_2\hbar\sqrt{\Omega\omega_0} + b_2\hbar\omega_0,\label{s:fit2}
\end{align}
then we can use a least-square fit to the numerical data shown in Fig.~\ref{s:spec}(b) to extract the four parameters
\begin{align}
	a_1 = 7.383, \qquad b_1 = -1.146, \qquad
	a_2 = 7.123, \qquad b_2 = -1.037,
\end{align}
resulting in the red dashed curves shown in the Figure, which agree very well with the numerically calculated results.

\section*{Singlet basis states and Schrieffer-Wolff transformation}

As explained in the main text, the system is modeled using a Hubbard-like Hamiltonian,
\begin{equation}\label{s:hubbard}
{H} = \sum_i^3\left( H_1^{(i)} - V_i{n}_i\right) + \sum_{\langle i,j\rangle}U_c{n}_i{n}_j - \sum_{\langle i,j\rangle,\eta}\frac{t_{ij}}{\sqrt{2}}{c}_{i,\eta}^\dagger{c}_{j,\eta},
\end{equation}
where  $H_1^{(i)}$ is the many-particle Hamiltonian in Eq.~(\ref{s:h1}) acting on the electrons in dot $i$,
the $V_i$ describe the gate-tunable offset voltages on the three dots,
${n}_i=\sum_\eta{c}_{i,\eta}^\dagger{c}_{i,\eta}$ is the electron number operator for dot $i$,
$U_c$ characterizes the interdot electrostatic coupling between the electrons,
the $t_{ij}$ describe tunneling between neighboring dots $i$ and $j$, and ${c}_{i,\eta}$ annihilates an electron on dot $i$ with spin $\eta$.

We assume six electrons occupying a linear array of three quantum dots.
Our qubit is defined in the singlet subspace of the (1,4,1) charge configuration, the lowest four singlet states being
\begin{align}
	\ket 0 = {}&{} \ket{S_{\alpha} S_{(13)}},
	\label{s:ket0}
	\\
	\ket 1 = {}&{} \frac{1}{\sqrt 3}\left[\ket{T_{\beta}^{0} T_{(13)}^{0}} - \ket{T_{\beta}^{-} T_{(13)}^{+}} - \ket{T_{\beta}^{+} T_{(13)}^{-}} \right],
	\label{s:ket1}
	\\
	\ket 2 = {}&{} \ket{S_{\beta} S_{(13)}},
	\label{s:ket2}
	\\
	\ket 3 = {}&{} \ket{S_{\gamma} S_{(13)}},
	\label{s:ket3}
\end{align}
where $\ket{S_{(13)}}$ and $\ket{T_{(13)}}$ indicate pairing in a singlet or triplet state of the two electrons in the outer dots, and $\ket{S_{\alpha,\beta,\gamma}}$ and $\ket{T_{\beta}}$ are the lowest four-particle singlets and triplet on the central dot, see the main text and above.
In the absence of interdot tunneling, i.e., for $t_{12}=t_{23}=0$, the energies of these four states are
\begin{align}
	E_0 = {} & {} 2E_{D_0^{(1,1)}} + E_{S_\alpha^{(4,2)}} - V_1 - 4V_2 - V_3 + 8U_c, \\
	E_1 = {} & {} 2E_{D_0^{(1,1)}} + E_{T_\beta^{(4,2)}} - V_1 - 4V_2 - V_3 + 8U_c, \\
	E_2 = {} & {} 2E_{D_0^{(1,1)}} + E_{S_\beta^{(4,2)}} - V_1 - 4V_2 - V_3 + 8U_c, \\
	E_3 = {} & {} 2E_{D_0^{(1,1)}} + E_{S_\gamma^{(4,2)}} - V_1 - 4V_2 - V_3 + 8U_c.
\end{align}
Here, $E_{D_0^{(n,i)}}$ denotes the lowest doublet eigenenergy of the many-particle Hamiltonian (\ref{s:h1}) for the case of $n$ electrons on dot $i$, where we assumed that dots 1 and 3 are identical.
Note that these energies do not include the Zeeman energy, i.e., they represent the case $E_{\rm Z} = 0$~\cite{SNote1}.
Since the corresponding state $\ket{D_0^{(1,1)}}$ has the exact orbital configuration $(0,0)^1$ this means that we simply have $E_{D_0^{(1,1)}} = \hbar\Omega^{(1)}$, where the superscript (1) indicates that we have to use the orbital energy of dot 1 (we assume the lateral dots to have a slightly smaller size than the central one).
Similarly, $E_{S_\alpha^{(4,2)}}$ is the eigenenergy of the lowest four-particle singlet state on the central dot, such as investigated above, etc.

We would now like to introduce the effect of the tunnel coupling between the dots, to leading order in the coupling parameters $t_{l,r}$ denoting the tunneling coupling between the leftmost and rightmost two dots, respectively.
For that purpose, we need to consider virtual transitions to the neighboring charge configurations (1,3,2), (2,3,1), (1,5,0), and (0,5,1).
(The charge states (2,2,2) and (0,6,0) also border the (1,4,1) region, but a transition to one of these requires \emph{two} tunneling events.)
The states that are directly coupled to the four basis states (\ref{s:ket0}--\ref{s:ket3}) are
\begin{align}
    \ket{4} = {} & {} \ket{\{ D_0^{(1,1)} D_0^{(3,2)}\}_S S_0^{(2,3)} }, \\
    \ket{5} = {} & {} \ket{\{ D_0^{(1,1)} D_1^{(3,2)} \}_S S_0^{(2,3)}}, \\
    \ket{6} = {} & {} \ket{\{ D_0^{(1,1)} D_0^{(5,2)} \}_S}, \\
    \ket{7} = {} & {} \ket{\{ D_0^{(1,1)} D_1^{(5,2)} \}_S},
\end{align}
where we used the same notation as above, i.e., $\ket{S_0^{(2,3)}}$ denotes the ground state singlet formed by two electrons in dot 3, in addition to which we used $\ket{\{D_aD_b\}_S}$ to denote the spin singlet formed by the two doublets $\ket{D_a}$ and $\ket{D_b}$.
We see that states $\ket{4}$ and $\ket{5}$ are (1,3,2) states and $\ket{6}$ and $\ket{7}$ (1,5,0) states, and there are thus four more states, $\ket{8}$--$\ket{11}$, that are exactly the same but with the dot indices 1 and 3 interchanged.
The three- and five-particle doublet ground states $\ket{D_0^{(3,2)}}$ and $\ket{D_0^{(5,2)}}$ have a dominating orbital configuration of $(0,0)^2(0,1)^1$ and $(0,0)^2(0,1)^2(0,-1)^1$, respectively, and since the splitting to the first excited doublet states with main configurations $(0,0)^2(0,-1)^1$ and $(0,0)^2(0,1)^1(0,-1)^2$ is relatively small we included them in the perturbation theory; these first excited doublets are indicated with a subscript 1.
Using the same notation as before, the energies of these virtual states follow straightforwardly as
\begin{align}
	E_4 = {} & {} E_{D_0^{(1,1)}} + E_{D_0^{(3,2)}} + E_{S_0^{(2,3)}} - V_1 - 3V_2 - 2V_3 + 9U_c, \\
	E_5 = {} & {} E_{D_0^{(1,1)}} + E_{D_1^{(3,2)}} + E_{S_0^{(2,3)}} - V_1 - 3V_2 - 2V_3 + 9U_c, \\
	E_6 = {} & {} E_{D_0^{(1,1)}} + E_{D_0^{(5,2)}} - V_1 -5 V_2 + 5U_c, \\
	E_7 = {} & {} E_{D_0^{(1,1)}} + E_{D_1^{(5,2)}} - V_1 - 5 V_2 + 5U_c,
\end{align}
and the energies $E_{8,9,10,11}$ again by interchanging the dot indices 1 and 3.

We now include the tunnel couplings $t_{12} \equiv t_l$ and $t_{23} \equiv t_r$, and assuming that we are deep enough in the (1,4,1) region so that the energy differences to the other four charge states is much larger than the tunnel couplings we can evaluate the exchange effects perturbatively in $t_{l,r}$.
This is done using a Schrieffer-Wolff transformation, and gives up to second order in $t_{l,r}$ the effective (1,4,1) Hamiltonian
\begin{equation}\label{s:sw}
    {H}^{(1,4,1)} = \left(\begin{array}{cccc}
        0 & J_x & J_{0,2} & 0 \\
        J_x & E_{ST} + J_1 - J_0 & J_{1,2} & J_{1,3} \\
        J_{2,0} & J_{2,1} & E_{S_\beta^{(4,2)}} - E_{S_\alpha^{(4,2)}}+J_2-J_0 & J_{2,3} \\
        0 & J_{3,1} & J_{3,2} & E_{S_\gamma^{(4,2)}} - E_{S_\alpha^{(4,2)}}+J_3-J_0
    \end{array}\right),
\end{equation}
where we subtracted $E_0+J_0$ as a constant and defined $E_{ST} \equiv E_{T_\beta^{(4,2)}} - E_{S_\alpha^{(4,2)}}$.
The qubit regime we consider is where $(E_{S_{\beta,\gamma}^{(4,2)}} - E_{S_\alpha^{(4,2)}}) \gg E_{ST}$: the magnetic field is tuned not too far from the $S_\alpha T^0_\beta$-crossing in the four-electron central dot, so that the splitting between $\ket{S_\alpha}$ and $\ket{T_\beta^0}$ ($\sim 10~\mu$eV) is much smaller than the distance to the other two singlets $\ket{S_\beta}$ and $\ket{S_\gamma}$ ($\sim 0.5~$meV).
Since typical exchange energies are $J \sim 1~\mu$eV, we can, to first approximation, neglect the exchange-induced coupling of the qubit to the states $\ket{2}$ and $\ket{3}$; they would lead to small corrections of the order $\sim J^2 / (E_{S_{\beta,\gamma}^{(4,2)}} - E_{S_\alpha^{(4,2)}})$.
A second (more practical) reason to neglect these couplings is that these corrections are $\propto J^2 \propto t^4$, and for consistency one would then also have to perform the original Schrieffer-Wolff transformation to order $t^4$, now including the charge states (2,2,2) and (0,6,0) as well.

We then arrive at the effective qubit Hamiltonian
\begin{align}\label{s:hqb}
	H_\text{qb} = \frac{1}{2}(E_{ST} + J_z)\sigma_z + J_x \sigma_x,
\end{align}
where
\begin{align}
	J_z = {} & {} \ \frac{3t_r^2}{4}\left( \frac{1}{E_1-E_4} + \frac{1}{E_1-E_5}+ \frac{1}{E_1-E_6} + \frac{1}{E_1-E_{7}} \right)
	\nonumber\\&
	+ \frac{3t_l^2}{4}\left( \frac{1}{E_1-E_8} + \frac{1}{E_1-E_9}+ \frac{1}{E_1-E_{10}} + \frac{1}{E_1-E_{11}} \right)\nonumber\\
	{} & {} -\frac{t_r^2}{2} \left( \frac{1}{E_0-E_4} + \frac{1}{E_0-E_6} \right) - \frac{t_l^2}{2}\left( \frac{1}{E_0-E_8} + \frac{1}{E_0-E_{10}} \right),\\
	J_x = {} & {} \frac{\sqrt{3}t_r^2}{4\sqrt{2}}\left( \frac{1}{E_0-E_6} + \frac{1}{E_0-E_4}+ \frac{1}{E_1-E_6} + \frac{1}{E_1-E_{4}} \right)
	\nonumber\\&
	- \frac{\sqrt{3}t_l^2}{4\sqrt{2}}\left( \frac{1}{E_0-E_{8}} + \frac{1}{E_0-E_{10}}+ \frac{1}{E_1-E_{8}} + \frac{1}{E_1-E_{10}} \right).
\end{align}

The largest term in the Hamiltonian (\ref{s:hqb}) is
\begin{align}
	E_{ST} = \hbar\omega_c + (a_1-a_2)\hbar\sqrt{\Omega^{(2)}\omega^{(2)}_0} +(b_1-b_2)\hbar\omega^{(2)}_0,
\end{align}
in terms of the notation of Eqs.~(\ref{s:fit1}--\ref{s:fit2}), where $\Omega^{(2)}$ and $\omega_0^{(2)}$ are the oscillator frequencies (with and without magnetic field) of the central dot.
Here one can use $a_1 = \kappa\gamma_T^{(1)}$, $b_1 = \kappa^2\gamma_T^{(2)}$, $a_2 = \kappa\gamma_{S1}^{(1)}$, and $b_2 = \kappa^2\gamma_{S1}^{(2)}$ for $\kappa \lesssim 0.5$.

The exchange terms are relatively small and their approximate magnitude can be related to the width (and ``height'') $2\Delta$ of the stable (1,4,1) region in terms of the tuning parameters $V_d=\frac{1}{2}\left(V_3-V_1\right)$ and $V_m=\frac{1}{2}\left(V_1+V_3\right)-V_2$, respectively.
We assume that we can neglect the difference between $E_0$ and $E_1$ compared to the splitting to the other eight states.
For a given tuning $(V_m,V_d)$ the splitting to the states $\ket 4$, $\ket 6$, $\ket 8$, and $\ket{10}$ then equals the distance to the corresponding excited charge state in the charge stability diagram.
The energies of the four states that involve one excited orbital state, $\ket 5$, $\ket 7$, $\ket 9$, and $\ket{11}$, are higher in energy by $\hbar\omega_c$.
Assuming that the stable (1,4,1) region is roughly symmetric in $V_d$ and $V_m$, we then arrive at the approximate expressions
\begin{align}
	E_{0,1} - E_4 = {} & {} -\Delta+V_d+V_m,\\
	E_{0,1} - E_5 = {} & {} -\Delta+V_d+V_m-\hbar\omega_c,\\
	E_{0,1} - E_6 = {} & {} -\Delta-V_d-V_m,\\
	E_{0,1} - E_7 = {} & {} -\Delta-V_d-V_m-\hbar\omega_c,\\
	E_{0,1} - E_8 = {} & {} -\Delta-V_d+V_m,\\
	E_{0,1} - E_9 = {} & {} -\Delta-V_d+V_m-\hbar\omega_c,\\
	E_{0,1} - E_{10} = {} & {} -\Delta+V_d-V_m,\\
	E_{0,1} - E_{11} = {} & {} -\Delta+V_d-V_m-\hbar\omega_c,
\end{align}
With these approximations we find
\begin{align}
	J_z = {} & {} \frac{1}{2}
	\left( \frac{t_l^2\Delta}{(V_d-V_m)^2-\Delta^2}
	+ \frac{t_r^2\Delta}{(V_d+V_m)^2-\Delta^2}
	+  \frac{3t_l^2(\Delta+\omega_c)}{(V_d-V_m)^2-(\Delta+\omega_c)^2}
	+ \frac{3t_r^2(\Delta + \omega_c)}{(V_d+V_m)^2-(\Delta+\omega_c)^2}
	\right), \\
	J_x = {} & {} \sqrt{\frac{3}{2}}
		\left( \frac{t_r^2\Delta}{(V_d+V_m)^2-\Delta^2} - \frac{t_l^2\Delta}{(V_d-V_m)^2-\Delta^2}\right).
\end{align}
Assuming approximately equal tunnel couplings $t_{l,r}$ and qubit operation near the ``line'' where $V_d = 0$, we define $t = \frac{1}{2}(t_l+t_r)$ and $\delta t = t_l - t_r$, and expand the exchange energies to leading order in $\delta t$ and $V_d$,
\begin{align}
J_z \approx {} & {} -t^2\left[ \frac{\Delta}{\Delta^2 - V_m^2} + \frac{3(\Delta + \hbar\omega_c)}{(\Delta + \hbar\omega_c)^2 - V_m^2} \right], \\
J_x \approx {} & {} \frac{\sqrt 6 t \Delta}{\Delta^2 - V_m^2} \left[ \delta t + \frac{2 t V_m}{\Delta^2 - V_m^2}V_d \right].
\end{align}

\vspace{4ex}
\section*{Higher-order hyperfine interaction and dephasing}
Hosting the qubit in singlet states only results in having no direct coupling between the qubit states via the hyperfine interaction to lowest order.
Higher-order effects, however, can give rise to energy shifts of the qubit splitting that may lead to qubit dephasing.

We treat the hyperfine interaction between the electron spins and the spins of the many nuclei on a mean-field level, resulting in a Zeeman-like Hamiltonian
\begin{equation}
{H}_\text{hf} = \frac{g\mu_B}{2}\sum_{i}\mathbf{K}_i\cdot {\boldsymbol\sigma}_i,\label{s:hf}
\end{equation}
where $\mathbf{K}_i$ is the effective nuclear field acting on electron $i$ and $(\hbar/2){\boldsymbol\sigma}_i$ is the spin operator for electron $i$.
Due to the tiny nuclear magnetic moments, the nuclear-spin density matrix will be in a high-temperature mixed state for all experimentally relevant temperatures.
This results in random nuclear fields ${\bf K}_i$ that have zero mean and a standard deviation $\sigma_K \sim A/\sqrt N$, where $A$ is an effective material-dependent hyperfine coupling parameter and $N$ the number of nuclear spins the electron is coupled to.
For typical GaAs-based quantum dots $\sigma_K \sim 1$--5~mT.

The main contribution to higher-order hyperfine terms in the qubit subspace comes from other (non-qubit) spin states that are close in energy.
There are three triplets that are energetically close to the qubit states (labeled $\ket{T_{1,2,3}^0}$ in Fig.~3 in the main text), one lying below $\ket 0$ and two in between $\ket 1$ and the quintuplet $\ket Q$.
The energy splitting between the qubit states and the triplets are governed by exchange effects, where the splitting between $\ket{0}$ and $\ket{T_1^0}$ goes as $\sim t^4/\Delta^3$ and the splitting between $\ket{1}$ and $\ket{T_2^0}$ and $\ket{T_3^0}$ as $\sim t^2/\Delta$.
Thus, the hyperfine-induced shift in the qubit splitting is dominated by the coupling between $\ket{0}$ and the triplet states
\begin{align}
	\ket{T_1^{+}} = {} & {} \ket{S_\alpha T^{+}_{(13)}}, \\
	\ket{T_1^{0}} = {} & {} \ket{S_\alpha T^{0}_{(13)}}, \\
	\ket{T_1^{-}} = {} & {} \ket{S_\alpha T^{-}_{(13)}}.
\end{align}

We project the hyperfine Hamiltonian (\ref{s:hf}) to the subspace $\{\ket{0},\ket{T_1^{+}},\ket{T_1^{0}},\ket{T_1^{-}}\}$, yielding
\begin{equation}
{H}_\text{hf}
= \left(\begin{array}{cccc}
0 & \iota^-_{13} & \kappa_{13} & \iota^+_{13} \\
\iota^+_{13} & 0 & 0 & 0 \\
\kappa_{13} & 0 & 0 & 0 \\
\iota^-_{13} & 0 & 0 & 0
\end{array}\right),
\end{equation}
where we defined the gradients
\begin{align}
	\kappa_{13} = {} & {} \frac{g\mu_B}{2}\left(K^z_1 - K^z_3\right), \\
	\iota_{13}^\pm = {} & {} \frac{g\mu_B}{2\sqrt{2}}\left(K^\pm_1 - K^\pm_3\right)
	= \frac{g\mu_B}{2\sqrt{2}}\left(K^x_1 \pm iK^y_1 - K^x_3 \mp iK^y_3\right).
\end{align}
Using perturbation theory, we find that, to leading order in $K_i^{x,y,z}$ the energy shift of $\ket{0}$ is
\begin{align}
	\delta E_0 = {} & {} \frac{\kappa_{13}^2}{E_0-E_{T_1^0}} + \frac{\iota_{13}^+\iota_{13}^-}{E_0-E_{T_1^{+}}} + \frac{\iota_{13}^+\iota_{13}^-}{E_0-E_{T_1^{-}}},
\end{align}
where the $E_\nu$ are the unperturbed energy levels.
Due to the larger separation between $\ket 1$ and $\ket{T_{2,3}^0}$, the hyperfine-induced shift of $\ket{1}$ is much smaller.
The shifts caused by the coupling to the polarized triplets $\ket{T_1^{\pm}}$, i.e.\ the two last terms, can be reduced by tuning the Zeeman energy, which can be done independently from $\omega_c$ by tilting the total externally applied field.
The energy of $\ket{T_1^0}$ is not affected by the Zeeman effect, and the shift caused by the coupling to this state is thus solely determined by exchange effects.

To obtain a very rough estimate for the scale of the dephasing time caused by these higher-order hyperfine fields, we consider their effect on the Rabi oscillations when the system is driven resonantly.
Following the approach of Ref.~\cite{SSala2017} we find the estimate $T_2^* \sim A_q\hbar(E_0-E_{T_1^0})^2/\sigma_{K}^4$, which indeed predicts a shorter dephasing time the closer the state $\ket{T_1^0}$ is to $\ket{0}$, i.e., a shorter hyperfine-induced $T_2^*$ at the sweet spot than in the resonant-exchange regime.

To illustrate, we calculate numerically the probability $\langle P_1(\tau)\rangle$ of finding the qubit in $\ket 1$ after initializing in $\ket 0$ and resonantly driving either $V_d$ or $\delta t$ for a time $\tau$, focusing on the two cases illustrated in Fig.~3(a,b) in the main text.
In Fig.~\ref{fig:supp:dephasing} we show the resulting time-dependent probabilities for (a) driving $V_d$ in the resonant-exchange regime, and (b) driving $\delta t$ at the sweet spot, after averaging over 2500 random nuclear configurations with $g\mu_B K_i^{x,y,z}$ each taken from a normal distribution with mean zero and $\sigma_{K} = 0.07\ \mu$eV (corresponding to 3~mT for $|g| = 0.4$).
We used the same parameters as in the main text, with driving amplitudes $\tilde V_d=10\ \mu$eV and $\delta\tilde t=2\ \mu$eV.
The results show that the hyperfine-induced dephasing, even at the sweet spot, is small compared to the Rabi period.

We finally compare this with the rough estimate $T_2^* \sim A_q\hbar(E_0-E_{T_1^0})^2/\sigma_{K}^4$.
Using the fact that the Rabi period is given by $T_{\rm Rabi} = h/A_q$, we find for the approximate number of coherent Rabi oscillations that should be visible
\begin{align}
	n_{\rm coh} \equiv \frac{T_2^*}{T_{\rm Rabi}} \sim \frac{1}{2\pi} \frac{A_q^2(E_0-E_{T_1^0})^2}{\sigma_{K}^4} =
	\frac{1}{2\pi} \frac{h^2(E_0-E_{T_1^0})^2}{T_{\rm Rabi}^2\sigma_{K}^4}.
\end{align}
For the resonant-exchange regime of Fig.~\ref{fig:supp:dephasing}(a) we read off $T_{\rm Rabi}\approx 25$~ns and we find $E_0-E_{T_1^0} \approx 0.80~\mu$eV which yields $n_{\rm coh} \sim 120$, whereas at the sweet spot, see Fig.~\ref{fig:supp:dephasing}(b), we have $T_{\rm Rabi}\approx 19$~ns and $E_0-E_{T_1^0} \approx 0.32~\mu$eV giving $n_{\rm coh} \sim 30$.
We see that these estimates are indeed roughly agreeing with the dephasing observed in Fig.~\ref{fig:supp:dephasing}.

\begin{figure}
	\centering
	\includegraphics[width=.75\textwidth]{dephasing.pdf}
	\caption{Time dependent probability $\langle P_1(\tau)\rangle$ after averaging over 2500 random nuclear field configurations, (a) in the resonant-exchange regime and (b) at the sweet-spot.
	We used the same parameters as in the main text, with $\tilde V_d = 10\ \mu$eV, $\delta\tilde t= 2\ \mu$eV, and $g\mu_B K_i^{x,y,z}$ taken from a normal distribution with zero mean and $\sigma_{K} = 0.07\ \mu$eV.}
	\label{fig:supp:dephasing}
\end{figure}

\end{NoHyper}