text:
\section{Introduction}
\label{sec:intro}

The top quark is the heaviest known elementary particle and the only quark that decays before hadronisation, which allows direct access to its properties in production and decay. Measurements of top-quark properties, predicted by the Standard Model (SM), provide important input to test theoretical calculations and have the potential to reveal deviations from the SM predictions. One of the relevant properties is related to the slight difference between the rapidity distributions of top quarks and top antiquarks produced in pairs (\ttbar). This asymmetry, referred to as charge asymmetry, is defined in proton--proton collisions as follows~\cite{Aguilar-Saavedra:2014vta, Aguilar-Saavedra:2014kpa, Jung:2011zv, Diener:2009ee}:
\begin{equation*}
\Ac = \frac{N(|y_{t}|>|y_{\bar{t}}|) - N(|y_{t}|<|y_{\bar{t}}|)}{N(|y_{t}|>|y_{\bar{t}}|) + N(|y_{t}|<|y_{\bar{t}}|)},
\end{equation*}

where $N$ is the number of events and $y_{t}$ ($y_{\bar{t}}$) the rapidity of the top quark (top antiquark). The production of \ttbar events is predicted to be symmetric under the exchange of top quark and
antiquark, i.e.\ $\Ac= 0$, at leading-order (LO) accuracy in perturbative QCD. However, at next-to-leading order (NLO), quark--antiquark-initiated \ttbar production is asymmetric in the top-quark rapidity distribution, owing to interference between processes with initial- and final-state gluon emission and between the Born diagram and the box diagram at $O(\alphas^4)$. The total asymmetry from the sum of all effects is expected to be positive~\cite{PhysRevD.86.034026}.

Previous measurements of the asymmetry in \ttbar production by the ATLAS and CMS Collaborations at centre-of-mass energies ($\sqrt{s}$) of 7, 8 and 13~\TeV~\cite{TOPQ-2012-17,TOPQ-2014-16,CMS-TOP-11-030,CMS-TOP-12-033,TOPQ-2013-09,TOPQ-2015-08,CMS-TOP-12-010,CMS-TOP-15-009,TOPQ-2015-10,CMS-TOP-13-013,CMS-TOP-17-014} agree with the SM expectation. The most recent measurement of the inclusive and differential charge asymmetry at 13~\TeV\ by the ATLAS Collaboration~\cite{arXiv:2208.12095} reported evidence for a non-zero asymmetry in \ttbar production (measured as $\Ac=0.0068\pm0.0015$). While \Ac at the LHC corresponds to a central--forward asymmetry, the \ttbar production asymmetry manifests itself as a forward--backward asymmetry at the Tevatron. Early measurements of this asymmetry showed deviations from NLO QCD predictions, particularly at large values of the \ttbar invariant mass~\cite{Aaltonen:2011kc}. However, more recent results by the CDF and D0 Collaborations~\cite{Abazov:2015fna,Aaltonen:2012it,Abazov:2014cca} are compatible with the improved SM predictions including NLO electroweak (EW) and higher-order QCD corrections~\cite{PhysRevLett.115.052001,PhysRevD.91.071502}.

\begin{figure}[ht]
\centering
\includegraphics[width=0.25\textwidth]{fig_01a.pdf}
\quad \quad
\includegraphics[width=0.30\textwidth]{fig_01b.pdf}
\caption{Example Feynman diagrams of \tty production contributing to the charge asymmetry.}
\label{fig:intro_diag_tty}
\end{figure}

The \ttbar charge asymmetry is diluted at the LHC owing to the large fraction of gluon--gluon-initiated \ttbar events, which are symmetric under the exchange of the top quark and antiquark. However, it is enhanced in other topologies where the fraction of quark--antiquark-initiated production is larger, such as in associated production of \ttbar with a photon (\tty)~\cite{Aguilar-Saavedra:2014vta, Aguilar-Saavedra:2014kpa}. Interference effects among QCD diagrams at NLO, similar to those in \ttbar events, are predicted in \tty production. However, the dominant contribution to the asymmetry in \tty arises from interference between QED initial-state radiation, Figure~\ref{fig:intro_diag_tty}(left), and final-state radiation, Figure~\ref{fig:intro_diag_tty}(right), which yields a larger asymmetry of negative sign. In addition, other QCD--EW higher-order contributions can have a sizeable effect on the observed asymmetry~\cite{Pagani:2021}. The overall asymmetry in \tty at $\sqrt{s}= 13$~\TeV\ is expected to have a negative value, of 1\%--2\% depending on the phase space, according to SM predictions~\cite{Schulze:2019ttyAc, Pagani:2021} and can be modified by \enquote{beyond-the-SM} contributions. Contributions from an $s$-channel colour octet or a $Z'$ boson would, for instance, result in a smaller \Ac absolute value~\cite{Aguilar-Saavedra:2014vta}.

These sources of asymmetry are only present in the \tty events where the photon is radiated from an initial-state parton or one of the top quarks (hereafter referred to as \textit{\tty production}). The asymmetry is diluted by \tty events where the photon arises from any of the charged decay products of the \ttbar system (\textit{\tty decay} in the following). Therefore, only \tty production events are considered as signal in this analysis.

This \paperornote presents the first measurement of the charge asymmetry of the top-quark pairs in \tty production in \ttbar single-lepton final states, which have one high-\pt lepton and at least four jets, two of which are $b$-jets. It is performed using the full 139~\ifb\ data set recorded with the ATLAS detector between 2015 and 2018 at $\sqrt{s}=13~\TeV$, referred to as \RunTwo. In order to extract the asymmetry, the top quarks are reconstructed using a kinematic likelihood fit. The separation between signal and background processes is enhanced using a neural network (NN) approach. The output distribution of the NN is used to define two regions, one enriched in background events and one in signal events. The \Ac value is determined by means of a maximum-likelihood fit to the distribution of the difference of absolute rapidities of the top quark and antiquark. This is also referred to as \enquote{maximum-likelihood unfolding}.

The \paperornote is organised as follows. The ATLAS detector is briefly introduced in Section~\ref{sec:detector}. The simulation of signal and background processes is summarised in Section~\ref{sec:modelling}. The event reconstruction,  selection, and estimation of the background processes are presented in Sections~\ref{sec:selection} and~\ref{sec:background}. The systematic uncertainties are described in Section~\ref{sec:systematics}. The analysis strategy is discussed in Section~\ref{sec:analysis}, followed by the result in Section~\ref{sec:result}. Finally, a summary is given in Section~\ref{sec:conclusion}.

\section{ATLAS detector}
\label{sec:detector}
\newcommand{\AtlasCoordFootnote}{
ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre
of the detector and the $z$-axis along the beam pipe. The $x$-axis points from the IP to the centre of the
LHC ring, and the $y$-axis points upwards. Cylindrical coordinates $(r,\phi)$ are used in the transverse plane,
$\phi$ being the azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of the polar angle
$\theta$ as $\eta = -\ln \tan(\theta/2)$, and the rapidity as $y = (1/2)[(E+p_z)/(E-p_z)]$.
Angular distance is measured in units of $\Delta R \equiv \sqrt{(\Delta\eta)^{2} + (\Delta\phi)^{2}}$.}

The ATLAS~\cite{PERF-2007-01,ATLAS-TDR-19,Abbott:2018ikt} detector is a multipurpose detector with a
forward--backward symmetric cylindrical geometry with respect to the LHC beam axis.\footnote{\AtlasCoordFootnote}
The innermost layers consist of tracking detectors in the pseudorapidity range
$|\eta|<2.5$. This inner detector (ID) is surrounded by a thin superconducting solenoid
that provides a \SI{2}{\tesla} axial magnetic field. It is enclosed by the electromagnetic and
hadronic calorimeters, which cover $|\eta|<4.9$. The outermost layers of ATLAS consist of an external
muon spectrometer within  $|\eta|<2.7$, incorporating three large toroidal magnetic assemblies with eight coils each.
The field integral of the toroids ranges between 2.0 and 6.0$\,$Tm for most of the acceptance.
The muon spectrometer includes precision tracking chambers and fast detectors for triggering.
A two-level trigger system~\cite{TRIG-2016-01} reduces the recorded event rate to an average of \SI{1}{\kilo\hertz}. An extensive software suite~\cite{ATL-SOFT-PUB-2021-001} is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

\section{Simulation of signal and background processes}
\label{sec:modelling}
Monte Carlo (MC) event generators were used to estimate the contributions from the expected signal and background processes. The response of the ATLAS detector was simulated~\cite{SOFT-2010-01} with \GEANT~\cite{Agostinelli:2002hh}.
A fast simulation (\textsc{AtlFast-II}), which relies on a parameterisation of the calorimeter response~\cite{ATL-PHYS-PUB-2010-013}, was used in samples employed to estimate uncertainties related to the \ttbar and \tWy modelling. The additional \pp~collisions in the same or neighbouring bunch crossings, referred to
as pile-up, were generated with \PYTHIA[8.186]~\cite{Sjostrand:2007gs} using a set of tuned parameters called the A3 tune~\cite{ATL-PHYS-PUB-2016-017} and the \NNPDF[2.3lo] parton distribution function (PDF) set~\cite{Pumplin:2002vw}.

The signal \tty production events were simulated with \MGNLO[2.7.3]~\cite{Alwall:2014hca} as a $2\to 3$ process at NLO accuracy in QCD. The interference effects between initial-state and final-state photon radiation were considered in this process. The final-state top quarks in the \tty production sample are on-shell and were decayed at LO using \MADSPIN~\cite{Frixione:2007zp,Artoisenet:2012st} to preserve spin correlations.

The background \tty decay events, where the photon arises from any of the decay products of the top quarks or one of the on-shell top quarks, were simulated with the same version of \MGNLO but at LO precision as a $2 \to 2$ process followed by the decay of the top quarks, also simulated at LO precision. Both samples were generated using the \NNPDF[3.0nlo]~\cite{Ball:2014uwa} PDF set and interfaced to \PYTHIA[8.240]~\cite{Sjostrand:2014zea}, which used the A14 tune~\cite{ATL-PHYS-PUB-2014-021} and the \NNPDF[2.3lo] PDF set to model the parton shower, hadronisation, fragmentation and underlying event. The renormalisation and factorisation scales were set to $0.5\times \sum_i \sqrt{m^2_i+p^2_{\text{T},i}}$, where $m_i$ and $p_{\text{T},i}$ are the masses and transverse momenta of the particles generated from the matrix element (ME) calculation. Photons are required to have $\pT > 15~\GeV$ and to be isolated according to a smooth-cone hadronic isolation criterion with $\delta_0=0.1$, $\epsilon_{\gamma}=0.1$ and $n=2$, defined in Ref.~\cite{Frixione:1998jh}, which avoids infrared divergences. The top-quark mass in the \tty sample and all other samples involving top quarks was set to \SI{172.5}{\gev} and the decays of bottom and charm hadrons were simulated using the \EVTGEN[1.6.0] program~\cite{Lange:2001uf}.

The \tty production sample is normalised to the NLO cross section given by the MC simulation, while the normalisation of the \tty decay sample is corrected by a NLO/LO inclusive $K$-factor of 1.5. This $K$-factor was derived by comparing the normalisation of the sum of the NLO \tty production sample and the LO \tty decay sample with the normalisation of a LO inclusive $2 \to 7$ \tty sample corrected with the $K$-factor obtained in Ref.~\cite{TOPQ-2017-14} using the calculation described in Ref.~\cite{Melnikov:2011ta}.

The \ttbar events were simulated at NLO accuracy in QCD using \POWHEGBOX[v2]~\cite{Nason:2004rx,Frixione:2007vw,Alioli:2010xd} and the \NNPDF[3.0nlo] PDF set. The parton shower was  generated with \PYTHIA[8.230] using the A14 tune~\cite{ATL-PHYS-PUB-2012-003}. The \ttbar events are normalised to a cross-section value calculated with the \TOPpp[2.0] program at next-to-next-to-leading order (NNLO) in perturbative QCD, including soft-gluon resummation to next-to-next-to-leading-logarithm order (see Ref.~\cite{Czakon:2011xx:mod} and references therein).

The \tWy events were generated at LO accuracy with the \MGNLO[2.7.3] generator in the five-flavour scheme. To simulate this process, two complementary samples were generated; one as a $2\rightarrow3$ process assuming a stable top quark and the other as a $2 \to 2$ process, where the photon is radiated from any other charged final-state particle. To avoid infrared divergences, the photon was required to have $\pT > 15~\GeV$ and $|\eta| < 5.0$ and to be separated by $\Delta R>0.2$ from any parton. Both samples make use  of the \NNPDF[2.3lo] PDF set and were interfaced to \PYTHIA[8.212] for parton showering using the A14 tune.

Single-top-quark $s$- and $t$-channel production and inclusive $tW$ production were simulated at ME level at NLO in QCD with \POWHEGBOX[v2] and the \NNPDF[2.3lo] PDF set. The event generator was interfaced to \PYTHIA[8.230], which used the A14 tune. The events are normalised to the NNLO cross section~\cite{Kidonakis:2010tc,Kidonakis:2010ux,Kidonakis:2011wy}.

Events with $W\gamma$~and $Z\gamma$~final states (with additional jets) were simulated with \SHERPA[2.2.8]~\cite{Bothmann:2019yzt} at NLO in QCD using the \NNPDF[3.0nnlo] PDF set. The samples are normalised to the cross sections given by the MC simulation. The \SHERPA generator performs all steps of the event generation, from the hard process to the observable particles.
Events with inclusive $W$- and $Z$-boson production in association with additional jets were simulated with \SHERPA[2.2.1]~\cite{Bothmann:2019yzt,Hoeche:2009rj} at NLO in QCD.
The \NNPDF[3.0nlo] PDF set was used together with a dedicated tune provided by the \SHERPA authors.
The samples are normalised to the NNLO cross section in QCD~\cite{STDM-2015-03}.

Diboson processes, $\mathit{WW}$, $\mathit{WZ}$ and $\mathit{ZZ}$, were generated with \SHERPA[2.2.2] (leptonic decays) and~2.2.1 (non-leptonic final states) at LO in QCD. The \NNPDF[3.0nnlo] PDF set was used with a dedicated tune provided by the \SHERPA authors.
The samples are normalised to NLO cross sections in QCD~\cite{Campbell:1999ah}.

Events with a \ttbar pair and an associated $W$~or $Z$~boson ($t\bar{t}V$) were simulated at NLO at the ME level with \MGNLO using the \NNPDF[3.0nlo] PDF set. The ME generator was interfaced to \PYTHIA[8.210], for which the A14 tune was used in conjunction with the \NNPDF[2.3lo] PDF set. The samples are normalised to NLO in QCD and electroweak theory~\cite{deFlorian:2016spz}.

A procedure was applied to remove the overlap between the samples in which events were generated at ME level without explicitly including a photon in the final state and the dedicated samples where photons were included in the ME-level event-generation step (\tty and \Wty final states, as well as $W\gamma$ and $Z\gamma$ final states with additional jets). Events in the inclusive samples are discarded if they contain a parton-level photon that fulfils $\pT(\gamma)>\SI{15}{\gev}$ and $\Delta R(\gamma, \ell) > 0.2$, where $\pT(\gamma)$ is the transverse momentum of the photon and $\Delta R(\gamma, \ell)$ is the angular distance between the photon and any charged lepton.

Corrections to the pile-up profile, the trigger, reconstruction and selection efficiencies, and the energy scales and resolutions, are applied to the MC simulation samples to improve the description of the data.

\section{Event reconstruction and selection}
\label{sec:selection}

The analysis uses a data set that passes stringent quality requirements and corresponds to an integrated luminosity of 139~\ifb{} collected with the ATLAS detector during \RunTwo of the LHC. Events are required to have at least one primary vertex reconstructed from at least two associated tracks, and only events where at least one single-electron~\cite{TRIG-2018-05} or single-muon~\cite{TRIG-2018-01} trigger was fired are selected.

Muons are reconstructed by combining a track in the muon spectrometer with a track in the ID system.
The reconstruction, identification and calibration methods are described in Ref.~\cite{MUON-2018-03}. The muon track is also required to originate from the primary collision vertex. Muons are required to be isolated according to measurements of nearby track \pT and calorimeter energy.
Only muons with calibrated $\pT > \SI{25}{\GeV}$ and $|\eta| < 2.5$ and passing `medium' quality requirements are considered.

Electrons are reconstructed from energy deposits in the electromagnetic calorimeter (ECAL) associated with reconstructed tracks in the ID system. The origin of the electron track also has to be compatible with the primary vertex. Electrons are identified with a combined likelihood technique~\cite{EGAM-2018-01} using a `tight' working point, and are required to be isolated  according to measurements of nearby calorimeter energy and track \pT. Electrons are calibrated with the method described in Ref.~\cite{EGAM-2018-01} and are selected if they fulfil $\pt>\SI{25}{\GeV}$ and $|\eta_{\mathrm{clus}}| < 2.47$, excluding the ECAL barrel/endcap transition region $1.37 < |\eta_{\mathrm{clus}}| < 1.52$, where $\eta_{\mathrm{clus}}$ refers to the pseudorapidity of the calorimeter energy cluster associated with the electron.

Photons are reconstructed from energy deposits in the central region of the ECAL~\cite{EGAM-2018-01}. Photons are required to fulfil tight identification and isolation requirements. The latter is defined as $\left. \ET^{\mathrm{iso}} \right|_{\Delta R < 0.4} < 0.022 \cdot \ET(\gamma) + \SI{2.45}{\GeV}$ in conjunction with $\left. \pT^{\mathrm{iso}} \right|_{\Delta R < 0.2} < 0.05 \cdot \ET(\gamma)$, where $\ET^{\mathrm{iso}}$ refers to the calorimeter isolation within a cone of size $\Delta R = 0.4$ around the direction of the photon candidate and $\pT^{\mathrm{iso}}$ is the track isolation within a $\Delta R = 0.2$ cone~\cite{PERF-2017-02}.
Photons are required to have transverse energy $\ET > \SI{20}{\GeV}$ and $|\eta_{\mathrm{clus}}| < 2.37$, excluding the calorimeter transition region. They are separated into two categories, one where the cluster is not matched to any reconstructed track in the ID system (unconverted photons) and the other where the cluster is matched to one or two reconstructed tracks that are consistent with originating from a photon conversion and, in addition, a conversion vertex can be found (converted photons).

Jets are reconstructed using the anti-$k_{t}$ algorithm~\cite{Cacciari:2008gp} in the \Fastjet implementation~\cite{Fastjet} with a distance parameter $R=0.4$. Their reconstruction is performed on particle-flow objects~\cite{PERF-2015-09}. The jet energy scale and jet energy resolution are calibrated using an energy- and $\eta$-dependent calibration scheme based on simulation with in situ corrections obtained from data~\cite{PERF-2016-04}. Only jets with $\pT > \SI{25}{\GeV}$ and $|\eta| < 2.5$ are considered in the analysis.
Jets with a large contribution from pile-up vertices are identified with the \emph{jet vertex tagger} (JVT)~\cite{ATLAS-CONF-2014-018} and rejected.

Jets arising from \bquark hadronisation, referred to as \bjets, are identified using the DL1r \btag algorithm~\cite{ATLAS:2022qxm}, which is based on an artificial deep neural network combining information from other algorithms using track impact parameters and secondary vertices, and a multi-vertex reconstruction algorithm. The flavour-tagging efficiency for \bjets, as well as for $c$-jets and light-flavour jets, is calibrated as described in Ref.~\cite{FTAG-2018-01}. The working point used to select the \bjets corresponds to a selection efficiency of 77\% in simulated \ttbar events.

The magnitude of the reconstructed missing transverse momentum (\met)~\cite{ATLAS-CONF-2018-023,PERF-2016-07} is calculated from the negative vector sum of the \pT of all calibrated physics objects and the remaining unclustered energy, also called the \emph{soft term}. This term is estimated from low-\pt tracks associated with the primary vertex but not with any reconstructed object.

An overlap removal procedure is implemented to avoid the reconstruction of the same energy clusters or tracks as different objects. Electron candidates that share their track with a muon candidate are removed and jets within a $\Delta R = 0.2$ cone around any of the remaining electrons are excluded. If the distance between an electron and any remaining jet is $\Delta R < 0.4$, the electron is subsequently removed. In the next step, muon candidates within $\Delta R = 0.4$ of a jet are removed if the jet has more than two associated tracks, otherwise the jet is discarded. In the final step, photons within a $\Delta R = 0.4$ cone around any remaining electron or muon are excluded and then jets within a $\Delta R = 0.4$ cone around any remaining photon are removed.

Events are selected if they have exactly one electron or one muon that is matched to the corresponding trigger-level object. The \pt thresholds for the leptons are $25~\GeV$ in 2015 data, $27~\GeV$ in 2016 data, and $28~\GeV$ in 2017 and 2018 data, which are at least 1~\GeV\ above the \pt thresholds of the single-lepton triggers. This is done to avoid differences due to the calibration of the objects used in the trigger logic, and objects used in the physics analysis.
Only events containing exactly one reconstructed photon fulfilling the condition $\Delta R(\ell,\gamma)>0.4$ are considered in the measurement. Additionally, events where the invariant mass of the electron--photon system is within $5~\GeV$ of the $Z$-boson mass are rejected. The event is also required to have at least four jets and at least one of them must be \btagged.

The kinematic properties, in particular the rapidities of the top quark and antiquark, are determined by means of a constrained kinematic fitting algorithm, KLFitter~\cite{KLFitter:2014}, based on a maximum-likelihood approach applied to the four-momenta of the selected lepton and up to five leading jets (\pT-ordered) and \met, representing the transverse momentum of the neutrino. The likelihood is constructed as the product of transfer functions that relate the energies of the reconstructed objects and parton-level objects and Breit--Wigner distributions that match the selected objects to the $W$ bosons and the top quarks. The fit is constrained to reconstruct two $W$ bosons, each with a mass of 80.4~\GeV. In addition, the reconstructed masses of the top quark and antiquark are constrained to 172.5~\GeV. The combination of jets that gives the highest likelihood is selected and the jets are used to reconstruct the hadronically and leptonically decaying top quark or antiquark. In the case of leptonic top-quark decay, the invariant mass of the lepton--neutrino system is constrained to be the $W$-boson mass and a quadratic equation for the neutrino's longitudinal momentum is obtained. For real solutions of the equation, the solution that results in the top-quark mass closer to 172.5~\GeV\ is chosen, while in the case of complex solutions, the real part of the solution is considered. The fraction of top quarks that are reconstructed within $\Delta R=1.0$ of the parton-level top quarks is about 64\% for hadronically decaying top quarks and 72\% for leptonically decaying top quarks.

\section{Background estimation}
\label{sec:background}
Each background process is assigned to one of several categories depending on the origin of the reconstructed photon, whether it is a prompt photon or whether another object mimics a photon signature. The estimation of background events from different sources closely follows the methods employed in Ref.~\cite{TOPQ-2017-14}.

After the event selection, the largest background contribution is that of \tty decay events (about 30\% of the total number of events). The \emph{\catprompty background} category, which contains any other type of background process with a prompt photon, constitutes about 15\% of the selected events. Both contributions are estimated using MC simulation. Agreement between data and simulation of the \catprompty background was validated in dedicated regions. A validation region enriched in the $Z\gamma$ process is defined by selecting events with exactly one photon, two same-flavour opposite-sign leptons, fewer than four jets and no $b$-tagged jets. Events in the $W\gamma$ validation region fulfil the same lepton and photon requirements as in the signal region. The simulation of $W\gamma$ is known to underestimate the number of events with heavy-flavour jets in data. Thus, to define a region orthogonal to the signal region while selecting events with heavy-flavour content, the events in the validation region are additionally required to have fewer than four jets and at least one must pass the $b$-tagging working point with a selection efficiency of 85\% but fail the one with 70\% efficiency in simulated \ttbar events. The expected fractions of $Z\gamma$ and $W\gamma$ events in the corresponding control regions are about 95\% and 50\%, respectively.

Another significant contribution arises from processes with an electron mimicking a photon signature in the detector, referred to as \emph{e-fake}. This background contribution amounts to 16\% of the total number of selected events and it is estimated from data by applying a tag-and-probe method to $Z\rightarrow e^+e^-$ events. Two control regions are defined in order to determine the \efake rate in data and simulation: one region contains events with an electron--positron pair and the other, enriched in \efake events, contains events with an electron and a photon satisfying the object selection criteria described in Section~\ref{sec:selection}. Additionally, the invariant mass of the pair of objects is required to be in the range [40, 140]~\GeV, and their angular separation in $\phi$ must exceed 2.62~rad to suppress the contributions from events where the photon is radiated from the electron. Background contributions not originating from $Z$-boson events are subtracted in data with a fit of the invariant mass
distribution. The $Z \rightarrow e^+e^-\gamma$ contribution with prompt photons where one of the electrons is not reconstructed or identified is subtracted using simulation. The \efake rate is obtained as the ratio of the event yield after background subtraction in the e-fake-enriched control region to the yield in the electron--positron control region. The calculated ratio is binned in the \pt (three bins) and $|\eta|$ (four bins) of the photon in the $e\gamma$ events and either the electron or the positron in the $e^+e^-$ events, selected randomly to avoid biasing the selection, and separately for converted and unconverted photons. Scale factors are calculated to correct the e-fake background estimate in the signal region, based on a comparison between the e-fake photon rates obtained using either data or simulation. The systematic uncertainties in the scale factors account for possible mismodelling of the signal and background processes in the fit. The values of the scale factors vary from 0.8 to 1.4 with uncertainties between 5\% and 20\%.

The background contribution from events where the photon signature arises from hadronic energy depositions in the ECAL or from hadron decays such as $\pi^0 \to \gamma\gamma$, generically referred to as \emph{\cathfake}, constitutes about 7\% of the events. The \cathfake background is estimated from data by using the so-called ABCD method. Three orthogonal regions enriched with \cathfake photon events are defined by inverting the photon isolation selection and the requirements on four variables related to the shower shape in the first layer of the ECAL, which are part of the photon tight identification criteria. They are chosen because of their small correlation with the photon isolation and their power to discriminate between prompt and \cathfake photons. Events are selected for regions A and B if their photon fails at least two out of four identification requirements, while satisfying all other identification criteria, and pass or fail the isolation requirements, respectively. Region C contains events where the photons fail the isolation requirements but satisfy the tight identification criteria. Additionally, the sum of the \pt of all tracks within $\Delta R=0.2$ of the photon is required to be larger than 3~\GeV\ to further suppress the prompt-photon contribution in regions B and C. The \cathfake background contribution in the signal region is measured as the product of the numbers of events in regions A and C divided by the number of events in region B. The estimate is corrected for the correlation between the criteria, which is obtained using MC simulations. The scale factors are obtained separately for converted and unconverted photons and as a function of the photon \pt (two bins) and $|\eta|$ (four bins). The considered sources of systematic uncertainty in the \cathfake background contribution include the modelling of the \ttbar process, which contributes about 90\% of the \cathfake events, the shower shapes, and the normalisation uncertainties of the background processes. The scale factors range from 0.6 to 1.5, with uncertainties of 30\%--60\%.

The contribution from events with a non-prompt or misidentified lepton, referred to as \emph{lepton fake}, is obtained using the data-driven approach referred to as the matrix method~\cite{ATLAS:2022swp}. Events are separated into two categories that are based on tighter or looser lepton identification and isolation requirements, and thus enriched in events with real leptons or non-prompt/fake leptons, respectively. The contribution in the signal region is estimated from the data events passing loose lepton selection requirements, corrected by a weight that depends on the real- and fake-lepton efficiencies obtained from the two event categories described above. The efficiencies are parameterised as a function of the lepton kinematic properties. The uncertainties are estimated by using different parameterisations and tighter control regions, and by including the normalisation uncertainty of the prompt-lepton background. This background contribution amounts to around 1\% of all selected events, dominated by events with a misidentified electron.

\section{Systematic uncertainties}
\label{sec:systematics}
The precision of the obtained asymmetry \Ac is affected by several sources of systematic uncertainty, arising from detector effects or theoretical assumptions, as well as uncertainties due to the limited number of events in the MC simulations. The different sources of systematic uncertainty, discussed in the following, affect the event yields, the distribution shape of the observable of interest, or both. The sources of systematic uncertainty affecting the shape of the distribution typically have a larger impact on the precision of the result because global normalisation uncertainties cancel out in the ratio of event yields that defines \Ac.

The experimental systematic effects include uncertainties in the integrated luminosity and the simulation of pile-up events, as well as effects related to the reconstruction and identification of the physics objects in the analysis.

The uncertainty in the total integrated luminosity is estimated to be 1.7\%~\cite{ATLAS-CONF-2019-021}, using the LUCID-2 detector \cite{LUCID2} for the primary luminosity measurements.  The uncertainty associated with the modelling of pile-up is determined by varying the pile-up reweighting in the simulation within its uncertainties.

The photon and lepton identification and isolation efficiencies, momentum scale and resolution~\cite{PERF-2017-03,PERF-2015-10}, and lepton trigger efficiencies in simulation are corrected using scale factors to better describe the corresponding values in data. These corrections, which typically depend on \pT and $\eta$, are varied within their uncertainties to estimate the corresponding systematic uncertainty.

The jet energy scale (JES) uncertainty is derived from a combination of simulations, test-beam data and in situ measurements~\cite{PERF-2016-04}. Contributions from jet-flavour composition, $\eta$-intercalibration, punch-through, single-particle response, calorimeter response to different jet flavours, and pile-up are also taken into account, yielding a total of 30 uncorrelated JES uncertainty subcomponents, of which 29 are non-zero in a given event depending on the type of simulation used. The jet energy resolution in simulation is smeared by its corresponding uncertainty~\cite{ATL-PHYS-PUB-2015-015} split into eight uncorrelated sources. The uncertainty associated with the JVT discriminant for pile-up jet rejection is obtained by varying the efficiency correction factors.

The uncertainties in the $b$-jet tagging calibration are determined separately for $b$-jets, $c$-jets and light-flavour jets~\cite{PERF-2016-05,FTAG-2020-08,ATLAS-CONF-2018-006}. For each jet category, the uncertainties are decomposed into several uncorrelated components.

The uncertainty in \MET arises from the propagation of the energy scales and resolutions of photons, leptons and jets, and the modelling of its soft term~\cite{PERF-2016-07}.

The signal and background modelling uncertainties include those owing to the choice of QCD scales, parton shower, amount of QCD initial-state radiation (ISR), and PDF set.

The effect of the QCD scale uncertainty for each of the \tty, \Wty and \ttbar processes is evaluated independently by separately halving and doubling the renormalisation and factorisation scales relative to the default scale choice. The uncertainty from the parton shower and hadronisation for those processes is estimated by comparing the nominal simulated samples interfaced with \PYTHIA[8] with alternative samples interfaced to \HERWIG[7]~\cite{Bahr:2008pv,Bellm:2015jjp}. Uncertainties due to the value of \alphas used in the ISR parton shower modelling are estimated by comparing the nominal \tty, \Wty and \ttbar simulations with alternative samples simulated with higher or lower radiation parameter settings in the A14 tune, controlled by the \textit{var3c} parameter implemented in \PYTHIA[8].
An additional ISR uncertainty is obtained for the \ttbar process by comparing the nominal sample with an additional one with the \hdamp parameter, which controls the \pt of the first additional emission, varied by a factor of two~\cite{ATL-PHYS-PUB-2018-009}. The uncertainty in the PDFs for the signal and background \tty processes is estimated using the 30 PDF variations of the PDF4LHC15 prescription~\cite{Butterworth:2015oua}. The PDF variations are propagated by using alternative generator weights and each of them is considered as a separate nuisance parameter in the fit.

For the e-fake, h-fake, and lepton-fake background contributions, the total uncertainties associated with the corrections obtained using data are propagated to the final result. A normalisation uncertainty of 20\% is assigned to the \tty decay process, based on the estimated uncertainty in the NLO $K$-factor~\cite{TOPQ-2017-14}, and a 50\% normalisation uncertainty is assigned to $W\gamma$, based on the differences between data and simulation observed in the dedicated control region, and to the minor background processes contributing to the prompt-photon category, i.e.\ single top quark, \ttV, diboson, and $Z\gamma$.

\section{Signal discrimination}
\label{sec:analysis}
A multivariate analysis using a neural network is performed to further separate the \tty signal from the background processes. The NN is fully connected and consists of three hidden layers. The first two layers consist of 96 nodes and are followed by a batch normalisation layer. The third layer has 16 nodes. The hidden layers use a parametric ReLU activation function, while the output node uses a sigmoid activation function. The training is performed with Keras~\cite{chollet2015keras} with the TensorFlow~\cite{tensorflow2015-whitepaper} backend with binary cross-entropy as a loss function. The overall set of events is split into a training set (68\%), a set (17\%) for validation of the NN configuration, used in a 5-fold cross-validation, and a testing set (15\%). The event weights are applied to the events in the training, testing and validation samples. The NN uses a total of 21 variables related to the kinematic properties of individual objects, such as the photon \pt and $\eta$, event shape variables (e.g.\ \met and the scalar sum of the \pt of the jets in the event), the number of $b$-tagged jets, the pseudo-continuous binned $b$-tagging discriminant~\cite{FTAG-2018-01}, the photon conversion type, and invariant masses and angular separations of different objects in the event (e.g.\ the invariant mass and $\Delta R$ of the lepton or the photon and the closest $b$-tagged jet). Example variables from among those with the most discriminating power are shown in Figure~\ref{fig:prefit_variables_1}. The MC simulation describes the data within the uncertainties. The largest contributions to the MC uncertainty band are the normalisation uncertainties associated with the backgrounds with prompt photons.

\begin{figure}[ht]
\centering
\includegraphics[width=0.3\textwidth]{fig_02a.pdf}
\includegraphics[width=0.3\textwidth]{fig_02b.pdf}
\includegraphics[width=0.3\textwidth]{fig_02c.pdf}
\caption{Distributions of photon \pt (left), angular separation of the lepton and the closest $b$-tagged jet (middle), and transverse mass of the leptonically decaying $W$ boson (right) before the fit. The uncertainty band includes all experimental and modelling systematic uncertainties (cf.\ Section~\ref{sec:systematics}) added in quadrature. Overflow events are included in the last bin of each distribution. The lower part of the plot shows the ratio of the data to the prediction.}
\label{fig:prefit_variables_1}
\end{figure}

The resulting NN discriminant output, $O_{\mathrm{NN}}$, shown in Figure~\ref{fig:prefit_variables_NN}, is used to divide the events into a background-enriched region and a signal-enriched region, defined by $O_{\mathrm{NN}} < 0.6$ and $O_{\mathrm{NN}} \geq 0.6$, respectively. This threshold was optimised to give the smallest expected uncertainty in \Ac. The observed and expected signal and background event yields in the two regions are summarised in Table~\ref{tab:prefityields}. The slight underestimate of the data by the SM prediction, observed in Figure~\ref{fig:prefit_variables_1}, is reflected in Figure~\ref{fig:prefit_variables_NN} at large values of $O_{\mathrm{NN}}$ because it is expected to partially come from the normalisation of the \tty production simulation, which is a free parameter in the profile likelihood unfolding described in the following.

\begin{figure}[ht]
\centering
\includegraphics[width=0.45\textwidth]{fig_03.pdf}
\caption{Distribution of the NN output discriminant before the fit. The uncertainty band includes all experimental and modelling systematic uncertainties (cf.\ Section~\ref{sec:systematics}) added in quadrature. The lower part of the plot shows the ratio of the data to the prediction.}
\label{fig:prefit_variables_NN}
\end{figure}

\begin{table}
\centering
\caption{Event yields before the profile likelihood unfolding after the full selection in the two regions defined by the NN discriminant value. The quoted uncertainties correspond to all statistical and systematic uncertainties (cf.\ Section~\ref{sec:systematics}) added in quadrature.}
\label{tab:prefityields}
\sisetup{table-format=4.0(3)}
\footnotesize
\begin{tabular}{l|r|r}
\hline
& {$O_{\mathrm{NN}} < 0.6$} & {$O_{\mathrm{NN}} \geq 0.6$}\\
\hline
\tty prod (signal)   	& 6660 	$\pm$  \phantom{0}350 	& 6910 	$\pm$  \phantom{0}340 \\
\tty decay   		& 14\,100 	$\pm$ 3100 			& 1900 	$\pm$  \phantom{0}560 \\
h-fake $\gamma$   	& 3400 	$\pm$ 1400 			& 790 	$\pm$  \phantom{0}360 \\
e-fake $\gamma$   	& 6420 	$\pm$  \phantom{0}860 	& 1480 	$\pm$  \phantom{0}260 \\
Prompt $\gamma$  & 6400 	$\pm$ 2000 			& 1300 	$\pm$  \phantom{0}400 \\
Lepton fake   		& 410 	$\pm$  \phantom{0}110 	& 57 		$\pm$  \phantom{00}35 \\
\hline
Total  & \numRF{37405.9}{3} $\pm$ \numRF{4543.05}{2} & \numRF{12438.5}{3} $\pm$ \numRF{1055.71}{2} \\

\hline
Data   & 38\,527\phantom{~$\pm$~0000} & 13\,763\phantom{~$\pm$~0000} \\

\hline
\end{tabular}
\end{table}

\section{Results}
\label{sec:result}
The value of \Ac is extracted from the $|y_{t}| - |y_{\bar{t}}|$ distribution in a fiducial region defined at particle level. The top quark and antiquark are defined at parton level in the MC simulation after final-state radiation but before decay. The fiducial region at particle level is defined by applying selection requirements similar to those at reconstruction level to the stable particles after the event generation and before the detector simulation. The fiducial phase space is defined by requiring exactly one photon, exactly one electron or muon, and at least four jets, of which at least one must be a $b$-jet, defined as follows. Photons are required to not originate from a hadron decay, to have $\et > \SI{20}{\GeV}$ and $|\eta| < \SI{2.37}{}$, and to be isolated such that the sum of transverse momenta of all charged particles surrounding the photon within $\Delta R \leq 0.2$ must be less than 5\% of its own \pt. Muons and electrons must have $\pt>25$~\GeV\ and $|\eta| < 2.5$, and must not originate from hadron decays. The momenta of nearby photons, within a $\Delta R = 0.1$ cone, are added to the lepton before applying the selection. Jets are clustered with the anti-$k_{t}$ algorithm with a radius parameter of $R=0.4$. All stable particles are considered in the clustering, except for the selected electrons, muons and photons, and the neutrinos originating from the top quarks. Jets are required to have $\pt > 25~\GeV$ and $|\eta| < 2.5$. A particle-level jet is identified as a $b$-jet if a hadron with $\pT > \SI{5}{\GeV}$ containing a $b$-quark is matched to the jet through a ghost-matching method~\cite{Cacciari:2008gn}. Jets within $\Delta R = 0.4$ of lepton or isolated photon candidates are removed.

The \Ac value is obtained by means of a simultaneous maximum-likelihood unfolding of the $|y_{t}| - |y_{\bar{t}}|$ distributions in the two regions defined by the NN output discriminant. The efficiency of selecting and reconstructing an event that is generated in the fiducial phase space is about 30\% in the two regions, while the fraction of events that fulfil the selection at reconstruction level but fail the particle-level requirements is about 20\%. The fraction of events that are reconstructed in the $|y_{t}| - |y_{\bar{t}}|$ bin where they were generated is approximately 75\%. The parameters of interest, which float freely in the fit, are the signal strength of the bin $|y_{t}| - |y_{\bar{t}}|>0$ and \Ac, which replaces the signal strength of the other bin, using the following expression: $\Ac = (\mu_{+}T_{+} - \mu_{-}T_{-})/(\mu_{+}T_{+} + \mu_{-}T_{-})$. The signal strength $\mu_{+}$ $(\mu_{-})$ is defined as the ratio of the measured cross section to the expected value given by the SM simulation in the bin $|y_{t}| - |y_{\bar{t}}|>0~(<0)$ at generator level. The variable $T_{+/-}$ represents the number of events at particle level in the corresponding bin. No regularisation is applied.  The systematic uncertainties are taken into account via nuisance parameters in the likelihood function. They are symmetrised, taking half of the difference between the upward and downward variations as the uncertainty. If both variations have the same sign, the average of the difference between each variation and the nominal value is chosen as the two-sided uncertainty. If only one variation is available (e.g.\ the uncertainty from the parton shower or the PDF uncertainties) the difference from the nominal value is taken as both the upward and downward variations for the corresponding source. In addition, a pruning procedure is implemented in order to remove the smallest systematic uncertainties.

In addition to the NLO/LO $K$-factor applied to correct the normalisation of the \tty decay process, the $|y_{t}| - |y_{\bar{t}}|$ template is reweighted to account for the asymmetry in \ttbar production. The weight is obtained at parton level by considering the central value of the prediction for the inclusive \ttbar asymmetry, $\Ac^{\ttbar} = 0.0064^{+0.0005}_{-0.0006}$, calculated at NNLO accuracy in QCD with EW corrections at NLO~\cite{PhysRevD.98.014003,arXiv:2208.12095}. The weight is then propagated to the distribution at reconstruction level. For consistency, the MC templates obtained with the NLO \ttbar samples are also reweighted to match that asymmetry value. To probe the robustness of the method, and in particular to verify that the unfolding procedure does not bias the result towards the asymmetry in the signal MC simulation, the measurement was repeated with pseudodata. Several pseudodata sets were obtained by reweighting the $|y_{t}| - |y_{\bar{t}}|$ distribution of the signal \tty production sample to correspond to different \Ac values and by adding the background contributions. The unfolding procedure was repeated using the nominal simulation. The resulting values of the asymmetry were in good agreement with the true \Ac asymmetry of each pseudodata set and do not indicate any bias.

\begin{figure}[ht]
\centering
\includegraphics[width=0.45\textwidth]{fig_04a.pdf}
\includegraphics[width=0.45\textwidth]{fig_04b.pdf}
\caption{The distributions of $|y_{t}| - |y_{\bar{t}}|$ after the fit in the two regions defined by the NN output. Underflow and overflow events are included in corresponding bins of the distributions. The uncertainty band represents the total post-fit uncertainties. Correlations among uncertainties are taken into account as determined in the fit. The lower part of the plot shows the ratio of the data to the prediction.}
\label{fig:dy_postfit}
\end{figure}

The post-fit $|y_{t}| - |y_{\bar{t}}|$ distributions in the two regions are shown in Figure~\ref{fig:dy_postfit}. The \Ac of the \tty production process is expected to have a negative sign, while the asymmetry of the background contributions with a top-quark pair, i.e.\ \tty decay and \ttbar, is expected to be positive as discussed in Section~\ref{sec:intro}. Good agreement is observed between the data and the prediction after the fit. As a result of the fit, a few nuisance parameters are slightly constrained: the uncertainties in the normalisation of the \tty decay and $W\gamma$ backgrounds are reduced by 30\% and 15\%, respectively. In all cases, the best-fit values of the nuisance parameters are well within one standard deviation of their initial values.

The asymmetry is found to be $\Ac=-0.003 \pm 0.029 = -0.003 \pm 0.024(\mathrm{stat}) \pm 0.017(\mathrm{syst})$, assuming the SM \ttbar charge asymmetry of $\Ac^{\ttbar}=0.0064$~\cite{PhysRevD.98.014003}. The systematic uncertainty is derived from its squared value, calculated as the difference of the squares of the total uncertainty and the statistical uncertainty, obtained from a fit without systematic uncertainties. The \Ac value is compatible with the value obtained from the \MGNLO MC simulation in the same phase space, $\Ac=-0.014\pm0.001(\mathrm{scale})$. The precision of the result is limited by the statistical uncertainty. The impact of the different sources of systematic uncertainty, grouped in categories, is summarised in Table~\ref{tab:syst_impact}. The most relevant sources of systematic uncertainty are the MC statistical uncertainty of the prompt-photon background and the experimental systematic sources related to jets and \MET. The dependence of the measured \tty \Ac on $\Ac^{\ttbar}$ is estimated by repeating the measurement for different values of the \ttbar asymmetry in the range between 0 and $2\times \Ac^{\ttbar}$. The dependence is found to be linear and can be parameterised as $\Ac = -0.57\times \Ac^{\ttbar}+0.0005$.

\begin{table}[htbp]
\centering
\caption{
Summary of the impact of the systematic uncertainties on \Ac grouped into different categories. The quoted uncertainties are obtained by repeating the fit with certain sets of nuisance parameters fixed to their post-fit values, and calculating the squared uncertainties as the difference of the squares of the full-fit and repeated-fit uncertainties. The category \textit{Other experimental} includes uncertainties associated with leptons, pile-up and luminosity.}
\label{tab:syst_impact}
\footnotesize
\begin{tabular}{l S}
\hline
\textbf{Total uncertainty} & 0.029 \\
\hline
\textbf{Statistical uncertainty} & 0.024 \\
\hline
\textbf{MC statistical uncertainties} & \\
\midrule
Background processes & 0.008 \\
\tty production & 0.004 \\
\midrule
\textbf{Modelling uncertainties} & \\
\midrule
\tty production modelling & 0.003 \\
Background modelling & 0.002 \\
Prompt background normalisation & 0.002 \\
\midrule
\textbf{Experimental uncertainties} & \\
\midrule
Jet & 0.009 \\
Fake-lepton background estimate & 0.005 \\
\MET & 0.005 \\
Fake-photon background estimates & 0.003 \\
Photon & 0.001 \\
$b$-tagging & 0.001 \\
Other experimental & 0.004 \\
\hline
\end{tabular}
\end{table}

\FloatBarrier

\section{Conclusion}
\label{sec:conclusion}

This \paperornote presents a measurement of the top-quark pair charge asymmetry in \tty events using \intlumitot of $pp$ collision data at a centre-of-mass energy of 13~\TeV\ collected by the ATLAS experiment at the LHC. The selected events have exactly one photon, one lepton, and at least four jets, of which at least one is $b$-tagged. The inclusive charge asymmetry yields $\Ac=-0.003 \pm 0.029 = -0.003 \pm 0.024(\mathrm{stat}) \pm 0.017(\mathrm{syst})$, which is compatible with the Standard Model prediction within the uncertainties. The precision is limited by the statistical uncertainty.

\printbibliography

\end{document}

redpajama_set_name: RedPajamaArXiv


text:
\section*{Supplemental Material on First-principles-based method for electron localization: Application to monolayer hexagonal boron nitride}
\section{Details of the TMDCA@DFT Self-Consistency}

To study the interplay of disorder and electron interactions, we utilize the Anderson-Hubbard Hamiltonian of the form
\begin{equation}
    H = H_0+\sum_{i\alpha\sigma} V_{i\sigma}^\alpha n_{i\sigma}^\alpha + U\sum_{i\alpha} n^\alpha_{i \uparrow} n^\alpha_{i \downarrow},
    \label{e.1}
\end{equation}
where $n_{i\sigma}^\alpha$ is the number operator, and $i$, $\alpha$, and $\sigma$ are site, orbital, and spin indices, respectively. The first term $H_0$ describes the single-particle Hamiltonian. This is obtained by downfolding the Kohn-Sham Bloch functions using a set of symmetry-adapted Wannier functions~\cite{PhysRevB.56.12847,*Mostofi2008685} that can accurately describe the states around the Fermi level. In our case, these functions are boron and nitrogen $s$, $p$, and $d$ orbitals of 18-band spin-restricted $H_0$ that accurately describes the DFT band structure in the energy interval of interest. The second term represents the disorder modeled by a binary site potentials $V_{i\sigma}^\alpha=V_i\in\{0,W\}$, where 0($W$) represent the absence (presence) of a vacancy, with the vacancy potential $W$ being a lot greater than the material bandwidth.  The vacancy is generated by a random disorder configurations using the probability mass functions $P(V_i=W)=\delta$ and $P(V_i=W)=0$ for B and N sites, respectively, where the averaged boron vacancy concentration $\delta$ satisfies the stoichiometry B$_{1-\delta}$N of the disordered material. The last term describes the Coulomb repulsion $U$ between two electrons occupying site $i$ incorporated self-consistently into our formalism
via interacting, non-local cluster self-energy ($\Sigma_c[{ \tilde{G}}](i,j\neq i)$), which is up to second order in the perturbation expansion of the interactions, $U^2$.

As explained in the main text, we aim to obtain the single-particle Green function and the associated density of states (DoS). Obtaining these ``exactly'' for Eq.~\ref{e.1} is, however, not feasible. We instead, obtain approximate Green functions based on the typical medium dynamical cluster approximation (TMDCA). The TMDCA is a mean-field approach with an intrinsic order parameter based on the typical DoS that ``properly'' characterizes disordered and/or interacting electron systems even in the proximity of electron localization. The TMDCA self-consistency is based on a formalism consistent with the generalized dynamical cluster theory approaches to correlated electron systems~\cite{PhysRevB.63.125102,*PhysRevB.61.12739,*Jarrell01}. Through a set of self-consistency equations, the TMDCA approach maps a disordered lattice onto a finite cluster embedded in a typical medium, as illustrated in Fig.1(A) of the main text and explained therein, and described in details below.

\begin{wrapfigure}{l}{0.2\textwidth}
\center{\includegraphics[trim = 0mm 5mm 5mm 15mm,width=.2\textwidth]{SOPT_SE_Diagrams.eps}}

\caption{The first and second-order Feynman diagrams of the interacting self-energy between sites $i$ and $j$.}
\label{Fig:sopt}
\end{wrapfigure}

The self-consistency procedure follows:  (i) After obtaining $H_0$ from a converged, density functional theory (DFT) calculations by downfolding the Kohn-Sham Bloch functions using a set of symmetry-adapted Wannier functions~\cite{PhysRevB.56.12847,*Mostofi2008685}, we make an initial guess of a hybridization function $\Gamma(\mathbf{K})$. The hybridization function measures the escape rate of electrons from the cluster to the effective medium. (ii) We Fourier transform the hybridization function to real space, $\Gamma_{n,m} = \sum_{\mathbf{K}} \Gamma(\mathbf{K})\exp[i \mathbf{K}\cdot(\mathbf{R}_n-\mathbf{R}_m)]$ and then for a given vacancy configuration, form the cluster Green function $G^c(V)=(E-\Gamma -H_0-V)^{-1}$. (iii) Using $G^c(V)$, we construct a Hartree-corrected cluster Green function
${\cal \tilde{G}}^{-1}_c(V,U) = G^c(V)^{-1} + \epsilon_d(U)$, where $\epsilon_d(U)=\mu -Un_i/2$ and $n_i$ $=$ $-1/\pi \int_{-\infty}^0  \Im {\cal \tilde{G}}(i,i,E) d E$ is the site occupancy at zero temperature. An additional internal self-consistency is used to converge $n_i$ and ${\cal \tilde{G}}_c$. This internal self-consistency ensures that $n_i$ is numerically the same as the one that could be obtained using the fully dressed Green function, but, at a fraction of the computation cost~\cite{PhysRevB.92.201114}. This is important as it is the most computational intensive part of the self-consistency. Also, it enables the systematic incorporation of crossing diagrams (for $N_c>1$) from both the disorder due to vacancy and electron interactions at equal footing~\cite{PhysRevB.92.201114}. (iv) For a given electron interaction strength and randomly generated disorder due to boron vacancy configuration $V$, the fully dressed cluster Green function $G^c(V,U)=(E-\Gamma -H_0 -V-\Sigma^\mathrm{Int})^{-1}$ is calculated, where $\Sigma^\mathrm{Int}$ is the second-order expansion of the electron interactions obtained self-consistently, as illustrated in the secondary loop in Fig.1(A) of the main text using the diagram shown in Fig.~\ref{Fig:sopt} herein. (v) With the fully dressed Green function, we calculate the cluster density of states $\rho^c=-\frac{1}{\pi} \Im G^c$ and average over many configurations to obtain the momentum-resolved, non-self-averaged typical DoS~\cite{PhysRevB.89.081107,PhysRevB.92.201114}
\begin{equation}
    \rho_{t}^c(\mathbf{K}) = \langle \rho^c_i\rangle_\mathrm{geom} \left\langle \frac{\rho^c(\mathbf{K})}{\frac{1}{N_c} \sum_{i}\rho_{i}^c}\right\rangle_\mathrm{arit},
    \label{e.2}
\end{equation}
where $\langle \rho^c_i\rangle_\mathrm{geom}=\exp\,\langle\ln\rho_i\rangle_\mathrm{arit}$. The second factor ensures that non-local fluctuations are captured while the system remains in a typical medium.  (vi) A cluster typical Green function $G_t^c(\mathbf{K})$ is then obtained from the Kramer-Kronig transform of Eq.\,\ref{e.2} as $ G_{t}^c(\mathbf{K},E)=\int d E' \displaystyle \rho_{typ}^c(\mathbf{K},E')/(E - E')$ and used to calculate the coarse-grained Green function
\begin{equation}
    \bar{G} (\mathbf{K}) = \frac{N_c}{N} \sum_{\mathbf{\tilde k}}  \bigg [\displaystyle G^c_{t} (\mathbf{K})^{-1} + \Gamma (\mathbf{K}) - H_0(\mathbf{k}) + \bar{H}_0(\mathbf{K}) + \mu \bigg ]^{-1}
    \label{e.3}
\end{equation}
where the overbar depicts cluster coarse-graining and $\mu$ is the Fermi level, which we obtained in step (iii).  (vii) A new hybridization function is obtained using linear mixing
\begin{equation}
    \Gamma_\mathrm{n}(\mathbf{K}) = (1-\zeta)\Gamma_\mathrm{o}(\mathbf{K})+\zeta\left[(G^c)^{-1}-\bar G^{-1}\right],
    \label{e.4}
\end{equation}
where $\Gamma_\mathrm{n}$ ($\Gamma_\mathrm{o}$) refers to the new (old) hybridization function and $\zeta$ is a mixing parameter.  (viii) $\Gamma(\mathbf{K})\equiv\Gamma_\mathrm{n}(\mathbf{K})$ is then used in (ii) to repeat the above procedure until the hybridization function converges to the desired accuracy. When this happens, $G_t^c\equiv\bar{G} $ within the computational error.

\section{Convergence of the TMDCA@DFT with Cluster Size}
We show in Fig.\,\ref{f.2} the typical density of states (TDoS) for various cluster sites $N_c=1$, 4, 8, and 16.  Figure~\ref{f.2}(A) depicts the TDoS in the noninteracting limit at the boron vacancy concentrations $\delta=5$ and $10\%$, while Fig.\,\ref{f.2}(B) shows the the TDoS at the boron vacancy concentrations $\delta=0.5$ and $1.0\%$ at the electron interaction strength $U=1.2$ eV.
\begin{figure}[b!]
    \includegraphics[trim = 0mm 0mm 0mm 0mm,scale=0.45,keepaspectratio,clip=true]{TDos_Nc.eps}
    \caption{The TMDCA@DFT results of $h$-BN obtained for $N_c$ of 1, 4, 8, and 16 at two boron vacancy concentrations $\delta$ for (A) noninteracting and (B) interacting cases.  Note that the electrochemical potential has been shifted to allow direct comparison of the spectra.}
    \label{f.2}
\end{figure}
The essence of this plot is to benchmark the convergence of the TMDCA@DFT method.  In both cases, the TDoS systematically converges as $N_c$ increases throughout almost the entire energy spectrum. This rather fast convergence ensures that with relatively small $N_c$, accurate results can be obtained, enabling simulation of many-orbital material as in our hexagonal boron nitride example. Even the $N_c=1$ cluster gives qualitative results, which could be important for computations of larger systems.  Note also that in the interacting limit, there are additional scattering processes leading to the emergence of deep-level impurity states.

This fast convergence also makes it difficult to carry out standard finite-size-scaling of the TDoS. However, since the TDoS profile at almost all energy is converged as a function of $N_c$, our calculations remain in the thermodynamic limit to within the computation error of our method. The TMDCA self-consistently calculates the crucial, critical order parameter describing the disorder-induced formation of electronic bound state, by focusing on the typical (most probable) local density of the states. The validity and even quantitative accuracy of the TMDCA applied to model systems have been established in a series papers. The TMDCA is benchmarked against exact numerical methods, e.g. Kernel polynomial method, Transfer matrix method, Exact diagonalization (see for, e.g., Refs.~\cite{PhysRevB.92.014209,*PhysRevB.90.094208,*PhysRevB.92.205111} for the Anderson model. Also, the TMDCA is benchmarked against the continuous time quantum Monte-Carlo simulations for interacting disordered electron systems using the Anderson-Hubbard model~\cite{PhysRevB.92.201114}.

\ifx\mcitethebibliography\mciteundefinedmacro
\PackageError{unsrtM.bst}{mciteplus.sty has not been loaded}
{This bibstyle requires the use of the mciteplus package.}\fi
\begin{mcitethebibliography}{10}

\bibitem{PhysRevB.56.12847}
N.~Marzari and D.~Vanderbilt, Phys. Rev. B {\bf 56}, 12847 (1997)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{Mostofi2008685}
A.~A. Mostofi, J.~R. Yates, Y.-S. Lee, I.~Souza, D.~Vanderbilt, and N.~Marzari,
  Comp. Phys. Commun. {\bf 178}(9), 685 (2008)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.63.125102}
M.~Jarrell and H.~R. Krishnamurthy, Phys. Rev. B {\bf 63}, 125102 (2001)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.61.12739}
M.~H. Hettler, M.~Mukherjee, M.~Jarrell, and H.~R. Krishnamurthy, Phys. Rev. B
  {\bf 61}, 12739 (2000)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{Jarrell01}
M.~Jarrell, T.~Maier, C.~Huscroft, and S.~Moukouri, Phys. Rev. B {\bf 64},
  195130 (2001)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.92.201114}
C.~E. Ekuma, S.-X. Yang, H.~Terletska, K.-M. Tam, N.~S. Vidhyadhiraja,
  J.~Moreno, and M.~Jarrell, Phys. Rev. B {\bf 92}, 201114 (2015)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.89.081107}
C.~E. Ekuma, H.~Terletska, K.-M. Tam, Z.-Y. Meng, J.~Moreno, and M.~Jarrell,
  Phys. Rev. B {\bf 89}, 081107 (2014)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.92.014209}
C.~E. Ekuma, C.~Moore, H.~Terletska, K.-M. Tam, J.~Moreno, M.~Jarrell, and
  N.~S. Vidhyadhiraja, Phys. Rev. B {\bf 92}, 014209 (2015)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.90.094208}
H.~Terletska, C.~E. Ekuma, C.~Moore, K.-M. Tam, J.~Moreno, and M.~Jarrell,
  Phys. Rev. B {\bf 90}, 094208 (2014)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem{PhysRevB.92.205111}
Y.~Zhang, H.~Terletska, C.~Moore, C.~Ekuma, K.-M. Tam, T.~Berlijn, W.~Ku,
  J.~Moreno, and M.~Jarrell, Phys. Rev. B {\bf 92}, 205111 (2015)\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\end{mcitethebibliography}

\end{widetext}

\end{document}