# Machine-learned tuning of artificial Kitaev chains
from tunneling-spectroscopy measurements

Jacob Benestad<EMAIL_ADDRESS>Center for Quantum Spintronics,
Department of Physics, Norwegian University of Science and Technology, 7491
Trondheim, Norway  Athanasios Tsintzis  Division of Solid State Physics and
NanoLund, Lund University, S-22100 Lund, Sweden  Rubén Seoane Souto  Instituto
de Ciencia de Materiales de Madrid, Consejo Superior de Investigaciones
Científicas (ICMM-CSIC), Madrid, Spain  Division of Solid State Physics and
NanoLund, Lund University, S-22100 Lund, Sweden  Martin Leijnse  Division of
Solid State Physics and NanoLund, Lund University, S-22100 Lund, Sweden  Evert
van Nieuwenburg  Lorentz Institute and Leiden Institute of Advanced Computer
Science, Leiden University, P.O. Box 9506, 2300 RA Leiden, The Netherlands
Jeroen Danon  Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, 7491 Trondheim, Norway

###### Abstract

We demonstrate reliable machine-learned tuning of quantum-dot-based artificial
Kitaev chains to Majorana sweet spots, using the covariance matrix adaptation
algorithm. We show that a loss function based on local tunnelling-spectroscopy
features of a chain with two additional sensor dots added at its ends provides
a reliable metric to navigate parameter space and find points where crossed
Andreev reflection and elastic cotunneling between neighbouring sites balance
in such a way to yield near-zero-energy modes with very high Majorana quality.
We simulate tuning of two- and three-site Kitaev chains, where the loss
function is found from calculating the low-energy spectrum of a model
Hamiltonian that includes Coulomb interactions and finite Zeeman splitting. In
both cases, the algorithm consistently converges towards high-quality sweet
spots. Since tunnelling spectroscopy provides one global metric for tuning all
on-site potentials simultaneously, this presents a promising way towards
tuning longer Kitaev chains, which are required for achieving topological
protection of the Majorana modes.

##  I Introduction

Realization and control of Majorana bound states (MBSs), which remain much
sought after due to their non-local and topological nature, are a topic of
intense study within condensed matter physics [1, 2, 3, 4, 5, 6]. The simple
Kitaev model gives one example of a topological material where MBSs can be
found [7], and recreating the ingredients of the Kitaev model with engineered
hybrid devices is a very promising avenue in this pursuit [8]. While a lot of
work has focused on making artificial chains using semiconductor nanowires
proximitized by superconductors [9, 10, 11, 12, 13, 14, 15, 16], device
imperfections, such as material disorder, seem currently to be a major
challenge for further progress in this direction [17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29].

An alternative is to rather use arrays of alternating normal and proximitized
quantum dots to build an artificial Kitaev chain, which could circumvent some
of the material disorder issues [30, 31, 32, 33, 34, 35, 36, 37]. Signatures
of so-called “Poor Man’s Majoranas” (PMMs), which are non-local but lack
topological protection, have been demonstrated in a two-site minimal
artificial Kitaev chain [38, 39, 40]. An important challenge faced in this
field is to make and tune devices that are long enough that they can host
topologically protected Majoranas [30]. Recent progress along these lines
includes control over three-site devices [41] and the observation of
signatures of tuning to a Majorana sweet spot in three-site devices [42].

To tune these artificial Kitaev chains to a sweet spot, where an even–odd
ground-state degeneracy coincides with maximal Majorana quality (which we will
define later) of the part of the corresponding excitation that lives on the
outer sites of the chain, it is necessary to tune the system parameters such
that crossed Andreev reflection (CAR) and elastic cotunneling (ECT) processes
via the proximitized dots are pairwise balanced in each link in the chain [43,
44]. In experiment, one faces the challenge of imperfect information about
said tuning parameters, due to microscopic details such as material
imperfections and varying gate lever arms. In practice, tuning has so far been
done by working with pairs of sites at a time, since there are too many
variables for navigating a global metric manually, where all dots are tuned
simultaneously. Such sequential manual tuning across the chain can be time-
consuming, and at the same time one must be careful that the overall tuning
remains correct each time a new pair is tuned [41, 42]. Pairwise tuning
becomes especially precarious for the case of strong tunnel couplings, where
multi-dot renormalization effects can be expected. With single CAR and ECT
processes thus depending on parameters all along the chain, it will be
difficult to tune the system while only considering pairs of dots. Thus,
tuning all dots simultaneously would be advantageous, and the fact that this
would be very difficult for a human to do motivates us to investigate machine-
learning methods to perform automated tuning.

Figure 1: Sketch of the artificial Kitaev chain systems considered in this
paper. In addition to the sensor dots marked $L/R$ on either side of the
Kitaev chain, the array consists of (a) $3$ dots or (b) $5$ dots, where every
other dot is proximitized by a superconductor (blue square) to create an
Andreev bound state that can mediate CAR and ECT processes between the normal
dots.

Such automatic tuning has already been proposed as a way to mitigate disorder
in bulk hybrid nanowire sytems, by using an Aharonov–Bohm interferometer to
get a metric of the Majorana quality that can be used as a global loss
function [45]. A simpler and more practical probe of Majorana quality could be
to perform a tunneling spectroscopy experiment while adding an extra sensor
dot at the end of the system [9, 46, 47], which was recently investigated in
the context of artifical Kitaev chains based on quantum dot arrays [48]. Work
has also been done on the machine-learned tuning of dot-based artificial
Kitaev chains, in the context of using generative neural networks to
automatically classify avoided crossings in charge stability diagrams [49].

In this paper, we demonstrate automatic tuning in simulated two- and three-
site artificial Kitaev chains using the covariance matrix adaptation
evolutionary strategy (CMA-ES) [50, 51, 52] with tunneling spectroscopy
through sensor dots at both ends of the chain as a metric to navigate
parameter space. We first apply the tuning algorithm to a minimal (two-site)
Kitaev chain, where locations of the sweet spots can be relatively easily
determined by numerical calculations for a given set of tight-binding
parameters [33], and we show how the algorithm indeed converges towards these
points. We then consider a three-site chain, where the locations of the sweet
spots are not known a priori. Also in this case, the algorithm converges to
points with vanishing even–odd ground-state splitting and very high Majorana
quality. We thus demonstrate the capability of the CMA-ES algorithm to
navigate a large tuning parameter space and efficiently locate Majorana sweet
spots, based on simple tunneling-spectroscopy data.

The rest of this paper is structured as follows: In Section II we introduce
the Hamiltonian for our quantum-dot-based artificial Kitaev model, and we
explain the basic working of the CMA-ES algorithm used for the automatic
tuning. In Section III we present our results for the automatic tuning, first
for a two-site model in Section III.1, and a three-site model in Section
III.2. Finally, we present a short conclusion in Section IV.

##  II Model

We consider a similar setup to that presented in Refs. [48, 34]: a linear
array of quantum dots where every other dot is proximitized by a
superconductor, and with the addition of extra sensor dots at each end of the
chain. To make the optimization more robust, we probe the low-energy physics
of the system first with the sensor dot added to the left side of the array
and then a second time with the sensor dot on the right side. In an
experiment, such freedom of changing sides for the sensor dot could be
achieved by defining extra dots on both ends of the Kitaev chain and then
“removing” the one that is not being used by lowering its outer barrier such
that the dot becomes incorporated as a part of the lead.

We describe the full array of dots with the Hamiltonian

$H^{(n,s)}=H_{\text{Kit}}^{(n)}+H_{\text{S}}^{(n,s)},$ (1)

where $s=L,R$ indicates on which side the sensor dot is located and $n$ labels
the number of “Kitaev-chain” sites, i.e., the number of normal (non-
proximitized) quantum dots.

The first term in Eq. 1 describes the Kitaev-chain part of the system, which
includes normal dots and dots proximitized by a superconductor, together with
spin-conserving and spin-non-conserving hopping between them,

$\displaystyle H_{\text{Kit}}^{(n)}$
$\displaystyle=\sum_{j=1}^{n}\sum_{\sigma}\varepsilon_{j\sigma}d_{j\sigma}^{\dagger}d_{j\sigma}+\sum_{j=1}^{n}\frac{U}{2}n_{j}(n_{j}-1)$
(2)
$\displaystyle+\sum_{j=1}^{n-1}\Big{[}\varepsilon_{Sj}n_{Sj}+\Delta\big{(}d_{Sj\uparrow}^{\dagger}d_{Sj\downarrow}^{\dagger}+d_{Sj\downarrow}d_{Sj\uparrow}\big{)}\Big{]}$
$\displaystyle+\sum_{j=1}^{n-1}\sum_{\sigma}\Big{[}t_{j,Sj}\,d_{j\sigma}^{\dagger}d_{Sj\sigma}+t_{Sj,j+1}\,d_{Sj\sigma}^{\dagger}d_{(j+1)\sigma}$
$\displaystyle+\beta_{\sigma}\big{(}t_{j,Sj}^{\text{SO}}\,d_{j\sigma}^{\dagger}d_{Sj\bar{\sigma}}+t_{Sj,j+1}^{\text{SO}}\,d_{Sj\sigma}^{\dagger}d_{(j+1)\bar{\sigma}}\big{)}+\text{H.c.}\Big{]},$

where $d_{j\sigma}^{\dagger}$ ($d_{Sj\sigma}^{\dagger}$) is the creation
operator for an electron with spin $\sigma=\uparrow,\downarrow$ on
(proximitized) dot $j$ ($Sj$), and
$n_{j}=d_{j\uparrow}^{\dagger}d_{j\uparrow}+d_{j\downarrow}^{\dagger}d_{j\downarrow}$
($n_{Sj}=d_{Sj\uparrow}^{\dagger}d_{Sj\uparrow}+d_{Sj\downarrow}^{\dagger}d_{Sj\downarrow}$)
is the number operator on that dot. The normal dots are described by the first
line in Eq. 2. Their on-site single-particle energy
$\varepsilon_{j\sigma}=\varepsilon_{j}+E_{Z}\,\delta_{\sigma,\downarrow}$ is
spin-dependent, where $E_{Z}$ is the Zeeman energy and $\varepsilon_{j}$ can
be tuned by nearby electrostatic gates $V_{j}$, and the charging energy
associated with a doubly occupied normal dot is given by $U$. The second line
in Eq. 2 describes the proximitized dots, labeled $Sj$. The proximity of the
superconductor is expected to result in a strong reduction of the effective
on-site $g$-factors, which we incorporate by using spin-independent potentials
$\varepsilon_{Sj}$, for simplicity. The superconductor will also efficiently
screen the on-site Coulomb interactions on these dots, which we thus neglect
here 111Including small Coulomb interactions on the proximitized dots does not
change the qualitative results [33]. The proximity-induced pairing potential
$\Delta$ is assumed to be real and equal for all proximitized dots. Finally,
the last two lines in Eq. 2 describe the coupling between the dots. The
parameters $t_{i,j}$ ($t_{i,j}^{\text{SO}}$) set the amplitude of spin-
conserving (spin-non-conserving) hopping between dots $i$ and $j$, where
$\beta_{\uparrow,\downarrow}=\pm 1$.

The second term in Eq. 1 describes the sensor dots and hopping between the
sensors and the Kitaev chain,

$\displaystyle H_{\text{S}}^{(n,s)}$
$\displaystyle=\sum_{\sigma}\varepsilon_{s\sigma}d_{s\sigma}^{\dagger}d_{s\sigma}+\frac{U}{2}n_{s}(n_{s}-1)$
(3)
$\displaystyle+\delta_{s,L}\sum_{\sigma}\Big{[}t_{L,1}\,d_{L\sigma}^{\dagger}d_{1\sigma}+\beta_{\sigma}t_{L,1}^{\text{SO}}\,d_{L\sigma}^{\dagger}d_{1\bar{\sigma}}+\text{H.c.}\Big{]}$
$\displaystyle+\delta_{s,R}\sum_{\sigma}\Big{[}t_{n,R}\,d_{n\sigma}^{\dagger}d_{R\sigma}+\beta_{\sigma}t_{n,R}^{\text{SO}}\,d_{n\sigma}^{\dagger}d_{R\bar{\sigma}}+\text{H.c.}\Big{]},$

using similar notation as above. In this case,
$\varepsilon_{s\sigma}=\varepsilon_{s}+E_{Z}\,\delta_{\sigma,\downarrow}$
while the Kronecker delta functions $\delta_{s,L}$ and $\delta_{s,R}$ choose
which sensor dot is included.

In the following, we will limit ourselves to the cases where the number of
sites is either $n=2$ or $n=3$. For the case with $n=2$, locations of the
sweet spots are known for the model we use [33] and we can thus use this case
as a test bed to assess the basic working of the tuning algorithm. For the
case with $n=3$, we expect a larger number of sweet spots, the locations of
which are not known a priori. This case will be used to investigate the
robustness of the algorithm when operating in a higher-dimensional parameter
space.

In an experiment, one would connect the chain to a source and drain lead (also
indicated in Fig. 1) and measure the differential conductance across the array
of dots as a function of $\varepsilon_{L}$ or $\varepsilon_{R}$ and the bias
voltage applied to the source and drain. The location of the lowest-bias
conductance peak then reveals the energy splitting between the lowest even and
odd parity states of the array, and its dependence on the sensor dot potential
contains information about the Majorana quality of the states on the outer
sites of the Kitaev chain [48, 46, 47]. The deviation from the ideal case
(zero splitting independent of $\varepsilon_{L}$ or $\varepsilon_{R}$) can
then be quantified into a loss function used by the tuning algorithm to
navigate parameter space, as we will detail in Sec. III. In our numerical
simulations, it suffices to diagonalize $H^{(n,s)}$ and extract the even–odd
ground state splitting directly from its eigenenergies as this corresponds to
the location of the lowest-bias peak [48], which saves considerable
computational overhead as compared to simulating a full differential-
conductance measurement.

To find the optimal set of tuning parameters $\varepsilon_{j,Sj}$ we use the
CMA-ES algorithm 222We use the cma Python library from pypi.org/project/cma,
which is a simple evolution strategy optimization technique, relying on
stochastic sampling rather than derivatives, making it suitable for problems
where the objective function to be minimized for instance has several local
minima or is non-smooth [51, 52, 50]. In our case, the limitation in sharpness
of the conductance curves in real tunneling-spectroscopy measurements due to
broadening from temperature and tunnel coupling to the leads [55] was what
motivated us to use a stochastic and derivative-free method: we can expect the
limited resolution of an actual experiment to make it difficult to accurately
evaluate derivatives of the loss function, likely rendering derivative-based
unsuitable for optimization of the loss. Additionally, we are dealing with
metrics containing several minima, where a purely derivative-based method may
get lost in a non-optimal local minimum.

The CMA-ES algorithm performs the optimization by iterative sampling from and
updating of a multivariate normal (MvN) distribution. For each generation $g$,
a number of $\lambda\geq 2$ offspring $\boldsymbol{\varepsilon}^{(g)}$, which
in this case are configurations of the gate voltages to tune, are randomly
drawn from a MvN distribution with mean $\boldsymbol{\mu}^{(g)}$ and
covariance matrix $\Sigma^{(g)}=(\sigma^{(g)})^{2}C^{(g)}$. The configurations
$\boldsymbol{\varepsilon}^{(g)}$ and means $\boldsymbol{\mu}^{(g)}$ are thus
($2n-1$)-dimensional vectors, e.g., for the two-site Kitaev chain
$\boldsymbol{\varepsilon}^{(g)}=\\{\varepsilon_{1}^{(g)},\varepsilon_{S1}^{(g)},\varepsilon_{2}^{(g)}\\}$.
The offspring are then ranked according to their “fitness” using an objective
(loss) function ${\cal L}(\boldsymbol{\varepsilon}^{(g)})$, and the MvN
distribution is updated by calculating a new position $\mu^{(g+1)}$,
“stepsize” $(\sigma^{(g+1)})^{2}$ and “search direction” $C^{(g+1)}$, based on
the $\eta$ best-ranked offspring and the previous history of changes to the
population. A diagram of the CMA-ES optimization loop is shown in Fig. 2. We
will consistently use $\lambda=56$ and $\eta=28$, based mainly on the
computational resources available, and set all other fine-tuning parameters of
the algorithm to the standard recommended values [50]. Specifically, the
algorithm achieves convergence if (i) the range of best loss values across the
last $10+\lceil{30(2n-1)/\lambda}\rceil$ generations or the range of all loss
values in the last generation is smaller than $10^{-12}$ or (ii) if the
standard deviation of all tuning parameters $\varepsilon_{j,Sj}$ becomes
smaller than $10^{-12}\sigma^{(0)}$ [50]. We also terminate the algorithm if
it exceeds $200$ generations.

Figure 2: Diagram of the CMA-ES algorithm interfaced to either a simulation or
an actual experiment where local differential conductances are measured. Each
generation can be divided into four stages: The first stage is to sample
$\lambda$ realizations from the MvN distribution to generate parameter
configurations. In the second stage, each of these configurations is used in
an experiment, where the dependence of the location of the low-bias
conductance peaks (or the numerically determined lowest even–odd splittings)
on the sensor dot voltage provides a measure for the population fitness. The
third stage is to rank all the parameter configurations according to this
fitness. In the fourth stage, the MvN mean value and covariance matrix are
updated based on a subset of the $\eta$ best-ranked candidates.

##  III Results

For the tuning of an artificial Kitaev chain, it has been suggested that
tracking the splitting between the lowest conductance peaks while sweeping the
sensor dot level $\varepsilon_{s}$ can provide a metric for the Majorana
quality [48]. Following this idea, we will investigate the performance of the
CMA-ES algorithm using the loss function

$f(\boldsymbol{\varepsilon}^{(g)})=\frac{1}{2}\sum_{s=L,R}\max_{\varepsilon_{s}}\left|\frac{E_{0,n,s}^{\text{odd}}-E_{0,n,s}^{\text{even}}}{\Delta}\right|,$
(4)

where $E_{0,n,s}^{\text{odd}}$ and $E_{0,n,s}^{\text{even}}$ are the lowest
eigenenergies of $H^{(n,s)}$ with total odd and even number of electrons,
respectively. The loss function $f(\boldsymbol{\varepsilon}^{(g)})$ thus
returns the maximum of the absolute energy splitting between the even and odd
ground states over a sweep of the sensor dot potential, averaged over the
cases where the sensor dot is attached to the left or right side of the chain.
The sweep range used for $\varepsilon_{s}$ should cover both spin-species
resonances, the locations of which could, e.g., be found roughly from
estimates of $E_{Z}$, $U$ and $\Delta$ [46].

In addition, we make use of an $L_{2}$ regularization to penalize on-site
potentials that grow large, as the amplitudes of CAR and ECT are inversely
proportional to the energy difference between the bound states on neighboring
dots [43] and the magnitude of these amplitudes ultimately sets the scale of
the topological gap in longer chains. Combining $L_{2}$ regularization with
the objective function Eq. 4, the final loss function we use to rank candidate
parameter configurations is

$\mathcal{L}(\boldsymbol{\varepsilon}^{(g)})=f(\boldsymbol{\varepsilon}^{(g)})+\alpha\bigg{\lVert}\frac{\boldsymbol{\varepsilon}^{(g)}}{\Delta}\bigg{\rVert}_{2},$
(5)

where $\alpha$ parameterizes the $L_{2}$ regularization. For such
regularization, one needs to establish some “origin” in parameter space, to
determine $\lVert\boldsymbol{\varepsilon}^{(g)}\rVert_{2}$. Identifying in an
experiment the true origin where $\boldsymbol{\varepsilon}^{(g)}=0$ can be
challenging, but small deviations from this point are not expected to pose an
immediate problem as long as the origin used is not too far away from high-
quality sweet spots.

Throughout, we will set all tunneling amplitudes to $t_{i,j}=0.5\,\Delta$ and
$t^{\text{SO}}_{i,j}=0.1\,\Delta$, except for the ones describing the coupling
to the sensor dot, where we use $t_{L,1}=t_{n,R}=0.25\,\Delta$ and
$t^{\text{SO}}_{L,1}=t^{\text{SO}}_{n,R}=0.05\,\Delta$. The Zeeman energy is
set to $E_{Z}=1.5\,\Delta$ and the onsite Coulomb interaction strength for the
normal dots to $U=5\,\Delta$. These values are on purpose chosen the same as
in Refs. [33, 48], allowing for a straightforward comparison of the results.
The mean of the initial distribution for the CMA-ES algorithm is always set to
the origin, with an initial covariance matrix of $C^{(0)}=\mathbb{I}$, and
$\sigma^{(0)}=0.5\,\Delta$.

In all instances, we let the algorithm search for the optimal set of tuning
parameters until it converges, using the loss function given in Eq. (5). We
then assess the resulting set of tuning parameters
$\boldsymbol{\varepsilon}^{\text{fin}}$ by evaluating (i) the resulting
even–odd ground state splitting in the Kitaev chain (without any sensor dot
attached) $\delta E_{eo}=|E^{\text{odd}}_{0,n}-E^{\text{even}}_{0,n}|$, as
well as (ii) the Majorana polarization (MP), given by [56, 57, 58, 34]

$M_{j}=\frac{\sum_{\sigma}(w_{j\sigma}^{2}-z_{j\sigma}^{2})}{\sum_{\sigma}(w_{j\sigma}^{2}+z_{j\sigma}^{2})},$
(6)

where

$\displaystyle w_{j\sigma}={}$
$\displaystyle{}\bra{O}d_{j\sigma}+d_{j\sigma}^{\dagger}\ket{E},$ (7)
$\displaystyle z_{j\sigma}={}$
$\displaystyle{}\bra{O}d_{j\sigma}-d_{j\sigma}^{\dagger}\ket{E},$ (8)

where $\ket{E}$ and $\ket{O}$ are the even and odd ground states of
$H_{\text{Kit}}^{(n)}$, respectively, and $E^{\text{even}}_{0,n}$ and
$E^{\text{odd}}_{0,n}$ the corresponding ground-state energies. With this
definition $-1\leq M_{j}\leq 1$, and $|M_{j}|=1$ corresponds to having a well-
behaved single Majorana mode localized at site $j$. In Kitaev chains with
finite Zeeman energy and Coulomb interactions, there exist in general no
points in parameter space that have $\delta E_{eo}=0$ and $|M_{1}|=|M_{n}|=1$,
which would indicate the existence of two “perfect” PMMs occupying the outer
dots of the chain. Instead, it is common to look for sweet spots that maximize
$|M_{1,n}|$ on the manifold in parameter space where $\delta E_{eo}=0$ [33,
34]. Below we will thus use $\delta E_{eo}$ and the average MP
$M\equiv(|M_{1}|+|M_{n}|)/2$ as measures for the success of the tuning
procedure.

###  III.1 Two-site Kitaev chain

Figure 3: (a) The difference between even and odd ground state energy of the
two-site Kitaev chain without a sensor dot attached and (b) the average
Majorana polarization at its outer sites as a function of $\varepsilon_{1,2}$
and $\varepsilon_{S1}$, cf. Ref. [33]. There are two sweet spots within the
given parameter range, indicated by the black stars. (c,d) The $60$ first
generations for two example runs of the CMA-ES algorithm for (c) a case that
converged to the left sweet spot shown in (a,b) and (d) a case that found the
right sweet spot. Solid lines show the population mean of the parameters and
the shaded regions give their marginal standard deviation.

As a first benchmark, we test the optimization routine on a simulated minimal
(two-site) Kitaev chain with $n=2$, where the sweet spots for our set of
system parameters have been numerically calculated before [33], to see if it
can reproduce the same results.

In Fig. 3(a,b) we show the even–odd ground state splitting $\delta E_{eo}$ and
average polarization $M$ as a function of $\varepsilon_{1,2}$ and
$\varepsilon_{S1}$, calculated numerically from $H_{\text{Kit}}^{(2)}$. As
expected, these plots are identical to the ones presented in Ref. [33], since
we use the same model and same parameters. The two sweet spots identified in
[33] are indicated with the two stars, and have $\delta E_{eo}=0$ and
$M=0.986$.

Fig. 3(c,d) display the first $60$ generations of two example runs out of in
total $50$ simulations for the CMA-ES optimization algorithm, using the loss
function given in Eq. 5 while setting $\alpha=0$. In the case of the minimal
Kitaev chain, the objective function landscape is simple enough that $L_{2}$
regularization is unnecessary; including a finite $\alpha$ would increase the
probability for the algorithm to discover the left sweet spot, as the right
sweet spot is further away from the origin $\bm{\varepsilon}=0$ and therefore
receives a larger punishment from the regularization. The solid lines in Fig.
3(c,d) show the evolution of the population mean of the three parameters and
the shaded regions represent their marginal standard deviation. We see that
from generation to generation the standard deviations of the population
decrease, eventually vanishing as the mean value approaches a minimum in the
loss.

Fig. 3(c) shows a case where the algorithm converged towards the left sweet
spot shown in Fig. 3(a,b), settling at the parameter configuration
$\varepsilon_{1}=\varepsilon_{2}=-0.171\,\Delta$ and
$\varepsilon_{S1}=-0.321\,\Delta$. The “target” values, corresponding to the
corresponding sweet spot identified in Ref. [33], are indicated at the right
vertical axis, where
$\varepsilon_{1}^{\ast}=\varepsilon_{2}^{\ast}=-0.151\,\Delta$ and
$\varepsilon_{S1}^{\ast}=-0.319\,\Delta$. In $41$ cases out of $50$, the
algorithm converged towards this sweet spot. In the remaining $9$ cases, the
algorithm converged towards the other sweet spot, an example of which is shown
in Fig. 3(d). Here the final tuning parameters are
$\varepsilon_{1}=\varepsilon_{2}=0.0622\,\Delta$ and
$\varepsilon_{S1}=0.615\,\Delta$, which is again close the true location of
the right sweet spot at
$\varepsilon_{1}^{\ast}=\varepsilon_{2}^{\ast}=0.0785\,\Delta$ and
$\varepsilon_{S1}^{\ast}=0.634\,\Delta$. The bias towards finding the left
sweet spot is due to the initial MvN being centered around
$\bm{\varepsilon}=0$, which is closer to the left sweet spot than the right
one. Convergence in these simulations was typically achieved after $80$–$100$
generations.

Thus we see that the algorithm always converges to tuning parameters that
indeed lie very close to the sweet spots found in Ref. [33], deviating from
the sweet spots by Euclidian distances
$\lVert\bm{\varepsilon}-\bm{\varepsilon^{\ast}}\rVert_{2}$ of $0.028\,\Delta$
and $0.030\,\Delta$ for the left and right sweet spots respectively. As a
result, the left sweet spot estimate found by the algorithm has $\delta
E_{eo}=2.7\times 10^{-3}\,\Delta$ and $M=0.967$ while the right sweet spot
estimate has $\delta E_{eo}=5.7\times 10^{-3}\,\Delta$ and $M=0.977$. In both
cases we thus find high values for the MP and low splittings $\delta E_{eo}$
333Although such residual $\delta E_{eo}$ would still translate into a
relatively short upper-bound time scale for “Majorana manipulation,” it lies
well within the resolution of an actual tunnelling-spectroscopy experiment,
likely making it impossible to resolve such small energy splitting anyway.
Furthermore, the automated tuning found here could serve as the starting point
for a finer search using other methods., and the two-site Kitaev chain results
indicate that the very simple tunneling spectroscopy metric together with the
CMA-ES algorithm can indeed tune on-site potentials $\bm{\varepsilon}$ to
their sweet spots.

###  III.2 Three-site Kitaev chain

To get an idea of how the automatic tuning will fare for longer chains, where
the existence and locations of sweet spots are not known, we also tested it
for the case of $n=3$, i.e., a five-dot array corresponding to a three-site
Kitaev chain. We again run the optimization $50$ times, to obtain statistics
on what solutions the algorithm tends to find. In this case we add an $L_{2}$
term with $\alpha=0.02$ as described by Eq. 5. The larger $\alpha$ is set, the
more biased the algorithm is towards finding the sweet spot that is closest to
the origin $\boldsymbol{\varepsilon}=0$. In an experiment, this parameter
should be tuned empirically, where in order to determine a suitable $\alpha$
one could run the algorithm first without any regularization to find the
typical order of magnitude for the objective function, and based on this, one
can estimate how much regularization is needed to keep the gate voltages
inside a given parameter range.

| $\varepsilon_{1}/\Delta$ | $\varepsilon_{S_{1}}/\Delta$ | $\varepsilon_{2}/\Delta$ | $\varepsilon_{S_{2}}/\Delta$ | $\varepsilon_{3}/\Delta$ | MP | $\delta E_{eo}/\Delta$ | Loss
---|---|---|---|---|---|---|---|---
$\bigtriangledown$ | $-0.10$ | $0.12$ | $0.14$ | $0.12$ | $-0.10$ | $0.998$ | $7.1\cdot 10^{-4}$ | $4.1\cdot 10^{-3}$
$\square$ | $0.19$ | $0.14$ | $-0.31$ | $-0.28$ | $-0.18$ | $0.979$ | $2.5\cdot 10^{-4}$ | $1.1\cdot 10^{-2}$
| $0.07$ | $0.20$ | $-0.40$ | $-0.32$ | $-0.18$ | $0.985$ | $8.7\cdot 10^{-4}$ | $1.2\cdot 10^{-2}$
Table 1: (left) Overview of the three configurations of automatically tuned
on-site potentials found in $50$ simulated experiments on a three-site Kitaev
chain, using the loss function Eq. 5 with $\alpha=0.02$. (right) The
corresponding Majorana polarization, even–odd ground state energy splitting,
and final loss. The number of occurrences of each configuration is shown in
the inset of Fig. 4(a).

The resulting automatically tuned potential configurations are summarized in
Table 1 and a histogram of their occurrences is shown in the inset in Fig.
4(a). All of the simulations converge to configurations with MPs over $0.95$
and even–odd splittings below $10^{-3}\,\Delta$. The exact sweet spot
locations should have exactly zero $\delta E_{eo}$, and while we do see some
finite value for all solutions, their small magnitude indicates that we get
very close to the sweet spots. We see that a large majority of the runs ($42$
out of $50$) ends up at a sweet spot with a MP of 0.998, which is
significantly larger than at the two-site sweet spot. Indeed, the addition of
extra sites is expected to increase the maximal MPs that can be achieved, as
compared to the minimal two-site chain. To illustrate the working of the
algorithm, we show the optimization process for one of these $42$ high-MP
simulations in Fig. 4. In Fig. 4(a) we show the CMA-ES population means and
standard deviations for each parameter across generations and in Fig. 4(b) we
plot the corresponding values for the MP $M$ and energy splitting $\delta
E_{eo}$.

We further note that for the 8 remaining cases the algorithm tuned the system
to non-symmetric parameter configurations with high MP and low $\delta
E_{eo}$, see Table 1, which looks odd at first sight. However, a feature of
the three-site chain is that two pairs of ECT and CAR processes must be
balanced rather than just one as in the minimal chain case. The balanced pairs
may differ from each other, allowing indeed for non-symmetrical configurations
of the tuned on-site potentials along the chain. In principle, these non-
symmetric configurations could thus also correspond to valid sweet spots.

In the three-site chain simulations, convergence was achieved almost always
within $140$–$200$ generations due to the change in the loss function becoming
smaller than the set tolerance, two exceptions being experiments that
terminated as a result of reaching the $200$-generation limit. We finally note
that, as one increases the system size to longer and longer chains (without
drastically increasing $\lambda$), it can be expected that the algorithm will
to a lesser extent find one sweet spot most of the time, and rather discover
several different minima. This is a result of the sample statistics becoming
more and more sparse as the dimensionality of the parameter space becomes
large. In the opposite limit (small dimensionality and many samples
$\lambda$), we expect the algorithm to become more and more deterministic,
approaching a situation where one is essentially checking all solutions within
the search area given by the MvN.

Figure 4: Development of the $100$ first generations for an example run of the
CMA-ES algorithm for the three-site Kitaev chain, where it converges to the
most commonly found parameter configuration. (a) Population mean of the
parameters (solid lines) and their marginal standard deviation (shaded
regions). (b) Majorana polarization $M$ (blue) and even–odd ground state
splitting in the Kitaev chain $\delta E_{oe}$ (red) for each of the first
$100$ generations. (Inset in a) Histogram of the occurrences of the three
different sweet spots discovered in $50$ simulations.

##  IV Conclusions

We have simulated short artificial Kitaev chains based on arrays of quantum
dots, through exact diagonalization of a simple model Hamiltonian, and we
showed that tunneling spectroscopy in the presence of additional sensor dots
at the ends of the chain could be used to automatically tune the system to a
Majorana sweet spot, for which we specifically used the CMA-ES algorithm in
this work. We tested a simple metric based solely on the maximum difference
between even and odd ground state energies as a function of the sensor dot
voltages, suggested by Ref. [48], together with $L_{2}$ regularization on the
gate voltages. For a two-site chain, the algorithm was able to find the known
sweet spots for this system, as identified in Ref. [33]. In the case of a
three-site chain, the algorithm was also capable of reliably tuning to
apparent sweet spots, provided sufficient $L_{2}$ regularization was applied.
We conclude that a simultaneous automatic tuning of all on-site potentials
using tunneling spectroscopy through sensor dots and the CMA-ES algorithm
presents a promising way forward for automatizing the search for Majorana
sweet spots in quantum-dot-based artificial Kitaev chains.

## Acknowledgements

We acknowledge funding via the QuantERA II Programme under Grant Agreement No.
101017733, the work being part of INTFELLES-Project No. 333990, which is
funded by the Research Council of Norway (RCN). This work was further funded
by the Swedish Research Council under Grant Agreement No. 2020-03412, the
Spanish CM Talento Program (Project No. 2022-T1/IND-24070), the Spanish
Ministry of Science, Innovation, and Universities through Grant No.
PID2022-140552NA-I00, the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme under Grant Agreement
No. 856526, Nanolund, and the Horizon Europe Framework Program of the European
Commission through the European Innovation Council Pathfinder Grant No.
101115315, QuKiT. Simulations were performed on resources provided by the NTNU
IDUN/EPIC computing cluster [60], and we thank the NTNU HPC group for their
technical support.

## References

  * Beenakker [2013] C. W. J. Beenakker, Search for Majorana fermions in superconductors, Ann. Rev. Cond. Matt. Phys. 4, 113 (2013).
  * Beenakker [2020] C. W. J. Beenakker, Search for non-Abelian Majorana braiding statistics in superconductors, SciPost Phys. Lect. Notes, 15 (2020).
  * Aasen  _et al._ [2016] D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon, M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flensberg, and J. Alicea, Milestones toward Majorana-based quantum computing, Phys. Rev. X 6, 031016 (2016).
  * Alicea [2012] J. Alicea, New directions in the pursuit of Majorana fermions in solid state systems, Rep. Prog. Phys. 75, 076501 (2012).
  * Alicea  _et al._ [2011] J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Non-Abelian statistics and topological quantum information processing in 1D wire networks, Nat. Phys. 7, 412 (2011).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Kitaev [2001] A. Y. Kitaev, Unpaired Majorana fermions in quantum wires, Phys.-Usp. 44, 131 (2001).
  * Flensberg  _et al._ [2021] K. Flensberg, F. von Oppen, and A. Stern, Engineered platforms for topological superconductivity and Majorana zero modes, Nat. Rev. Mater. 6, 944 (2021).
  * Deng  _et al._ [2016] M. T. Deng, S. Vaitiekėnas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus, Majorana bound state in a coupled quantum-dot hybrid-nanowire system, Science 354, 1557 (2016).
  * Deng  _et al._ [2012] M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff, and H. Q. Xu, Anomalous zero-bias conductance peak in a Nb–InSb nanowire–Nb hybrid device, Nano Lett. 12, 6414 (2012).
  * Finck  _et al._ [2013] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Anomalous modulation of a zero-bias peak in a hybrid nanowire-superconductor device, Phys. Rev. Lett. 110, 126406 (2013).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Signatures of Majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012).
  * Nichele  _et al._ [2017] F. Nichele, A. C. Drachmann, A. M. Whiticar, E. C. O’Farrell, H. J. Suominen, A. Fornieri, T. Wang, G. C. Gardner, C. Thomas, A. T. Hatke, P. Krogstrup, M. J. Manfra, K. Flensberg, and C. M. Marcus, Scaling of Majorana zero-bias conductance peaks, Phys. Rev. Lett. 119, 136803 (2017).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Helical liquids and Majorana bound states in quantum wires, Phys. Rev. Lett. 105, 177002 (2010).
  * Prada  _et al._ [2020] E. Prada, P. San-Jose, M. W. A. de Moor, A. Geresdi, E. J. H. Lee, J. Klinovaja, D. Loss, J. Nygård, R. Aguado, and L. P. Kouwenhoven, From Andreev to Majorana bound states in hybrid superconductor–semiconductor nanowires, Nat. Rev. Phys. 2, 575 (2020).
  * Liu  _et al._ [2012] J. Liu, A. C. Potter, K. T. Law, and P. A. Lee, Zero-bias peaks in the tunneling conductance of spin-orbit-coupled superconducting wires with and without Majorana end-states, Phys. Rev. Lett. 109, 267002 (2012).
  * Pikulin  _et al._ [2012] D. I. Pikulin, J. P. Dahlhaus, M. Wimmer, H. Schomerus, and C. W. J. Beenakker, A zero-voltage conductance peak from weak antilocalization in a Majorana nanowire, New J. Phys. 14, 125011 (2012), publisher: IOP Publishing.
  * Kells  _et al._ [2012] G. Kells, D. Meidan, and P. W. Brouwer, Near-zero-energy end states in topologically trivial spin-orbit coupled superconducting nanowires with a smooth confinement, Phys. Rev. B 86, 100503 (2012).
  * Zhang  _et al._ [2017] H. Zhang, Ö. Gül, S. Conesa-Boj, M. P. Nowak, M. Wimmer, K. Zuo, V. Mourik, F. K. de Vries, J. van Veen, M. W. A. de Moor, J. D. S. Bommer, D. J. van Woerkom, D. Car, S. R. Plissard, E. P. A. M. Bakkers, M. Quintero-Pérez, M. C. Cassidy, S. Koelling, S. Goswami, K. Watanabe, T. Taniguchi, and L. P. Kouwenhoven, Ballistic superconductivity in semiconductor nanowires, Nat. Commun. 8, 16025 (2017).
  * Reeg  _et al._ [2018] C. Reeg, O. Dmytruk, D. Chevallier, D. Loss, and J. Klinovaja, Zero-energy Andreev bound states from quantum dots in proximitized Rashba nanowires, Phys. Rev. B 98, 245407 (2018).
  * Woods  _et al._ [2019] B. D. Woods, J. Chen, S. M. Frolov, and T. D. Stanescu, Zero-energy pinning of topologically trivial bound states in multiband semiconductor-superconductor nanowires, Phys. Rev. B 100, 125407 (2019).
  * Pan and Das Sarma [2020] H. Pan and S. Das Sarma, Physical mechanisms for zero-bias conductance peaks in Majorana nanowires, Phys. Rev. Research 2, 013377 (2020).
  * Ahn  _et al._ [2021] S. Ahn, H. Pan, B. Woods, T. D. Stanescu, and S. Das Sarma, Estimating disorder and its adverse effects in semiconductor Majorana nanowires, Phys. Rev. Materials 5, 124602 (2021).
  * Das Sarma and Pan [2021] S. Das Sarma and H. Pan, Disorder-induced zero-bias peaks in Majorana nanowires, Phys. Rev. B 103, 195158 (2021).
  * Yu  _et al._ [2021] P. Yu, J. Chen, M. Gomanko, G. Badawy, E. P. A. M. Bakkers, K. Zuo, V. Mourik, and S. M. Frolov, Non-Majorana states yield nearly quantized conductance in proximatized nanowires, Nat. Phys. 17, 482 (2021).
  * Hess  _et al._ [2021] R. Hess, H. F. Legg, D. Loss, and J. Klinovaja, Local and nonlocal quantum transport due to Andreev bound states in finite Rashba nanowires with superconducting and normal sections, Phys. Rev. B 104, 075405 (2021).
  * Cayao and Burset [2021] J. Cayao and P. Burset, Confinement-induced zero-bias peaks in conventional superconductor hybrids, Phys. Rev. B 104, 134507 (2021).
  * Hess  _et al._ [2023] R. Hess, H. F. Legg, D. Loss, and J. Klinovaja, Trivial Andreev band mimicking topological bulk gap reopening in the nonlocal conductance of long Rashba nanowires, Phys. Rev. Lett. 130, 207001 (2023).
  * Sau and Sarma [2012] J. D. Sau and S. D. Sarma, Realizing a robust practical Majorana chain in a quantum-dot-superconductor linear array, Nat. Commun. 3, 964 (2012).
  * Leijnse and Flensberg [2012] M. Leijnse and K. Flensberg, Parity qubits and poor man’s Majorana bound states in double quantum dots, Phys. Rev. B 86, 134528 (2012).
  * Miles  _et al._ [2023] S. Miles, D. van Driel, M. Wimmer, and C.-X. Liu, Kitaev chain in an alternating quantum dot-Andreev bound state array, arXiv:2309.15777 (2023).
  * Tsintzis  _et al._ [2022] A. Tsintzis, R. S. Souto, and M. Leijnse, Creating and detecting poor man’s Majorana bound states in interacting quantum dots, Phys. Rev. B 106, L201404 (2022).
  * Tsintzis  _et al._ [2024] A. Tsintzis, R. S. Souto, K. Flensberg, J. Danon, and M. Leijnse, Majorana qubits and non-abelian physics in quantum dot–based minimal kitaev chains, PRX Quantum 5, 010323 (2024).
  * Luna  _et al._ [2024] J. D. T. Luna, A. M. Bozkurt, M. Wimmer, and C.-X. Liu, Flux-tunable Kitaev chain in a quantum dot array, arXiv:2402.07575 (2024).
  * Liu  _et al._ [2023a] C.-X. Liu, H. Pan, F. Setiawan, M. Wimmer, and J. D. Sau, Fusion protocol for Majorana modes in coupled quantum dots, Phys. Rev. B 108, 085437 (2023a).
  * Souto and Aguado [2024] R. S. Souto and R. Aguado, Subgap states in semiconductor-superconductor devices for quantum technologies: Andreev qubits and minimal Majorana chains, arXiv:2404.06592 (2024).
  * Dvir  _et al._ [2023] T. Dvir, G. Wang, N. van Loo, C.-X. Liu, G. P. Mazur, A. Bordin, S. L. D. ten Haaf, J.-Y. Wang, D. van Driel, F. Zatelli, X. Li, F. K. Malinowski, S. Gazibegovic, G. Badawy, E. P. A. M. Bakkers, M. Wimmer, and L. P. Kouwenhoven, Realization of a minimal Kitaev chain in coupled quantum dots, Nature 614, 445 (2023).
  * ten Haaf  _et al._ [2023] S. L. D. ten Haaf, Q. Wang, A. M. Bozkurt, C.-X. Liu, I. Kulesh, P. Kim, D. Xiao, C. Thomas, M. J. Manfra, T. Dvir, M. Wimmer, and S. Goswami, Engineering Majorana bound states in coupled quantum dots in a two-dimensional electron gas, arXiv:2311.03208 (2023).
  * Zatelli  _et al._ [2023] F. Zatelli, D. van Driel, D. Xu, G. Wang, C.-X. Liu, A. Bordin, B. Roovers, G. P. Mazur, N. van Loo, J. C. Wolff, A. M. Bozkurt, G. Badawy, S. Gazibegovic, E. P. A. M. Bakkers, M. Wimmer, L. P. Kouwenhoven, and T. Dvir, Robust poor man’s Majorana zero modes using Yu-Shiba-Rusinov states, arXiv:2311.03193 (2023).
  * Bordin  _et al._ [2024a] A. Bordin, X. Li, D. van Driel, J. C. Wolff, Q. Wang, S. L. D. ten Haaf, G. Wang, N. van Loo, L. P. Kouwenhoven, and T. Dvir, Crossed Andreev reflection and elastic cotunneling in three quantum dots coupled by superconductors, Phys. Rev. Lett. 132, 056602 (2024a).
  * Bordin  _et al._ [2024b] A. Bordin, C.-X. Liu, T. Dvir, F. Zatelli, S. L. D. t. Haaf, D. van Driel, G. Wang, N. van Loo, T. van Caekenberghe, J. C. Wolff, Y. Zhang, G. Badawy, S. Gazibegovic, E. P. A. M. Bakkers, M. Wimmer, L. P. Kouwenhoven, and G. P. Mazur, Signatures of Majorana protection in a three-site Kitaev chain, arXiv:2402.19382 (2024b).
  * Liu  _et al._ [2022] C.-X. Liu, G. Wang, T. Dvir, and M. Wimmer, Tunable superconducting coupling of quantum dots via Andreev bound states in semiconductor-superconductor nanowires, Phys. Rev. Lett. 129, 267701 (2022).
  * Liu  _et al._ [2023b] C.-X. Liu, A. M. Bozkurt, F. Zatelli, S. L. D. t. Haaf, T. Dvir, and M. Wimmer, Enhancing the excitation gap of a quantum-dot-based Kitaev chain, arXiv:2310.09106 (2023b).
  * Thamm and Rosenow [2023] M. Thamm and B. Rosenow, Machine learning optimization of Majorana hybrid nanowires, Phys. Rev. Lett. 130, 116202 (2023).
  * Prada  _et al._ [2017] E. Prada, R. Aguado, and P. San-Jose, Measuring Majorana nonlocality and spin structure with a quantum dot, Phys. Rev. B 96, 085418 (2017).
  * Clarke [2017] D. J. Clarke, Experimentally accessible topological quality factor for wires with zero energy modes, Phys. Rev. B 96, 201109 (2017).
  * Souto  _et al._ [2023] R. S. Souto, A. Tsintzis, M. Leijnse, and J. Danon, Probing Majorana localization in minimal Kitaev chains through a quantum dot, Phys. Rev. Research 5, 043182 (2023).
  * Koch  _et al._ [2023] R. Koch, D. van Driel, A. Bordin, J. L. Lado, and E. Greplova, Adversarial Hamiltonian learning of quantum dots in a minimal Kitaev chain, Phys. Rev. Applied 20, 044081 (2023).
  * Hansen [2016] N. Hansen, The CMA evolution strategy: A tutorial, arXiv:1604.00772 (2016).
  * Hansen and Ostermeier [2001] N. Hansen and A. Ostermeier, Completely derandomized self-adaptation in evolution strategies, Evol. Comput. 9, 159 (2001).
  * Hansen  _et al._ [2003] N. Hansen, S. D. Müller, and P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evol. Comput. 11, 1 (2003).
  * Note [1] Including small Coulomb interactions on the proximitized dots does not change the qualitative results [33].
  * Note [2] We use the cma Python library from pypi.org/project/cma.
  * Kouwenhoven  _et al._ [1997] L. P. Kouwenhoven, C. M. Marcus, P. L. McEuen, S. Tarucha, R. M. Westervelt, and N. S. Wingreen, Electron Transport in Quantum Dots, in  _Mesoscopic Electron Transport_, edited by L. L. Sohn, L. P. Kouwenhoven, and G. Schön (Springer Netherlands, Dordrecht, 1997) pp. 105–214.
  * Sedlmayr and Bena [2015] N. Sedlmayr and C. Bena, Visualizing Majorana bound states in one and two dimensions using the generalized Majorana polarization, Phys. Rev. B 92, 115115 (2015).
  * Sticlet  _et al._ [2012] D. Sticlet, C. Bena, and P. Simon, Spin and Majorana polarization in topological superconducting wires, Phys. Rev. Lett. 108, 096802 (2012).
  * Aksenov  _et al._ [2020] S. V. Aksenov, A. O. Zlotnikov, and M. S. Shustin, Strong Coulomb interactions in the problem of Majorana modes in a wire of the nontrivial topological class BDI, Phys. Rev. B 101, 125431 (2020).
  * Note [3] Although such residual $\delta E_{eo}$ would still translate into a relatively short upper-bound time scale for “Majorana manipulation,” it lies well within the resolution of an actual tunnelling-spectroscopy experiment, likely making it impossible to resolve such small energy splitting anyway. Furthermore, the automated tuning found here could serve as the starting point for a finer search using other methods.
  * Själander  _et al._ [2022] M. Själander, M. Jahre, G. Tufte, and N. Reissmann, EPIC: An energy-efficient, high-performance gpgpu computing research infrastructure, arXiv:1912.05848 (2022).


# On the contribution of the X-ray source to the extended nebular HeII
emission in IZw18

C. Kehrig11affiliationmark:  M.A. Guerrero11affiliationmark:  J. M.
Vílchez11affiliationmark:  G. Ramos-Larios22affiliationmark:
11affiliationmark:  Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta
de la Astronomía s/n, E-18008 Granada, Spain  22affiliationmark:  Instituto de
Astronomía y Meteorología, Dpto. de Física, CUCEI, Universidad de Guadalajara,
Av. Vallarta No. 2602, C.P. 44130, Guadalajara, Jalisco, Mexico

###### Abstract

Nebular Heii emission implies the presence of energetic photons (E$\geq$54
eV). Despite the great deal of effort dedicated to understanding Heii
ionization, its origin has remained mysterious, particularly in metal-
deficient star-forming (SF) galaxies. Unfolding Heii-emitting, metal-poor
starbursts at z$\sim$0 can yield insight into the powerful ionization
processes occurring in the primordial universe. Here we present a new study on
the effects that X-ray sources have on the Heii ionization in the extremely
metal-poor galaxy IZw18 (Z $\sim$3$\%$ Z⊙), whose X-ray emission is dominated
by a single high-mass X-ray binary (HMXB). This study uses optical integral
field spectroscopy, archival Hubble Space Telescope observations, and all of
the X-ray data sets publicly available for IZw18. We investigate the time-
variability of the IZw18 HMXB for the first time; its emission shows small
variations on timescales from days to decades. The best-fit models for the
HMXB X-ray spectra cannot reproduce the observed Heii ionization budget of
IZw18, nor can recent photoionization models that combine the spectra of both
very low metallicity massive stars and the emission from HMXB. We also find
that the IZw18 HMXB and the Heii-emission peak are spatially displaced at a
projected distance of $\simeq$ 200 pc. These results reduce the relevance of
X-ray photons as the dominant Heii ionizing mode in IZw18, which leaves
uncertain what process is responsible for the bulk of its HeII ionization.
This is in line with recent work discarding X-ray binaries as the main source
responsible for Heii ionization in SF galaxies.

###### Subject headings:

galaxies: dwarf — galaxies: individual (IZw18) — galaxies: ISM — X-rays:
binaries

##  1\. Introduction

The presence of nebular Heii emission points to the existence of a hard
radiation field as photon energies $\geq$ 54 eV are necessary to ionize He+.
Highly ionized systems, like Heii emitters, are expected to be more frequent
at high-z (e.g., Cassata, et al., 2013; Mainali, et al., 2018). At the same
time, observations of local, Heii-emitting star-forming (SF) galaxies suggest
that the lower the metallicity (Z), the larger the nebular Heii line
intensities (e.g., Guseva, et al., 2000; Senchyna, et al., 2017). This is in
line with the harder ionizing stellar spectra predicted at the lower
metallicities typical in the early universe (e.g., Tumlinson & Shull, 2000);
the first, hot metal-free Population III stars produce strong nebular Heii
lines (e.g., Yoon, et al., 2012; Visbal, et al., 2015). Still, despite much
observational and theoretical effort, including relevant advances in stellar
modeling, the source of Heii ionization remains puzzling, especially in metal-
poor galaxies (e.g., Eldridge, et al., 2017; Götberg, et al., 2018; Kehrig, et
al., 2018; Stanway & Eldridge, 2019; Kubátová, et al., 2019).

There is evidence for enhanced total X-ray luminosities per unit of star
formation rate ($L_{\rm X}$/SFR) and high-mass X-ray binary (HMXB) populations
in extremely metal-poor galaxies (Z/Z⊙ $<$ 0.1) compared to those of $\sim$Z⊙
(e.g., Prestwich, et al., 2013; Brorby, et al., 2016; Lebouteiller et al.,
2017; Schaerer, et al., 2019; Ponnada, et al., 2020). Furthermore, examples of
nearby Heii nebulae associated with HMXBs that appear to act as the only
ionizing source are reported (e.g., Garnett et al., 1991; Kaaret, et al.,
2004; Gutiérrez & Moon, 2014). Thus HMXBs have been proposed as one of the
leading mechanisms to ionize He+ in SF systems. However, it has been recently
questioned whether they can provide the required He+ ionizing flux (e.g.,
Kehrig, et al., 2018; Senchyna, et al., 2020; Saxena, et al., 2020).
Variability of the X-ray source could alleviate this conflict, but so far only
the extensive study of the puzzling Heii$\lambda$4686-emitting nebula N44C in
the LMC considered the effects of an on-off switching X-ray source in time
periods of years (see Nazé, et al., 2003).

IZw18 is among the most metal-poor (12+log(O/H)= 7.11 $\sim$ 3$\%$ solar111
Assuming a solar metallicity 12+log(O/H)⊙ = 8.69 (Asplund et al. 2009); e.g.,
Kehrig et al. 2016, K16) Heii-emitting SF galaxies known at z $\sim$ 0\. It
thus can be considered a nearby analog of distant metal-poor Heii emitters,
providing a unique laboratory to study the high-ionization phenomenon expected
to be common in the reionization era. Based on integral field spectroscopy
(IFS) data, Kehrig, et al. (2015, K15) unveiled the entire nebular
Heii$\lambda$4686-emitting region in IZw18 with a diameter of $\approx$5′′
($\approx$ 440 pc at the distance of 18.2 Mpc; Aloisi et al. 2007). K15 found
that only peculiar hot (nearly) metal-free massive stars could explain the
Heii ionization in IZw18 (see also Kehrig, et al., 2018, K18). The main goal
of this Letter is to investigate in detail the spatial correlation between the
Heii emission in IZw18 and its single HMXB (Thuan et al., 2004, T04), the
variability of the latter on short (days) and long (decades) timescales and
its expected He+ ionizing photon production to assess the contribution of the
X-ray emission to the nebular Heii$\lambda$4686 ionization.

Table 1X-ray Observations of IZw18

X-ray Observatory | Instrument | Data set | Date | Net Exposure Time | Count Ratea
---|---|---|---|---|---
|  |  |  | (ks) | (ks-1)
ROSAT | PSPC | rp600165n00 | 1992 04 30 | 16.2 | 7.2$\pm$0.7
Chandra | ACIS-I | 600108 | 2000 02 08 | 31.0 | 12.3$\pm$0.6
XMM-Newton | EPIC-pn | 0112520101 | 2002 04 10 | 24.2 | 90.4$\pm$2.0
| EPIC-MOS1 |  |  | 31.1 | 28.3$\pm$1.0
| EPIC-MOS2 |  |  | 30.9 | 28.4$\pm$1.0
XMM-Newton | EPIC-pn | 0112520201 | 2002 04 16 | 5.4 | 76.9$\pm$4.0
| EPIC-MOS1 |  |  | 10.4 | 26.8$\pm$1.7
| EPIC-MOS2 |  |  | 10.4 | 29.2$\pm$1.7
Swift | XRT | 00040765001 | 2011 03 25 | 2.0 | 2.1$\pm$1.1
Swift | XRT | 00040765002 | 2011 03 29 | 0.9 | 6.6$\pm$2.8
Swift | XRT | 00040765003 | 2012 03 23 | 0.4 | 9.2$\pm$5.0
Swift | XRT | 00040765004 | 2012 03 26 | 0.3 | 8.9$\pm$5.6
Swift | XRT | 00040765005 | 2012 08 29 | 2.6 | 6.7$\pm$1.6

(a) Count rates are computed within the nominal energy bands of each X-ray
observatory: ROSAT/PSPC (0.1-2.4 keV), Chandra/ACIS-I (0.3-10.0 keV), XMM-
Newton/EPIC (0.3-10.0 keV), and Swift/XRT (0.2-10.0 keV).

##  2\. Observations

This study combines archival X-ray and Hubble Space Telescope (HST)
observations with IFS data.

The IFS observations of IZw18, obtained on 2012 December, use the Potsdam
Multi-Aperture Spectrophotometer (PMAS; Roth et al., 2005) at the Calar Alto
Observatory (Spain). PMAS provides a field of view (FOV) of
$16^{\prime\prime}\times 16^{\prime\prime}$ with a sampling of
$1^{\prime\prime}\times 1^{\prime\prime}$ over $\sim$3640–7200 Å (see K15,K16
for details). The data have been used to obtain a continuum-subtracted
spectral map of the nebular Heii$\lambda$4686 emission line in IZw18.

IZw18 has been targeted by most modern X-ray observatories, including ROSAT,
Chandra, XMM-Newton and Swift. Table 1 lists the details of these X-ray
observations. The Chandra and XMM-Newton data were reprocessed using the
specific reduction packages Chandra Interactive Analysis of Observations (CIAO
v4.9, Fruscione et al., 2006) and XMM-Newton Science Analysis System (SAS
v17.0, Gabriel et al., 2004), respectively. The ROSAT and Swiftdatasets were
processed using generic HEASARC reduction packages. IZw18 is detected as a
point source in all these observations except in the Chandra observations,
with the highest spatial resolution ($\approx$0$\farcs$5) among the available
X-ray data, which suggest faint extended emission toward the south (T04).
Table 1 provides the count rates of these detections, which are below a
3$\sigma$ threshold in the shortest Swiftobservations.

We also use archival HST images of IZw18. The location of the hot, young
massive stars is probed using the HST STIS image obtained through the far-UV
filter F25SRF2 with $\lambda_{c}$=1457 Å (PI:T.M. Brown, prop. 9054, Brown et
al., 2002). The HST WFPC2/F658N image (PI:E. Skillman, prop. 6536, de Mello et
al., 1998) outlines the ionized gas emission in the H$\alpha$ line.

##  3\. Data Analysis

###  3.1. UV, optical and X-ray Imaging

Previous analyses of the spatial properties of the X-ray emission in IZw18
were presented by T04 using the same Chandra data set used here. They
suggested the presence of an HMXB spatially coincident with the bright
northwest (NW) SF knot of IZw18 and a fainter diffuse emission toward the
south. To investigate the spatial distribution of the X-ray emission from
IZw18, we have produced images in the soft 0.3-1.2 keV, medium 1.2-2.5 keV,
hard 2.5-8.0 keV, and broad 0.3-8.0 keV X-ray bands using the ACIS subpixel
event resolution capability, with subpixel size 1/4 their natural size of
0$\farcs$5\. Images were subsequently convolved with subpixel PSF images of
the same energy and adaptively smoothed with a Gaussian kernel size between
0$\farcs$5-2$\farcs$0.

To spatially correlate the different components of IZw18, we compare in detail
the distribution of the emission of X-ray with that of the hot young massive
stars and ionized gas from HST images, and with the PMAS Heii$\lambda$4686
map. The HST and Chandra images were registered using USNO stars in their FoV.
As noted by T04, the astrometric accuracy of the Chandra observations is
$\approx$0$\farcs$5 due to the limited number of optical counterparts of X-ray
sources. Similarly, the HST images and PMAS Heii$\lambda$4686 map were
registered using the HST image and PMAS map in the H$\alpha$ emission line.
Their relative registration accuracy is estimated to be within
$\approx$0$\farcs$5.

Figure 1.— Multiwavelength images of IZw18 and Chandra ACIS-I X-ray contours.
(Left panel) Color-composite Chandra ACIS-I X-ray picture in the soft 0.3-1.2
keV (red), medium 1.2-2.5 keV (green), and hard 2.5-8.0 keV (blue) energy
bands. (Center panel) Color-composite optical HST WFPC2/F658N H$\alpha$
emission line (red), HST STIS/FUV/F25SRF2 continuum at 1457 Å (green), and
X-ray Chandra ACIS-I 0.3-8.0 keV (blue) picture. (Right panel) PMAS map of
nebular Heii$\lambda$4686 emission. In this map, the pixel scale has been
regrid to 0$\farcs$1 applying a 2D linear interpolation, and then smoothed at
the original pixel scale using a $\sigma$=1′′ Gaussian kernel; intensity is
presented in the square root scale. The orientation, FoV, and scale of the
three panels match each other. The Heii emission peak in the right panel is
marked by a ’+’ symbol in the left and central panels. The white contours
correspond to the PSF of the X-ray point source, and display the intensity of
the 0.3-8.0 keV X-ray emission using a logarithmic scale from a 3$\sigma$
level to the intensity peak; the red contour shows the 3$\sigma$ level of the
soft 0.3-1.2 keV emission. The nominal FWHM of the Chandra ACIS PSF,
$\simeq$0$\farcs$5, is in between the first and second highest white contours,
whereas the Heii emission peak is observed well outside the latter.

Figure 1 (left) presents a color-composite picture of the soft, medium, and
hard X-ray band images. This image confirms the presence of a hard X-ray point
source at the position of the NW knot (T04) and reveals that the diffuse
emission toward the south is redder, i.e., the diffuse emission is softer than
the point source. Figure 1 (center) shows that the HMXB emission (blue) is
coincident with the young massive stellar cluster (green) centered in the NW
component of H$\alpha$ emission (red) of IZw18. The diffuse X-ray emission
extends south, along some holes in the H$\alpha$ emission. Figure 1 (right)
displays the extended nebular Heii$\lambda 4686$-emitting gas showing, for the
first time, that the Heii$\lambda 4686$ emission peak is displaced
$\approx$2$\farcs$5 ($\approx$200 pc projected at the distance of 18.2 Mpc)
from the south of the HMXB in IZw18. The Heii emission overlaps part of the
much more extended H$\alpha$ emission, whose peaks appear noncoincident (see
also Figure 2 (middle) panels from K15).

###  3.2. X-ray Variability

We used the different X-ray datasets listed in Table 1 to extract X-ray
spectra utilizing apertures encompassing the X-ray emission associated with
the NW SF knot. Background spectra were extracted from nearby source-free
regions. The source and background spectra were then processed to obtain
calibration files. The background-subtracted ROSAT/PSPC, Chandra/ACIS-I, and
XMM-Newton/EPIC X-ray spectra of IZw18 are shown in Figure 2. The Swiftspectra
are not shown given their limited count number.

Previous analyses of the X-ray spectrum of the IZw18 HMXB considered plasma
emission and power-law models for ROSAT (Martin, 1996) and Chandra data (T04),
and power-law and thin accretion disk around a Kerr black hole (BH) models for
the XMM-Newton data obtained on 2002 April 10 (2002-04-10-XMM-Newton; Kaaret &
Feng, 2013). The latter model is the most specific for HMXBs, but too complex
for the quality of some of the available datasets. To preserve the consistency
of the investigation of the HMXB X-ray variability, a power-law model will be
used for all datasets. For the fit, a fixed Galactic extinction along the line
of sight of IZw18 consistent with a hydrogen column density ($N_{\rm H}$) of
2.9$\times$1020 cm-2 (HI4PI Collaboration, 2016) is assumed. An additional
absorption component with the same chemical abundance as IZw18 is adopted to
describe the local absorption.

The results of our spectral modeling are summarized in Table 2. We first
started modeling the 2002-04-10-XMM-Newton spectra, which have the largest
count number. The power-law model (black, blue, and green histograms in Figure
2) results in an acceptable fitting ($\chi^{2}$/DoF=1.33). The best-fit power-
law photon index $\Gamma$=1.95 is typical of synchrotron emission, whereas the
$N_{\rm H}$ of the local component (=3.6$\times$1021 cm-2) is consistent with
that derived by Lelli et al. (2012) at the location of the X-ray source. We
then fitted the XMM-Newton/EPIC spectra extracted from the observations on
2002 April 16, using a similar absorbed power-law model. The spectra have
lower quality than those obtained just 6 days before given their shorter net
exposure times, but the parameters of the best-fit models are consistent
within their uncertainties with those obtained for the 2002-04-10-XMM-Newton
observations (see Table 2).

Figure 2.—  Background-subtracted ROSAT (top left), Chandra (top right), and
XMM-Newton (bottom) X-ray spectra of IZw18 (points) overplotted by best-fit
models (histograms). The error bars of each data point correspond to 1$\sigma$
statistical uncertainties. The best-quality XMM-Newton spectra are fitted
using an absorbed power law (black for EPIC-pn, blue for EPIC-MOS1, and green
for EPIC-MOS2). An absorbed irradiated disk model is also shown for the EPIC-
pn spectrum (red). The ROSAT and Chandra spectra are fitted using an absorbed
power law models (black), with the scaled best-fit power law model of the
2002-04-10-XMM-Newton observations also overplotted (purple). Lower panels
show the residuals of the best-fit models where the data points’ colors
correspond to those of the best-fit models.  Table 2X-ray Spectral Fittings
and Luminosities of IZw18

X-ray Observatory | Date | NH | $\Gamma$ |  $\chi^{2}$/DoF | $f_{\rm 0.3-2.4keV}$ | $f_{\rm 0.5-8.0keV}$ | $L_{\rm 0.3-2.4keV}$ | $L_{\rm 0.5-8.0keV}$
---|---|---|---|---|---|---|---|---
and Instrument |  | (1021 cm-2) |  |  | (10-14 erg cm-2 s-1) | (1038 erg cm s-1)
ROSAT PSPCa | 1992 04 30 | 16${}^{+34}_{-5}$ | 4.2${}^{+4.3}_{-2.0}$ | 12.7/14 = 0.91 | 3.8 | 6.5 | 400${}^{+20}_{-150}$ | 127${}^{+6}_{-60}$
|  | 3.6 | 1.95 | 16.5/16 = 1.03 | 7.5 | 8.9 | 30.0${}^{+4.2}_{-3.3}$ | 35.5${}^{3.8}_{-4.4}$
Chandra ACIS-I | 2000 02 08 | 0.5${}^{+1.3}_{-0.5}$ | 1.96${}^{+0.28}_{-0.21}$ | 15.2/20 = 0.76 | 3.7 | 6.3 | 19.4${}^{+0.8}_{-1.5}$ | 26.9$\pm$2.1
XMM-Newton EPIC | 2002 04 10 | 3.6$\pm$0.6 | 1.95$\pm$0.07 | 158.7/122 = 1.33 | 13.5 | 28.0 | 94.3$\pm$1.5 | 131.2${}^{+1.9}_{-4.0}$
XMM-Newton EPIC | 2002 04 16 | 3.2${}^{+1.1}_{-0.9}$ | 1.84$\pm$0.13 | 33.6/38 = 0.88 | 12.1 | 26.8 | 48.4$\pm$2.1 | 107.2$\pm$4.4
Swift/XRT | 2011 03 25 | 3.6 | 1.95 | $\dots$ | 4.5 | 9.2 | 31$\pm$17 | 44$\pm$23
Swift/XRT | 2011 03 29 | 3.6 | 1.95 | $\dots$ | 14 | 29 | 98$\pm$42 | 135$\pm$60
Swift/XRT | 2012 03 23 | 3.6 | 1.95 | $\dots$ | 20 | 41 | 135$\pm$70 | 190$\pm$100
Swift/XRT | 2012 03 26 | 3.6 | 1.95 | $\dots$ | 19 | 39 | 130$\pm$80 | 180$\pm$110
Swift/XRT | 2012 08 29 | 3.6 | 1.95 | $\dots$ | 14 | 30 | 98$\pm$23 | 138$\pm$33

(a) First and second rows show the results of ROSAT spectral modeling based on
a free power-law fit and using the best-fit model of the 2002-04-10-XMM-Newton
spectrum, respectively.

Figure 3.— Long-term (top) and short-term (bottom) X-ray variability of the
X-ray emission from IZw18. The intrinsic $L_{\rm X,0.5-8~{}keV}$ and
corresponding error bars are derived from the models detailed in Table 2. The
ROSAT red data point corresponds to a free power-law fit, whereas the black
data point uses the best-fit parameters of the fit to the 2002-04-10-XMM-
Newton spectrum. The bottom panels present expanded light-curves of the
Chandra and XMM-Newton observations. The dashed horizontal lines in these
panels correspond to the average X-ray luminosities during the observations.

Next we proceed to fit the ROSAT/PSPC and Chandra/ACIS-I spectra with this
same absorbed power-law model. The parameters of the best-fit models (black
histograms in the top panels of Figure 2) are notably different from those
obtained for the XMM-Newton/EPIC spectra. As for the ROSAT/PSPC spectrum, the
data quality is not adequate to reliably constrain the model parameters (see
their 1$\sigma$ uncertainties in Tab. 2) and indeed the scaled best-fit model
of the 2002-04-10-XMM-Newton spectra (purple histogram in the top left panel
of Figure 2) is consistent with the ROSAT/PSPC spectrum. On the other hand,
the Chandra ACIS-I spectrum is much softer than the scaled best-fit model of
the 2002-04-10-XMM-Newton spectra (purple histogram in the top right panel of
Figure 2), i.e., the internal extinction during the Chandra observation was
much lower than that during the XMM-Newton observations. The Chandra spectrum
also implies the lowest $L_{\rm X}$. Finally, we have adopted the best-fit
model of the 2002-04-10-XMM-Newton spectra to derive the $L_{\rm X}$ from the
count rate of the Swift observations. These are listed in Table 2, where the
uncertainties are derived from the 1$\sigma$ error bar of the count rate.

The long-term evolution of $L_{\rm X}$ in the energy band 0.5-8.0 keV ($L_{\rm
X,0.5-8~{}keV}$; see Table 2) is shown in the upper panel of Figure 3.
Apparently, the IZw18 HMXB switches from a low state with $L_{\rm
X,0.5-8~{}keV}$ a few times 1039 erg s-1 to an upper state with $L_{\rm
X,0.5-8~{}keV}$ of (1–2)$\times$1040 erg s-1. These high $L_{\rm X}$ and
spectral changes are consistent with those of ultraluminous X-ray sources, as
discussed by Kaaret & Feng (2013), although T04 argued that HMXBs in low
metallicity galaxies tend to have higher $L_{\rm X}$. The transition from one
state to the other occurs in timescales shorter than 2-3 yr, as suggested by
the increase in luminosity from the Chandra to the XMM-Newton observations,
and from the different Swift observations, although the latter are affected by
large uncertainties in the luminosity estimates.

The short-term variability is shown in the lower panels of Figure 3.
Lightcurves in the energy band 0.5-8.0 keV with a bin size of 4.0 and 1.2 ks
have been extracted from the Chandra/ACIS-I and XMM-Newton/EPIC-pn
observations, respectively. The count rates and their 1$\sigma$ uncertainties
have been converted to X-ray luminosities using the average $L_{\rm
X,0.5-8~{}keV}$ at each epoch given in the last column of Table 2. For each
data set, the X-ray emission is mostly consistent with its mean value, which
may be suggesting small fluctuations. In this sense, the X-ray level of the
two XMM-Newton observations, obtained $\sim$6 days apart, are mostly
consistent with each other.

This leads us to conclude that the IZw18 HMXB has a relatively constant
emission in short (hours to days) timescales, with variations of factors of
only a few in longer (a few years) timescales. There are hints that the larger
X-ray fluxes are associated with an increase of $N_{\rm H}$ toward this
source, as typically seen in X-ray sources fed by accretion such as HMXBs.

##  4\. Discussion

The nebular Heii$\lambda$4686 emission observed in some SF galaxies has often
been attributed to X-ray emission from HMXBs (e.g. Thuan & Izotov, 2005;
Schaerer, et al., 2019). This hypothesis has been critically tested by
Senchyna, et al. (2020, S20) using photoionization grids accounting for the
combined contribution of HMXBs and young massive stars for metallicities from
Z=0.02 to Z=0.001 ($\simeq$5$\%$ solar), and provide predictions for the
Heii/H$\beta$ ratio that correspond to a range of values of the $L_{\rm
X,0.5-8~{}keV}$/SFR ratio. We estimate a SFR of 0.14 M⊙ yr-1 for IZw18222The
IZw18 integrated $F(H\alpha)$=4.36$\times$10-13erg s-1 cm-2 (K16) was
translated to $L(H\alpha)$ assuming a distance of 18.2 Mpc. based on its
$L(H\alpha)$=1.73$\times$1040erg s-1(K16) and the widely used SFR-L(H$\alpha$)
relation from Kennicutt (1998). Therefore, taking the $L_{\rm X,0.5-8~{}keV}$
value in Table 2 derived from the 2002-04-10-XMM-Newton data set, this implies
an observed $L_{\rm X,0.5-8~{}keV}$/SFR of 9$\times$1040 erg s-1/(M⊙ yr-1).
According to S20’s predictions, this $L_{\rm X,0.5-8~{}keV}$/SFR value results
in Heii/H$\beta$ between 0.0004 to 0.0025, regardless of the metallicity
assumed in the models, a value well below the observed Heii/H$\beta$ of
0.023$\pm$0.004 in IZw18 (K16). To reinforce this result, we have computed the
total Heii-ionizing photon flux $Q$(Heii) from Kerr BH and irradiated disk
models of the X-ray emission detected in the 2002-04-10-XMM-Newton
observation. The thin accretion disk around a Kerr BH model including a
Comptonization component provides an extrapolation to UV wavelengths $<$ 228 Å
similar to that of the multicolor disk models used by S20. The best fit for
this model has very similar parameters to those reported by Kaaret & Feng
(2013), with an acceptable fitting quality ($\chi^{2}$/DoF=0.90). The
intrinsic (unabsorbed) $\log Q$(Heii) for this model is $\approx$48.4 photon
s-1, i.e.,$\approx$50 times lower than that of $\log Q$(Heii)=50.1 photon s-1
reported by K15. Meanwhile, the emission model of a disk irradiated by the
Compton tail (Gierliński et al., 2009) is of particular interest, as the
physical processes of disk irradiation provide a meaningful description of the
expected UV spectrum to constrain more reliably $Q$(Heii) (Berghea & Dudik,
2012). The best fit for this model has a slightly better fit quality
($\chi^{2}$/DoF=0.93) with a local $N_{\rm H}$=1.5$\pm$0.8$\times$1021 cm-2,
inner disk temperature 0.55$\pm$0.25 keV, photon index 2.9$\pm$0.6, and
irradiated flux fraction of 2.3$\%$ for the (other parameters were fixed at
typical values according to; Gierliński et al., 2009). The resulting intrinsic
(unabsorbed) $\log Q$(Heii) is $\approx$49.1 photon s-1, which is still 10
times lower than observed. It is worthwhile noting that the $L_{\rm X}$
derived from these models is very similar to that derived from a power-law
model.

Furthermore, the HMXB is coincident with the NW stellar cluster of IZw18,
whereas the nebular Heii emission (diameter $\approx$440 pc) shows a shell-
like structure extending south of the HMXB; a $\approx$2.′′5 (200 pc) offset
is observed between the Heii emission peak and the HMXB (see Section 3.1).
Interestingly, the diffuse X-ray emission to the south of the HMXB delineates
the southernmost edge of the Heii emission. A spectral analysis of this
diffuse emission is not feasible, but its photon median energy, $\simeq$0.7
keV, is consistent with that found in superbubbles (e.g., DEM L152 in the LMC;
Jaskot et al., 2011). The “true” spatial distribution of the diffuse X-ray
emission might certainly fill the cavity of ionized gas suggested by its Heii
and H$\alpha$ emission, but the bright X-ray emission from the HMXB prevents
the detection of its northern rim. Assuming that the X-ray-emitting plasma
responsible of the diffuse emission had similar temperature as that of the
typical superbubble DEM L152 (0.37 keV; Jaskot et al., 2011), the count rate
from the diffuse emission of 0.39$\pm$0.11 counts ks-1 ($\sim$30 times smaller
than that of the HMXB) implies $L_{\rm X,0.5-8~{}keV}$ $\sim$1038 erg s-1.
This $L_{\rm X}$ is consistent, at the high end, with that of LMC superbubbles
(Dunne et al., 2001). The $Q$(Heii) of this diffuse soft X-ray emission is at
least 4 orders of magnitude lower than that of the hard X-ray point source.

Contrary to our results, previous studies claim that emission from the hard
X-ray point source in IZw18 may be significant for its Heii ionization (e.g.,
Lebouteiller et al., 2017; Heap et al., 2019; Schaerer, et al., 2019).
Lebouteiller et al. (2017) and Heap et al. (2019), however, do not consider
the nebular Heii budget in their modeling, which is key to studying the origin
of the bulk of the Heii ionizing photons (see, e.g., K15; K18; Kehrig et al.,
2020). On the other hand, Schaerer, et al. (2019) assume that the nebular Heii
emission in IZw18 is fully ionized by its HMXB, this way overlooking the
energy-dependent cross section of He+ which must be considered for the
estimation of the intrinsic X-ray ionizing power (see above). Moreover, these
works do not take into consideration the distinct spatial locations of the
X-ray and the Heii emissions in IZw18 unveiled here.

##  5\. Summary and concluding remarks

Cosmic dawn marks a phase transition of the universe during which Population
III-stars and their hosts put an end to the dark ages. The nature of these
first sources is highly unconstrained. The hard spectra of Population III and
nearly Z=0 stars can produce intense Heii emission. However, the origin of
nebular Heii emission in metal-poor SF galaxies, near and far, is still
unclear (e.g., K18; Saxena, et al., 2020).

IZw18 is a very-metal-poor (Z $\sim$ 3$\%$ Z⊙), Heii-emitting SF system, which
is a unique local laboratory for probing the conditions dominating in distant
Heii emitters. Here we combine optical IFS observations with archival HST and
X-ray data of IZw18 to obtain a new view on the effects of X-ray photons on
its nebular Heii-ionization.

We confirm the previously reported HMXB which dominates the IZw18 X-ray
emission, although a softer diffuse emission is also present. The observation
period covered by the X-ray datasets allows us to study the X-ray variability
of the IZw18 HMXB for the first time. Its $L_{\rm X,0.5-8~{}keV}$ is
unchanging in hours-days while it can vary by factors of a few in time-scales
of years. There is thus no observational evidence for episodes of highly
enhanced X-ray emission that, eventually, could power the Heii emission in
IZw18.

We find that the best-fit Kerr BH and irradiated disk models yield,
respectively, $Q$(Heii) $\sim$ 50 and 10 times lower than the observed value.
Using photoionzation grids incorporating both stellar spectra and HMXB models
(S20), we checked that our derived $log(L_{\rm X,0.5-8~{}keV}/SFR$) $\sim$
40.9 corresponds to Heii/H$\beta$ values $\sim$ 9–58 times lower than the
integrated value of Heii/H$\beta$ measured in IZw18 (K16). Finally, we report
for the first time a substantial spatial separation ($\approx$200 pc) between
the HMXB and Heii emission peak.

Although a direct constraint to the $Q$(Heii) produced by an HMXB in extremely
metal-poor galaxies, such as IZw18, has not yet been reported, and one cannot
rule out the contribution of X-rays to the He+ ionization there, all the
points above argue against the HMXB being the main source responsible for the
bulk of the nebular HeII emission in IZw18. Recent findings obtained by
stacking the X-ray emission of Heii-emitting SF galaxies at z $\sim$ 3,
discard X-ray binaries as the dominant producers of He+ ionizing photons in
distant SF galaxies (Saxena, et al., 2020). Along with the detailed study of
the temporal variations and high-resolution mapping of the X-ray emission in
IZw18 presented here, these studies highlight that the main source of Heii
ionization is still puzzling both in nearby and distant SF systems (see also
K15; K18).

## Acknowledgement

We thank the referees for helpful comments and suggestions. This work has been
partially funded by projects AYA2016-79724-C4-4-P and PID2019-107408GB-C44
from the Spanish PNAYA. C.K., M.A.G., and J.V.M. acknowledge financial support
from the State Agency for Research of the Spanish MCIU through the “Center of
Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía
(SEV-2017-0709). M.A.G. acknowledges support from MCIU grant
PGC2018-102184-B-I00 co-funded with FEDER funds. G.R.-L. acknowledges support
from CONACyT (grant 263373) and PRODEP (México). This research has made use of
data and/or software provided by the High Energy Astrophysics Science Archive
Research Center (HEASARC), which is a service of the Astrophysics Science
Division at NASA/GSFC.

## References

  * Aloisi et al. (2007) Aloisi A., et al., 2007, ApJ, 667, L151
  * Asplund et al. (2009) Asplund M., et al., 2009, ARA&A, 47, 481
  * Berghea & Dudik (2012) Berghea, C. T. & Dudik, R. P. 2012, ApJ, 751, 104
  * Brorby, et al. (2016) Brorby M., et al., 2016, MNRAS, 457, 4081
  * Brown et al. (2002) Brown T. M., et al., 2002, ApJ, 579, L75
  * Cassata, et al. (2013) Cassata P., et al., 2013, A&A, 556, A68
  * de Mello et al. (1998) de Mello D. F., et al., 1998, ApJ, 507, 199
  * Dunne et al. (2001) Dunne, B. C., Points, S. D., & Chu, Y.-H. 2001, ApJS, 136, 119
  * Eldridge, et al. (2017) Eldridge J. J., et al., 2017, PASA, 34, e058
  * Fruscione et al. (2006) Fruscione A., et al. 2006, Proc. SPIE, 62701V
  * Gabriel et al. (2004) Gabriel, C., Denby, M., Fyfe, D. J., et al. 2004, Astronomical Data Analysis Software and Systems (ADASS) XIII, 759
  * Garnett et al. (1991) Garnett D. R., et al., 1991, ApJ, 373, 458
  * Gierliński et al. (2009) Gierliński, M., Done, C., & Page, K. 2009, MNRAS, 392, 1106
  * Götberg, et al. (2018) Götberg Y., et al., 2018, A&A, 615, A78
  * Guseva, et al. (2000) Guseva N. G., Izotov Y. I., Thuan T. X., 2000, ApJ, 531, 776
  * Gutiérrez & Moon (2014) Gutiérrez C. M., Moon D.-S., 2014, ApJL, 797, L7
  * HI4PI Collaboration (2016) HI4PI Collaboration, Ben Bekhti, N., Floer, L., et al. 2016, A&A, 594, A116
  * Heap et al. (2019) Heap S. R., Hubeny I., Lanz T. M., 2019, ASPC, 519, 267
  * Jaskot et al. (2011) Jaskot A. E., et al. 2011, ApJ, 729, 28
  * Kaaret, et al. (2004) Kaaret P., Ward M. J., Zezas A., 2004, MNRAS, 351, L83
  * Kaaret & Feng (2013) Kaaret, P., & Feng, H. 2013, ApJ, 770, 20
  * Kehrig, et al. (2015) Kehrig C., et al., 2015, ApJL, 801, L28
  * Kehrig, et al. (2016) Kehrig C., et al., 2016, MNRAS, 459, 2992
  * Kehrig, et al. (2018) Kehrig C., et al., 2018, MNRAS, 480, 1081
  * Kehrig et al. (2020) Kehrig C., et al., 2020, MNRAS, 498, 1638
  * Kennicutt (1998) Kennicutt R. C., 1998, ApJ, 498, 541
  * Kubátová, et al. (2019) Kubátová B., et al., 2019, A&A, 623, A8
  * Lebouteiller et al. (2017) Lebouteiller V., et al., 2017, A&A, 602, A45
  * Lelli et al. (2012) Lelli F., et al. 2012, A&A, 537, A72
  * Mainali, et al. (2018) Mainali R., et al., 2018, MNRAS, 479, 1180
  * Martin (1996) Martin, C. L. 1996, ApJ, 465, 680
  * Nazé, et al. (2003) Nazé Y., et al., 2003, A&A, 408, 171
  * Ponnada, et al. (2020) Ponnada S., Brorby M., Kaaret P., 2020, MNRAS, 491, 3606
  * Prestwich, et al. (2013) Prestwich A. H., et al., 2013, ApJ, 769, 92
  * Roth et al. (2005) Roth M. M., et al., 2005, PASP, 117, 620
  * Saxena, et al. (2020) Saxena A., et al., 2020, MNRAS, arXiv:2003.13800
  * Schaerer, et al. (2019) Schaerer D., Fragos T., Izotov Y. I., 2019, A&A, 622, L10
  * Senchyna, et al. (2017) Senchyna P., et al., 2017, MNRAS, 472, 2608
  * Senchyna, et al. (2020) Senchyna P., et al., 2020, MNRAS, 494, 941
  * Stanway & Eldridge (2019) Stanway E. R., Eldridge J. J., 2019, A&A, 621, A105
  * Thuan et al. (2004) Thuan T. X., et al., 2004, ApJ, 606, 213
  * Thuan & Izotov (2005) Thuan T. X., Izotov Y. I., 2005, ApJS, 161, 240
  * Tumlinson & Shull (2000) Tumlinson J., Shull J. M., 2000, ApJ, 528, L65
  * Visbal, et al. (2015) Visbal E., Haiman Z., Bryan G. L., 2015, MNRAS, 450, 2506
  * Yoon, et al. (2012) Yoon S.-C., Dierks A., Langer N., 2012, A&A, 542, A113