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

The existence of a non-luminous component of matter in the universe, dark matter (DM), is
supported by compelling astrophysical evidence~\cite{Zwicky:1933gu, DMreview}.
The abundance of dark matter has been precisely determined
from global
fits of cosmological parameters to a variety of observations~\cite{Komatsu:2010fb,Planck:2015fie}.
Nevertheless, the nature of DM remains largely unknown.

In this paper, models where the DM candidate is a
weakly interacting massive particle (WIMP)~\cite{Steigman:1984ac} are considered.
In proton--proton ($pp$) collisions at the Large Hadron Collider (LHC), pair-produced WIMP DM
does not interact with the detector
and it can be detected only if produced in association with Standard Model (SM) particles,
leading to signatures with missing transverse momentum.
A wide range of experimental searches
are focused on WIMP candidates at the LHC~
\cite{EXOT-2018-06,CMS-EXO-19-003,CMS-EXO-19-003,CMS-EXO-19-003,CMS-EXO-19-003,SUSY-2018-34,EXOT-2018-63,SUSY-2016-18,SUSY-2016-18,ATLAS:2021shl,EXOT-2017-32}.
All recent searches are based on simplified benchmark models
documented in the LPCC Dark Matter Working Group whitepapers~\cite{Boveia:2016mrp,Albert:2017onk,LHCDarkMatterWorkingGroup:2018ufk,Abercrombie:2015wmb}.
Benchmark models are chosen to have a minimal number of additional parameters relative to the SM.
This paper focuses on simplified models characterised by the introduction of a spin-0 particle mediator~\cite{Blanke:2020bsf,Westhoff:2019pzy,LHCDarkMatterWorkingGroup:2018ufk,Haisch:2018bby,Haisch:2016gry,Haisch:2015ioa,Abercrombie:2015wmb,Abdallah:2015ter,Buckley:2014fba}. In this case, fermionic DM particle pairs
are produced through the exchange of a colour-neutral scalar or
pseudoscalar mediator (denoted by $\phi$ or $a$, respectively).
In the Minimal Flavour
Violation~\cite{DAmbrosio:2002vsn} assumption,
the interaction between any new neutral spin-0
state and SM quarks is proportional to the fermion masses via
Yukawa-type couplings and it is also treated as a free parameter of the model by means of a multiplicative factor $g_q$.
Following Ref.~\cite{Buckley:2014fba}, couplings to leptons and $W$/$Z$ bosons, as well as explicit $\phi$--$H$ or $a$--$H$ couplings of dimension four to the SM Higgs boson,
are set to zero. The coupling of the mediator to the dark sector, $g_{\chi}$, is not assumed to be proportional to the mass of
the DM candidates and is treated simply as a free parameter.

The dominant production modes for such colour-neutral mediators are
loop-induced gluon-fusion and associated production of the mediator with a top quark pair.
Figures~\ref{fig:diagrams}(a) and~\ref{fig:diagrams}(b) show the two dominant production diagrams for the associated production of the mediator with a top quark pair.
As discussed in Refs.~\cite{Haisch:2016gry,Backovic:2015soa}, the diagram in Figure~\ref{fig:diagrams}(b) can be interpreted as the radiation of the mediator from a top quark. This process is enhanced for small energies and leads to a production cross section for scalar mediators about one order of magnitude larger than the one for pseudoscalar mediators for masses below the top quark threshold.
The associated production of spin-0 mediators with a single top quark
also has a sizeable, albeit non-dominant, cross section~\cite{Pinna:2017tay,Pani:2017qyd,Haisch:2018bby}, especially for higher-mass mediators.
The primary signal of interest in this paper is the associated production of a mediator particle with a pair of top quarks (DM+$\ttbar$),
although sensitivity is also retained in the case of single top quark production (collectively referred to as DM+$t$).
The relevant processes for DM+$t$ are shown in Figures~\ref{fig:diagrams}(c) and~\ref{fig:diagrams}(d).
The relative contribution of the two processes depends on the parameter space that is considered~\cite{Pinna:2017tay}.

\begin{figure}[ht!b]
\begin{center}
\begin{subfigure}[t]{0.22\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_01a.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}[t]{0.22\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_01b.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}[t]{0.22\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_01c.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}[t]{0.22\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_01d.pdf}
\caption{}
\end{subfigure}
\caption{Representative diagrams for spin-0 mediator associated production with (a) and (b) a top quark pair (DM+$\ttbar$),
(c) a single top quark and a $W$ boson (DM+$tW$) or (d) a single top quark and one (or more) jet(s) (DM+$tj$).
}
\label{fig:diagrams}
\end{center}
\end{figure}

This paper presents a statistical combination of three searches targeting events with two top quarks and invisible particles,
considering
either zero- (tt0L~\cite{SUSY-2018-12}), one- (tt1L~\cite{SUSY-2018-07}) or two-lepton (tt2L~\cite{SUSY-2018-08}) final states,
using 139~fb$^{-1}$ of $pp$ collisions data recorded by the ATLAS detector~\cite{PERF-2007-01, ATLAS-TDR-2010-19} at a centre-of-mass energy of $\sqrt{s}
= 13~\TeV$.
The tt0L analysis is extended and improved relative to that in Ref.~\cite{SUSY-2018-12},
benefitting from the Run-2 improvements in the trigger selection of jets containing
$b$-hadrons ($b$-jets) in ATLAS~\cite{ATLAS:2021piz} that increase the sensitivity in models with low-mass mediators.
The three analyses select independent (orthogonal) datasets which can be statistically combined to boost the sensitivity
because they seek signal events in similar parts of the parameter space.
Searches specifically targeting the associated
production of a single top quark and invisible particles \cite{ATLAS:2020yzc} are not orthogonal to
the searches presented here. They are therefore not included and their combination is left as a future development.

The data are found to be in agreement with the SM background prediction in all selections considered in
this paper. The statistical combination is used
to set 95\% confidence level (CL) constraints on the simplified DM
models.
The specific case where the mediator corresponds to the SM $125~\GeV$ Higgs boson~\cite{HIGG-2012-27}
is also considered when interpreting the results presented in this paper. It is referred to as the \hinv\ interpretation in the following. In
the SM, the invisible Higgs boson branching ratio, \brinv,  is 0.12\% from $H \rightarrow ZZ \rightarrow  4\nu$ decays~\cite{deFlorian:2016spz},
and higher branching ratios to invisible particles are predicted by Higgs--dark-matter portal models~\cite{Shrock:1982kd,Choudhury:1993hv,ArkaniHamed:1998vp,Eboli:2000ze,Godbole:2003it,Antoniadis:2004se,Davoudiasl:2004aj,Patt:2006fw,Kanemura:2010sh,Djouadi:2011aa,Ghosh:2012ep,Belanger:2013kya,Curtin:2013fra,Zaazoua:2021xls}.
Results on the invisible branching ratio obtained from the statistical combination of previous ATLAS searches conducted with the Run-1 and partial Run-2 dataset reported an observed (expected) limit of $0.26$ ($0.17$)~\cite{HIGG-2018-54}, while CMS reported $0.19$ ($0.15$)~\cite{CMS-HIG-17-023} at $95\%$ CL.
Recent updates by the ATLAS and CMS Collaborations in the vector-boson-fusion channel report improved observed (expected) upper limits of $0.145~(0.103)$ \cite{ATLAS:2022yvh} and $0.18~(0.10)$~\cite{PhysRevD.105.092007}, respectively, using the full Run-2 dataset.

The paper is structured as follows. The experimental dataset and the simulated event samples are presented in Section~\ref{sec:samples}, the new and previously published searches are introduced in Section~\ref{sec:analysis}, and their statistical combination, including the treatment of correlated systematic uncertainties is discussed in Section~\ref{sec:comb}. Results in terms of exclusion limits are presented in Section~\ref{sec:results}, followed by the conclusions in Section~\ref{sec:conclusion}. Finally, additional details of the extended 0-lepton channel analysis are discussed in Appendix~\ref{sec:lowmetapp}.

\section{Data and simulated event samples}
\label{sec:samples}

The dataset used in the analyses described in this paper consists of
$pp$ collision data recorded by the ATLAS detector~\cite{PERF-2007-01,ATLAS-TDR-2010-19} at a centre-of-mass energy of $\sqrt{s}
= 13~\TeV$ with stable beam conditions.
The ATLAS detector is a multipurpose particle detector with a forward--backward symmetric cylindrical geometry and nearly full coverage in solid angle.\footnote{
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 upward.
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)$.
The angular distance between two objects in $\eta$--$\phi$ space is defined by $\Delta R \equiv \sqrt{(\Delta\eta)^2+(\Delta\phi)^2}$.}
It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2~T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer.
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.

All collision events considered in this paper are required to have at
least one reconstructed interaction vertex with a minimum of two
associated tracks, each having $\pt>500$~\MeV.
In events with multiple vertices, the one with the highest sum of
squared transverse momenta of associated tracks is
chosen as the primary vertex~\cite{ATL-PHYS-PUB-2015-026}.
Electrons~\cite{EGAM-2018-01}, muons~\cite{PERF-2015-10}, $\tau$-leptons~\cite{PERF-2013-06} and jets are reconstructed by combining the
signals from the different components of the ATLAS
detector.
Small-radius jets are reconstructed from energy deposits in the calorimeters using
the anti-$k_t$ jet
algorithm~\cite{paper:antikt,Cacciari:2011ma} with a radius parameter
of $R = 0.4$. Reclustered large-radius jets \cite{Nachman:2014kla,ATLAS-CONF-2017-062} are reconstructed
with the same algorithm, using a radius parameter of $R=1.2$ and with small-radius jets as input, unless otherwise specified.
Multivariate algorithms  are used to
identify
small-$R$ jets with $\pt >20\;\GeV$ containing $b$-hadrons
($b$-jets)~\cite{PERF-2012-04,ATL-PHYS-PUB-2016-012}.
This is referred to as $b$-tagging.
The missing transverse momentum $\vec{p}_{\text{T}}^{\text{miss}}$ (with
magnitude $\met$) is calculated from the negative vector sum of the transverse
momenta  of electrons, muons
and jet candidates, and an additional `soft term'~\cite{PERF-2016-07}
which includes tracks found in the tracking system
that originate from the primary vertex but are not
associated with any reconstructed object.

Depending on the analysis channel, events are selected by lepton triggers~\cite{TRIG-2018-05,TRIG-2018-01}, $\met$ triggers~\cite{TRIG-2019-01} or $b$-jet triggers~\cite{ATLAS:2021piz}. The uncertainty in the combined integrated luminosity is 1.7\%~\cite{ATLAS-CONF-2019-021}, obtained using the LUCID-2 detector~\cite{LUCID2} for the primary luminosity measurements.
Events accepted by lepton and $\met$ triggers are required to meet the standard ATLAS data-quality assessment
criteria~\cite{DAPR-2018-01} to ensure that all subdetector systems were functioning normally. This leads to an integrated luminosity of
3.2~fb$^{-1}$, 33.0~fb$^{-1}$, 44.3~fb$^{-1}$ and 58.5~fb$^{-1}$ in 2015, 2016, 2017 and 2018 data-taking, respectively,
with a total integrated luminosity of 139.0~fb$^{-1}$.
Events accepted by $b$-jet triggers are required to meet additional criteria ensuring that the online beam-spot position measurement is valid,
which leads to the exclusion of the 2015 data and to a reduced integrated luminosity of 24.6~fb$^{-1}$, 43.7~fb$^{-1}$
and 57.7~fb$^{-1}$ in 2016, 2017 and 2018 data-taking, respectively, with a total integrated luminosity of 126.0~fb$^{-1}$.
The $b$-jet trigger chains considered in this paper require at least four jets, among which two jets are $b$-tagged
by the online version of the $b$-tagging boosted decision tree algorithm~\cite{ATL-PHYS-PUB-2016-012}.

Dedicated Monte Carlo (MC) simulated event samples are used to aid in the
estimation of the background from SM processes and to model the
dark matter or invisible Higgs boson signal. All simulated events were processed through a
simulation of the ATLAS detector and its response~\cite{SOFT-2010-01,Agostinelli:2002hh}.
The simulated events are reconstructed with the same algorithms as used for data.
They contain a realistic modelling of additional $pp$ collisions in the same or neighbouring bunch crossings (pile-up), obtained
by overlaying minimum-bias events simulated using the soft QCD processes of \PYTHIA[8.186]~\cite{Sjostrand:2014zea,Sjostrand:2007gs}
with the \NNPDF[2.3lo] set of parton distribution functions (PDFs)~\cite{Ball:2012cx} and the
A3~\cite{ATL-PHYS-PUB-2016-017} set of tuned parameters.

Standard Model processes producing one or two top quarks
were modelled using \POWHEGBOX[v2]~\cite{Frixione:2007nw,Nason:2004rx,Frixione:2007vw,Alioli:2010xd}
and normalised to cross sections calculated to next-to-next-to-leading-order (NNLO) plus next-to-next-to-leading-logarithm (NNLL)~\cite{Czakon:2011xx} accuracy in QCD.
Processes involving the production of one (two) vector boson(s) were simulated using \SHERPA[2.2.1]~(2.2.2)~\cite{Gleisberg:2008ta,Gleisberg:2008fv,Cascioli:2011va,Schumann:2007mg,Hoeche:2012yf} and
normalised to cross sections at NNLO~\cite{Catani:2007vq} (NLO) accuracy in QCD. The $t\bar{t}{+}V$ ($V=W,Z$), $tZ$ and $tWZ$ processes were modelled by
\MGNLO[2.3.3]~\cite{Alwall:2014hca} and normalised to NLO cross-section accuracy in QCD~\cite{Alwall:2014hca,deFlorian:2016spz,Alwall:2014hca}.
Lastly, $t\bar{t}H$ background events, with visible Higgs boson decay modes, were generated by \POWHEGBOX[v2].

Signal samples for dark matter produced in association with a \ttbar
pair (DM+\ttbar) were generated using a leading-order (LO) matrix element, with up to one
extra parton, in the \MGNLO[2.6.7]~\cite{Alwall:2014hca}
generator interfaced to \PYTHIA[8.244] and using the CKKW-L merging algorithm~\cite{Lonnblad:2001iq}.
The five-flavour scheme \NNPDF[3.0nlo]~\cite{Ball:2014uwa} PDF set was used.
The top quark decay was simulated using \MADSPIN~\cite{Artoisenet:2012st}.
Signal cross sections for this process were calculated to NLO QCD accuracy using the same version of
\MADGRAPH, as suggested in Ref.~\cite{Backovic:2015soa}.
Models with a DM particle mass of 1~\GeV\ and $g_q = g_{\chi} = 1$ are considered, while the masses of the mediators range from 10~\GeV\ to 400~\GeV.
The typical scaling factors from LO to NLO in these calculations range between 1.25 and 1.35,
depending on the mediator mass and whether it is a scalar or a pseudoscalar boson.
Signal samples for dark matter produced in association with a single top
quark were generated using the same settings as for the DM+\ttbar\
samples. For these signal models, the $tW$ and the $tj$ processes were generated separately.
Each one was normalised to the LO cross section predicted by the model and then
the samples were combined. No extra partons were generated from the matrix element in this case.

Signal samples modelling Higgs boson production in association with a \ttbar pair \((t\bar{t}H)\) were generated using
the \POWHEGBOX[v2]~\cite{Hartanto:2015uka}
generator at NLO with the \NNPDF[3.0nlo]~\cite{Ball:2014uwa} PDF set.
The \(tWH\) signal samples were produced with \MGNLO[2.6.2] in the five-flavour scheme with the \NNPDF[3.0nnlo] PDF set. The top quark and $W$ boson decays were handled by \MADSPIN to account for spin correlations among the decay products. The overlap of the \(tWH\) process with \(t\bar{t}H\) at NLO was removed by using a diagram removal technique~\cite{Frixione:2008yi,Demartin:2016axk}.
Contributions arising from $tH$ production in the \(t\)-channel and \(s\)-channel are below 0.1\% and are therefore neglected.
In these samples, the Higgs boson decays via $ZZ^*$ into neutrinos, and events are normalised using
the total cross section at NLO QCD and electroweak accuracy recommended by the LHC Higgs Cross
Section Working Group~\cite{deFlorian:2016spz,Demartin:2015uha,Demartin:2016axk}.

For the SM background samples, except those generated using \SHERPA,
the \EVTGEN[1.2.0]~\cite{Lange:2001uf} program was used to simulate
the properties of the $b$- and $c$-hadron decays, and \PYTHIA[8.186] with the A14
tune~\cite{ATL-PHYS-PUB-2014-021} was used for the parton showering,
hadronisation model and underlying event. For the \SHERPA\ generator, the default \SHERPA~\cite{Gleisberg:2008ta,Gleisberg:2008fv,Cascioli:2011va,Schumann:2007mg,Hoeche:2012yf} configuration recommended
by its authors was used.

The modelling of the response of the various ATLAS subdetectors was performed using $\GEANT$~\cite{Agostinelli:2002hh}
for all the background MC samples, while all signal MC samples were simulated using a faster simulation based on a parameterisation of the calorimeter response and $\GEANT$ for the other detector systems~\cite{SOFT-2010-01}.

\section{Experimental signatures and analysis strategy}
\label{sec:analysis}

Each of the analyses requires the presence of at least one $b$-tagged jet and
\met in the event.
To further suppress events
where the \met\ originates from mismeasurements, the ratio of the \met to
its resolution is used to construct the event's missing transverse
momentum significance~\cite{ATLAS-CONF-2018-038}, $\mathcal{S}$.

The tt0L analysis targets decay chains with no leptons in the final state. The tt0L analysis published
in Ref.~\cite{SUSY-2018-12} focuses on a set of signal regions (SR), referred to as \noLhigh\ in this paper, defined
so as to obtain maximal significance for the targeted signal events. It
uses the $\met$ triggers to select events with large missing transverse momentum
and requires at least one highly
energetic, hadronically decaying top quark candidate. An additional set of SRs,
referred to as \noLlow, is presented for the first time in this paper and extends the previously published results by relying on a combination of $\met$ and $b$-tagged jet triggers to retain events with lower-momentum jets that fail one of the
\noLhigh\ analysis criteria.

The two other analyses considered in this paper target leptonic final states. The \oneL\ analysis~\cite{SUSY-2018-07} selects events with exactly one
lepton ($e$ or $\mu$) and it is based on the $\met$ triggers. The \twoL\ analysis~\cite{SUSY-2018-08} targets events with exactly two opposite-charge leptons ($e$ or $\mu$) in the final states selected with dilepton triggers.

Common event-quality criteria and object reconstruction and identification definitions
for leptons, jets, $b$-tagged jets and \met
are applied in all analyses considered in the combination.
Minimum \pt requirements of $4.5\;\GeV$ and $20\;\GeV$ are applied to leptons and jets, respectively.
The lepton multiplicity requirements guarantee that the 0-lepton, 1-lepton and 2-lepton channels are by construction
non-overlapping.
Further kinematics, isolation and identification requirements are applied to each object in the individual analysis channels and
optimised for each topology. These details, together with the details of the reconstruction algorithms, quality requirements and efficiencies can be found in
Refs.~\cite{SUSY-2018-12,SUSY-2018-07,SUSY-2018-08} and in Appendix~\ref{app:var}.
For all analyses in this paper, background-enriched selections (control regions, CR) are defined so as to allow the data to aid in the
estimation of the dominant SM backgrounds, and validation regions (VR) are
used to verify the robustness of these estimates.
A dedicated background estimation strategy was developed for each channel~\cite{SUSY-2018-12,SUSY-2018-07,SUSY-2018-08},
using independent control regions for all dominant SM processes.

Observed and expected
event yields in the signal and control regions are used in Poisson
probability functions to build
likelihood functions, which are combined in a profile likelihood fit. A
profile likelihood ratio is employed in the $\mathrm{CL_s}$ method to exclude at 95\% CL the
signal-plus-background hypothesis for the signal models considered~\cite{Baak:2014wma,Cowan:2010js,Junk:1999kv,Read:2002hq}.

\subsection{Description of analysis channels}
\label{sec:analyses_detail}

\paragraph{The \noL\ analysis}

The experimental signature targeted in this channel consists of at least four jets, two of which are \btagged, and large missing transverse momentum ($\met > 160\;\GeV$).  Events with electrons, muons or $\tau$-leptons are rejected.
The \noLhigh\ selection consists of events with high missing transverse momentum (\met trigger, $\met > 250\;\GeV, \mathcal{S}> 14$) where at least one of the reconstructed large-radius jets ($R=1.2$) has a mass consistent with one produced from a boosted hadronic top quark decay.
This selection corresponds to signal region selections SRA and SRB in Ref.~\cite{SUSY-2018-12}, which were originally optimised to be sensitive to
high-mass supersymmetric partners of the top quark and are divided into three categories
(TT, TW and T0), depending on whether the subleading large-radius jet ($R=1.2$)
has an invariant mass consistent with a top quark, a $W$ boson or
neither. Given the signature and kinematic similarity between the signal considered in this paper and
top squark pairs decaying into a top quark and a neutralino, these regions can also be used to
constrain DM+$\ttbar$ models.
SRA and SRB are orthogonal because of a requirement on the
$\chi^2$-based stransverse mass variable, $m_{\mathrm T2, \chi^2}$, and are statistically combined.
This stransverse mass~\cite{Lester:1999tx,Barr:2003rg} variable is designed to reconstruct the
mass of two heavy particles produced in an event and decaying symmetrically into
an invisible particle and a top quark, as expected in supersymmetric top quark topologies. It relies on a $\chi^2$-based method to identify the hadronically decaying top
quark candidates.
The background estimation in SRA and SRB is aided by means of
dedicated control regions for all dominant SM
processes: $t\bar{t}Z$, $Z$+jets, \ttbar, single top quark in the $tW$ channel, and $W$+jets. Event
yields with updated jet calibrations~\cite{JETM-2018-05} were recalculated with
respect to Ref.~\cite{SUSY-2018-12} and are presented in
Figure~\ref{fig:srplots_0L}. The yield change due to the new calibration is between
6\% and 15\% and is reflected in the signal predictions as well. This is because the new calibration decreases the
contribution of events with artificially increased missing transverse momentum, due to mismeasured jets in the analysis, by reducing the
tails of the  missing transverse momentum significance distribution.
This translates to a slightly higher signal-to-background ratio in the
\noLhigh\ SRs.

The \noLlow\ selection is newly added in this paper with the aim of improving the sensitivity of the 0-lepton channel
for dark matter models by selecting final states with lower missing transverse momentum
and/or lower-momentum objects.
The \noLlow\ selection is therefore designed to maximise the sensitivity to
DM+\ttbar signals with low mediator masses
($m(\phi),m(a)<100~\GeV$). Events are selected by a combination of $\met$ triggers and $b$-jet triggers.
Events selected with $\met$ triggers must fulfil $\met > 250\;\GeV$
to ensure the triggers are fully efficient. In addition, they are required to have either $\mathcal{S}<14$
or no large-radius jets consistent with highly energetic top quark candidates to ensure orthogonality with the \noLhigh\ selections.
Events selected with $b$-jet triggers must fulfil the requirement that the missing transverse momentum be between 160 and 250~\GeV.
The lower bound suppresses the multi-jet background contamination, while
the upper bound ensures orthogonality with \noLhigh. The online $b$-tagged jet
candidates must match the offline $b$-tagged jet candidates within a cone of $\Delta R = 0.2$.

Three signal regions are defined, SR0X,
SRWX and SRTX, according to the mass of the heaviest large-radius jet, which is used to infer the \pt\ of the most energetic
top quark in the event. SR0X requires
that no large-radius jets (0) are present, while in SRWX and SRTX the
mass of the heaviest large-radius jet has to be respectively lower
or higher than 130~\GeV, i.e.\ in the neighbourhood of the $W$ boson
(W) or top quark (T) mass.

An optimal categorisation of the events in which a top quark is produced
is particularly important in the \noLlow\ analysis because of
the high top quark background rate and its large variation with respect to the
top quark transverse momentum. In order to maximise the correlation of the large-radius jet mass
with the \pt of the most energetic top quark,  jet reconstruction in the \noLlow\ analysis uses a
smaller radius parameter value ($R=1.0$) than in the \noLhigh\ analysis.
This is  associated with newly
introduced lower selections on the large-radius jet transverse
momentum ($\pt>200\;\GeV$), and invariant mass
($m_{\textrm{large-radius jet}} > 40\;\GeV$).
A requirement rejecting high-mass, large-radius jets with
radius parameter $R=1.2$, as defined in the \noLhigh\ analysis, is also applied to ensure orthogonality between the two analyses.

The final states targeted in the \noLlow\ analysis are less energetic than in the \noLhigh\ analysis, such that the decay
products of the top quarks are expected to be less collimated. Hence no subleading large-radius jet is required (X).

The full list of requirements for the three \noLlow\ signal regions is reported in
Table~\ref{tab:lowmet_sr}.
Three angular separation variables are used: $\Delta\phi_{\min}(\vec{p}_{\text{T,1-4}},\vec{p}_{\text{T}}^{\text{miss}})$
to reduce the contamination from multi-jet events,
\drbb\ to reduce the contamination from $b$-jets from gluon splitting as present in $Z$+jets events,
and $\Delta R_{\min}$(large-radius jet, $b$-tagged jets) to enhance the fraction of events with hadronic top quark decays, as in the signal.
The SM background originating from the top quark is reduced by using the \coshmax\ variable. This variable aims to identify events
where a leptonically decaying $W$ boson is the source of all \met\ in the event because the lepton is not reconstructed.
The hyperbolic cosine of the pseudorapidity difference between the missed $W$ boson candidate and each of the two $b$-tagged jets selected in the event
is estimated, and the \coshmax\ variable is defined to be the larger of the two values:
\begin{align}
\coshmax = \max\{\cosh\Delta\eta^{1}_{W,b}, \cosh\Delta\eta^{2}_{W,b}\}. \nonumber
\end{align}
The value of $\cosh\Delta\eta_{W,b}$
is estimated by solving the kinematics of the top quark decay,
assuming $m_{W} \ll \pt^{W}\cdot\cosh\eta_{W} \sim \met\cdot\cosh\eta_{W}$, given that $\met > 160~\GeV$ and $\cosh\eta_{W}\ge1$.
For events with \ttbar\ pairs decaying fully hadronically, the approximation $\pt^{W} \sim \met$ is not valid, resulting in unphysical values, e.g.\ $\coshmax < 1$,  while events with high \coshmax\ values are likely to contain a top quark decaying leptonically and are excluded from the signal regions.
Since the signal events have two hadronically decaying top quarks, while most of the backgrounds have one or none, a $\chi^2$-based variable, \chitwotop,
is used to identify events with a hadronically decaying top quark pair (Appendix~\ref{app:var}).
The \chitwotop\ variable is constructed using up to six jets in the
event and using the $b$-tagging classification when appropriate to
reconstruct the two $W$ bosons and the two top quarks in the
event. Invariant mass constraint terms for  these particles are
used to build the $\chi^2$. All possible jet assignments are tested and the one with the lowest $\chi^2$ is chosen.
The event distributions for the \coshmax\ and the \chitwotop\ variables are each shown in Figure~\ref{fig:key_discriminants}
after applying all the SR0X selection requirements except those on the variable presented in the distribution.
The transverse momentum of the \ttbar system ($\pt^{\ttbar}$) constructed with the \chitwotop\ method is compared with the \met.
In signal events these two quantities are expected to have similar values, so events with $\pt^{\ttbar}/\met$ values outside a window containing unity are excluded.

\begin{figure}[p]
\begin{center}
\begin{subfigure}[t]{0.96\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_02a.pdf}
\caption{}
\label{fig:srplots_0L}
\end{subfigure}

\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_02b.pdf}
\caption{}
\label{fig:srplots_1L}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_02c.pdf}
\caption{}
\label{fig:srplots_2L}
\end{subfigure}
\caption{(a) Post-fit signal region yields for the \noLhigh\ and the \noLlow\ analyses.
The bottom panel shows the statistical significance~\cite{ATL-PHYS-PUB-2020-025} of the difference
between the SM prediction and the observed data in each region.
The definitions of SRAT0, SRATW, SRATT, SRBT0, SRBTW and SRBTT can be found in Ref.~\cite{SUSY-2018-12}.
Representative post-fit distributions are presented for (b) the \oneL\ and (c) the \twoL\ analyses: each bin of these distributions,
starting from the red arrow in (c), corresponds to a single SR included in the fit.
In the \noLlow\ analysis, `\ttbar\ (other)' represents \ttbar\ events without extra jets or events with extra light-flavour jets.
In the \twoL\ analysis, `FNP' includes the contribution from fake/non-prompt lepton background arising from jets
(mainly $\pi/K$, heavy-flavour hadron decays and photon conversion) misidentified as leptons, estimated in a purely data-driven way.
`Other' includes contributions from $\ttbar{+}W$, $tZ$ and $tWZ$ processes, and also \ttbar (semileptonic) for the \oneL\ analysis.
The total uncertainty in the SM expectation is represented with hatched bands and the expected distributions
for selected signal models are shown as dashed lines. }
\label{fig:srplots}
\end{center}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_03a.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_03b.pdf}
\caption{}
\end{subfigure}
\caption{Distributions of (a) \coshmax and (b) \chitwotop in SR0X events passing all the SR requirements
except those on the variable being presented (which are indicated by the arrows).
The contributions from all SM backgrounds are shown after the simultaneous profile likelihood  fit to all \noLlow\ CRs,
with the hatched bands representing the total uncertainty.
The category `\ttbar\ (other)' represents \ttbar\ events without extra jets or events with extra light-flavour jets.
`Other' includes contributions from $\ttbar{+}W$, $tZ$ and $tWZ$ processes.
The expected distributions for selected signal models are shown as dashed lines.
The underflow (overflow) events are included in the first (last) bin.
The bottom panels show the ratio of the observed data to the total SM background prediction,
with the hatched area representing the total uncertainty in the background prediction.}
\label{fig:key_discriminants}
\end{center}
\end{figure}

\begin{table}[htbp]
\centering
\scriptsize
\caption{Selection criteria for the signal regions used in the \noLlow\ analysis.}
\label{tab:lowmet_sr}

\begin{tabular}{|>{\centering}p{0.3\textwidth-2\tabcolsep}|>{\centering}p{0.23\textwidth-2\tabcolsep}|>{\centering}p{0.23\textwidth-2\tabcolsep}|>{\centering\arraybackslash}p{0.23\textwidth-2\tabcolsep}|}
\hline
\textbf{Variables} & \textbf{SR0X} & \textbf{SRWX} & \textbf{SRTX} \\
\hline
$N_{\textrm{lepton}}$ & \multicolumn{3}{c|}{$=$ 0} \\
\hline
Orthogonalisation & \multicolumn{3}{c|}{$\met < 250~\GeV$ \hspace{0.01\textwidth} or \hspace{0.01\textwidth} $\mathcal{S} < 14$ \hspace{0.01\textwidth} or \hspace{0.01\textwidth} $m_{\textrm{large-radius jet}}^{R=1.2} < 120~\GeV$} \\
\hline
\multirow{2}{*}[0.6mm]{\met[GeV]} & \multicolumn{3}{c|}{$>$ 160} \vspace{-1.2mm}\\
& \multicolumn{3}{c|}{$<$ 250, when passing $b$-jet triggers} \\
\hline
$\mathcal{S}$ & \multicolumn{3}{c|}{$>$ 10} \\
\hline
$\Delta\phi_{\min}(\vec{p}_{\text{T,1-4}},\vec{p}_{\text{T}}^{\text{miss}})$ & $>$ 1.0 & \multicolumn{2}{c|}{$>$ 0.5} \\
\hline
\drbb & \multicolumn{3}{c|}{$>$ 1.2} \\
\hline
$N_{\textrm{large-radius jet}}$ & $=$ 0 & \multicolumn{2}{c|}{$>$ 0} \\
\hline
$m_{\textrm{large-radius jet}}$ [GeV] & --- & (40, 130) & $\geq$ 130 \\
\hline
$\Delta R_{\min}$(large-radius jet, $b$-tagged jets) & \multicolumn{2}{c|}{---} & $<$ 1.2 \\
\hline
\coshmax & $<$ 0.5 & $<$ 0.6 & $<$ 0.7 \\
\hline
\chitwotop & $<$ 4 & $<$ 6 & $<$ 8 \\
\hline
$\pt^{\ttbar}/\met$ & (0.7, 1.2) & \multicolumn{2}{c|}{(0.5, 1.2)} \\
\hline
\end{tabular}
\end{table}

The main contribution to the SM background in the signal regions originates from single-top-quark events in the $tW$ channel and  \ttbar\ events,
with a lepton missed by the reconstruction algorithms,
and $Z\to\nu\nu$ events ($Z$+jets, $\ttbar{+}Z$). Unlike in the \noLhigh\ SRs, \ttbar\ is the dominant contribution in the \noLlow\ selections, due to the lower \met\ requirements.
Dedicated control regions are used to aid the background estimation in the signal regions for all
dominant processes. These CRs are similar to the ones in the  \noLhigh\ analysis for single top quarks in the $tW$ channel, $Z$+jets and $\ttbar{+}Z$.
No control regions are defined for $W$+jets since it represents less than 5\% of the total background in all three \noLlow\ signal regions.
Due to the importance of the \ttbar\ background in the \noLlow\ selections, $\ttbar$ events with extra $b$-hadrons ($\ttbar{+}b$) are treated separately from \ttbar\ events without extra jets or events with extra light-flavour jets (referred to as \textit{other}) in the \noLlow\ analysis, and distinct control regions are defined for each of the two components. This is done because $\ttbar\text{+}b$ events are more likely to pass the signal region selections than \ttbar (other) events because of the presence of extra $b$-jets, which makes it more difficult to isolate and reconstruct the \ttbar\ system. In practice, the \ttbar\ (other) and $\ttbar{+}b$ event distributions exhibit a significant shape difference in the \coshmax\ variable. Moreover, the observed simulation mismodelling for $\ttbar{+}b$ events is larger than for \ttbar\ (other) events and can be corrected for by separating the two components.

The background estimates are validated in dedicated, non-overlapping,
validation regions, which require zero leptons and are orthogonal to the signal region selections.
In these regions, the background prediction agrees with the data to within one standard deviation ($1\sigma$).
More details are given in Appendix~\ref{app:bkg}.

The expected numbers of events are estimated in a simultaneous profile likelihood fit to all \noLlow\ CRs,
and are shown in Table~\ref{tab:lowmetyield}. The observed data are compatible with the prediction, agreeing to within $2\sigma$ in each signal region.

The results presented in this paper show the final combination of the \noLlow\ and \noLhigh\ analyses.
The details of this combination and the single-channel individual limits are discussed in Appendix~\ref{sec:lowmetresults}.

\begin{table}[htbp]
\centering
\scriptsize
\caption{Expected and observed numbers of events in SR0X, SRWX and SRTX. The background
yields and uncertainties are shown after the simultaneous profile likelihood  fit to all \noLlow\ CRs.
The category `\ttbar\ (other)' represents \ttbar\ events without extra jets or events with extra light-flavour jets.
`Other' includes contributions from $\ttbar{+}W$, $tZ$ and $tWZ$ processes.
The quoted background uncertainties include both the statistical and systematic contributions, while the signal uncertainties are purely statistical.}
\label{tab:lowmetyield}
\begin{tabular}{p{0.26\textwidth-2\tabcolsep}>{\centering}p{0.25\textwidth-2\tabcolsep}>{\centering}p{0.25\textwidth-2\tabcolsep}>{\centering\arraybackslash}p{0.25\textwidth-2\tabcolsep}}
\hline
\textbf{Process} & \textbf{SR0X} & \textbf{SRWX} & \textbf{SRTX} \\
\hline
Observed data & 60\leavevmode\phantom{.00 $\pm$ 0.00} & 74\leavevmode\phantom{.00 $\pm$ 0.00} & 36\leavevmode\phantom{.00 $\pm$ 0.00} \\
\hline
Expected SM events & 45\bblank\ $\pm$ 8\bblank & 59\bblank\ $\pm$ 6\bblank & 28\bblank\ $\pm$ 5\bblank \\
\hline
\ttbar\ (other) & 14\bblank\ $\pm$ 4\bblank & 15\bblank\  $\pm$ 4\bblank  & 9.4\blank\ $\pm$ 3.5  \\
$\ttbar{+}b$      & 10\bblank\ $\pm$ 7\bblank & 15.0\blank\ $\pm$ 3.1\blank & 7.2\blank\ $\pm$ 2.8  \\
Single-top      & 3.8\blank\ $\pm$ 3.0      & 4.3\blank\  $\pm$ 2.6       & 1.9\blank\ $\pm$ 1.5  \\
$Z$+jets        & 8.0\blank\ $\pm$ 1.6      & 12.1\blank\ $\pm$ 2.3\blank & 3.1\blank\ $\pm$ 0.8  \\
$W$+jets        & 1.6\blank\ $\pm$ 1.1      & 2.7\blank\  $\pm$ 2.1       & 0.6\blank\ $\pm$ 0.6  \\
$\ttbar{+}Z$      & 5.9\blank\ $\pm$ 1.0      & 7.8\blank\  $\pm$ 1.3       & 5.3\blank\ $\pm$ 1.1  \\
Diboson         & \blank0.28 $\pm$ 0.20     & 0.7\blank\  $\pm$ 0.4       & \blank0.30 $\pm$ 0.19 \\
Other           & \blank0.55 $\pm$ 0.15     & \blank0.88  $\pm$ 0.24      & \blank0.70 $\pm$ 0.22 \\
\hline
Pre-fit \ttbar      & 15\leavevmode\phantom{.00 $\pm$ 0.00} & 17\leavevmode\phantom{.00 $\pm$ 0.00} & 9.8\leavevmode\phantom{0 $\pm$ 0.0}   \\
Pre-fit $\ttbar{+}b$  & 7\leavevmode\phantom{.00 $\pm$ 0.0}   & 11.5\leavevmode\phantom{0 $\pm$ 0.00} & 5.6\leavevmode\phantom{0 $\pm$ 0.0}   \\
Pre-fit Single-top  & 7.1\leavevmode\phantom{0 $\pm$ 0.0}   & 8.2\leavevmode\phantom{0 $\pm$ 0.0}   & 3.6\leavevmode\phantom{0 $\pm$ 0.0}   \\
Pre-fit $Z$+jets    & 6.1\leavevmode\phantom{0 $\pm$ 0.0}   & 9.2\leavevmode\phantom{0 $\pm$ 0.0}   & 2.3\leavevmode\phantom{0 $\pm$ 0.0}   \\
Pre-fit $\ttbar{+}Z$  & 5.9\leavevmode\phantom{0 $\pm$ 0.0}   & 7.9\leavevmode\phantom{0 $\pm$ 0.0}   & 5.4\leavevmode\phantom{0 $\pm$ 0.0}   \\
\hline
Benchmark signal models \\
\hline
DM $m(\phi, \chi) = (10, 1)\,\GeV$ & 27.4\blank\ $\pm$ 2.4\blank & 33.2\blank\ $\pm$ 2.2\blank & 27.5\blank\ $\pm$ 2.2\blank \\
DM $m(a, \chi) = (50, 1)\,\GeV$    & 18.8\blank\ $\pm$ 1.3\blank & 22.6\blank\ $\pm$ 1.5\blank & 10.6\blank\ $\pm$ 1.0\blank \\
$H \rightarrow$ inv  ($\br=100$\%)              & 10.52       $\pm$ 0.34      & 17.1\blank\ $\pm$ 0.4\blank & 12.1\blank\ $\pm$ 0.4\blank \\
\hline
\end{tabular}
\end{table}

\paragraph{The \oneL\ analysis}

This analysis requires exactly one lepton ($e$ or $\mu$), at least
four jets, two of which must be \btagged, and $\met > 230\;\GeV$, and was
designed to target spin-0 DM models.
The \MET significance $\mathcal{S}$ must be above 15 and, only for this
analysis, it considers only jets and leptons in the
events and their resolution, as described in Ref.~\cite{SUSY-2013-15}.
A recursive variable-radius reclustering algorithm~\cite{SUSY-2016-16} is applied to the jets
to identify at least one large-variable-radius jet loosely consistent with a
top quark ($m_{\mathrm{top}^{\mathrm{reclustered}}} > 150\;\GeV$). The use of a variable-radius algorithm, instead of a
fixed-radius one, increases the acceptance of both highly boosted events and less boosted events when no explicit categorisation is
performed.
In addition, a requirement on the `topness' likelihood variable~\cite{Graesser:2012qy} is used to distinguish between the signal and dileptonic decays in SM \ttbar\ events where one of the leptons is misidentified or outside the acceptance. This variable quantifies how well each event satisfies the dileptonic \ttbar\ hypothesis, using the top quark and $W$ boson mass constraints and a requirement that the centre-of-mass energy of the event is minimised.
The \MET\ triggers were used to select
data that then populate this SR.  This region is
divided into four disjoint regions according to the azimuthal distance between
the \met and the lepton momentum,
$\Delta\phi(\vec{p}_{\textrm{T}}^{\textrm{miss}},\ell)$, which is presented in Figure~\ref{fig:srplots_1L}
and is found to be larger for pseudoscalar mediator models. The binning also maximises the sensitivity for
scalar mediator models, which are more similar to the
background but are characterised by a larger production cross section at low masses.
An additional requirement of
$\Delta\phi(\vec{p}_{\textrm{T}}^{\textrm{miss}},\ell) > 1.1$ is applied to suppress the SM background.
The dominant backgrounds, $t\bar{t}$ and $t\bar{t}Z$, are estimated
by means of dedicated CRs.

\paragraph{The \twoL\ analysis}

The last analysis considers events with two opposite-charge leptons
($e$ or $\mu$), at least one \btagged jet and large values of  \met  significance ($\mathcal{S} > 12$), exploiting events
collected with dilepton triggers.
Events are then separated into two categories depending on whether the two leptons have the same or
different flavour, and  in the same-flavour selection an additional
requirement of $|m_{\ell\ell}-m_Z|>20$~\GeV~is added to suppress
the \(Z\)+jets background. In this selection, the main discriminating
variable is the leptonic stransverse mass
$m_{\textrm{T2}}$~\cite{Barr:2003rg,Lester:1999tx}, which is used to
bound the individual masses of a pair of identical particles that are each
presumed to have decayed into one visible and one invisible
particle. This quantitiy is used to bound dileptonic top pair decays. To maximise the search sensitivity, the $m_{\textrm{T2}}$
spectrum is divided into six bins, starting from 110~\GeV. The
$m_{\textrm{T2}}$ distribution for selected events with two leptons
with the same flavour is presented in Figure~\ref{fig:srplots_2L}.
In this search, the main backgrounds are from $t\bar{t}$, $t\bar{t}Z$,
single-top-quark $tW$-channel, $Z$+jets, and diboson
processes. These backgrounds are estimated with MC simulations
and normalised with data in orthogonal CRs for
the dominant contributions ($t\bar{t}$ and $t\bar{t}Z$), while
the background arising from fake/non-prompt leptons is estimated
directly from the data.

\subsection{Orthogonalisation}

In order to combine the results of the different searches,
the searches are required to be statistically independent and
any possible overlaps of kinematic regions were investigated and removed as described in the following.
The three analysis channels are disjoint because of their requirements on
lepton multiplicity. The \noLhigh\ and \noLlow\ channels
are kept orthogonal by the requirements on the large-radius jet as
well as on the \met and its significance, $\mathcal{S}$. In addition, one of the $Z$+jets CRs
in the \noLhigh\ analysis, denoted by CRZAB-T0 in Ref.~\cite{SUSY-2018-12},
is not considered and a single control region, CRZAB-TTTW, is used to normalise
the $Z$+jets process in all SRs of the \noLhigh\ analysis. This has negligible
impact on the \noLhigh\ analysis results and it is done
to ensure orthogonality between the $Z$+jets CRs in the \noLhigh\ and \noLlow\
analyses, as those events are used to normalise the $Z$+jets background in the \noLlow\ analysis. To the same end,
the $Z$+jets CR in the \noLlow\ analysis only selects events with either
$N_{\textrm{large-radius jet}} < 2$ or subleading large-radius
jet mass~$< 60$~\GeV.

The CRs used to normalise the $t\bar{t}Z$ background overlap.
The three analysis channels share a common strategy to determine the amount of
$t\bar{t}Z$ (with \(\Zboson\rightarrow\nu\nu\)) background in their SRs.
The strategy is to construct CRs requiring three charged leptons in order to maximise their $t\bar{t}Z$ (with \(\Zboson\rightarrow\ell\ell\))
event content, which once determined can be scaled by the ratio of
\(\Zboson\rightarrow\nu\nu\) to \(\Zboson\rightarrow\ell\ell\) branching fractions.
These control regions differ only in minor selections adapted to the SR of each specific channel.
In the combination, the $t\bar{t}Z$ estimation is harmonised by using the
most inclusive CR$_{\ttbar Z}$, from the \twoL\ analysis~\cite{SUSY-2018-08}, as
a common CR across all channels. The fitted normalisation parameter value
obtained in the combination is consistent within 1\% with the one published in Ref.~\cite{SUSY-2018-07}.

\FloatBarrier
\section{Statistical combination and uncertainties}
\label{sec:comb}

The statistical combination of the analyses considered in this paper
consists of maximising a profile likelihood
ratio~\cite{Cowan:2010js} constructed from the product of the
individual analysis likelihoods:

\[
\Lambda(\alpha;\theta) = \frac{L\left(\alpha,\hat{\hat{\theta}}(\alpha)\right)}{L\left(\hat{\alpha},\hat{\theta}\right)}.
\]

The \(\alpha\) and \(\theta\) parameters represent, respectively, the parameter of interest
and the nuisance parameters. In the numerator, the nuisance parameters
are set to their profiled values \(\hat{\hat{\theta}}(\alpha)\), which
maximise the likelihood function for fixed values of the parameter of
interest \(\alpha\). In the denominator, both the parameter of
interest and the nuisance parameters are set to the values that
jointly maximise the likelihood: \(\hat{\alpha}\) and
\(\hat{\theta}\), respectively.

For the DM signal model interpretations, upper limits on the signal cross section are calculated following
the $\mathrm{CL_s}$ formalism, using the profile likelihood
ratio as a test statistic.
The parameter of interest is the overall signal strength, defined as a scale factor
multiplying the cross section predicted by the signal hypothesis, and it is bounded from below by zero.
The final result is provided as a ratio of the lowest excluded
signal cross section to the predicted cross section with all couplings set to unity.
For the \hinv\ signal model interpretation, the branching fraction \brinv is considered as the parameter of interest \(\alpha\), following the
implementation described in Refs.~\cite{Verkerke:2003ir,Moneta:2010pm}.

As described in Section~\ref{sec:analysis}, for each channel the estimation
of the dominant SM backgrounds is aided by means of dedicated control regions
that constrain free-floating normalisation factors for each
of these backgrounds.

Systematic
uncertainties are modelled in the likelihood function as nuisance
parameters $\theta$ constrained by Gaussian or log-normal probability
density functions~\cite{HIGG-2015-03}.

Three types of
sources of systematic uncertainty are considered: detector-related (experimental)
uncertainties, uncertainties related to the modelling of SM background
processes, and uncertainties related to the modelling of the signal
processes. Regarding the experimental and SM modelling uncertainties,
all details are given in Refs.~\cite{SUSY-2018-12,SUSY-2018-07,SUSY-2018-08} respectively
for the zero-, one- and two-lepton channels. The \noLlow\ channel includes the same
uncertainties as the \noLhigh\ channel and, in addition, uncertainties associated with the
$b$-jet trigger efficiencies. The typical size of these uncertainties is a few percent.
All analyses use common event-quality criteria and object
reconstruction and identification definitions. For this reason, all
experimental systematic uncertainties are treated as correlated across
channels in the statistical combination. The dominant sources of
experimental systematic uncertainty in the combination are the uncertainties
related to the jet energy scale and resolution, followed by either flavour-tagging
uncertainties or uncertainties related to the missing transverse momentum, depending on the analysis channel.

\begin{figure}
\centering
\includegraphics[width=.9\textwidth]{fig_04.pdf}
\caption{Summary of the total uncertainty in the background prediction for each SR
of the \noLlow, \noLhigh, \oneL, and \twoL\ analysis channels in the statistical combination. Their dominant contributions are indicated by individual lines. Individual uncertainties can be correlated, and do not necessarily add up in quadrature to the total background uncertainty.}
\label{g:SummUncertainties}
\end{figure}

Uncertainties in the modelling of the SM background processes in MC
simulation and their theoretical cross-section uncertainties are also
taken into account. All modelling uncertainties are treated as
uncorrelated across different channels as they probe different regions
of the available phase space.

Uncertainties related to the MC modelling of the DM signals include fragmentation and renormalisation scale uncertainties, and the uncertainties
related to the modelling of the parton shower. The impact
of these uncertainties varies from 10\% to 25\%. Uncertainties related to the $t\bar{t}H$ with \hinv\ signal modelling also include
fragmentation and renormalisation scale uncertainties, parton shower uncertainties and PDF uncertainties.
Among these, scale uncertainty effects, which are evaluated in the simplified template cross-section formalism~\cite{deFlorian:2016spz,HIGG-2020-23},
are the dominant contribution and range between 7\% and 17\%.
Signal modelling uncertainties are treated as fully correlated across analysis channels.

All sources of uncertainty in the SM backgrounds are summarised in Figure~\ref{g:SummUncertainties}.
In most of the SRs, the dominant systematic uncertainties are the ones related to theory predictions and
MC modelling, while jet uncertainties are the dominant experimental ones.
No significant difference from either the composition or the value of the total uncertainty presented in the published individual analyses is observed.

\section{Exclusion limits}
\label{sec:results}

Exclusion limits at 95\% CL are presented in
Figures~\ref{resscalar}~and~\ref{respseudo} for DM models with a spin-0 scalar or pseudoscalar mediator particle, respectively.
The three individual channels are also presented for comparison.
The \noL\ limits are the result of the statistical combination of
the \noLlow\ and \noLhigh\ SRs.
The \noLlow\ selection improves the expected scalar (pseudoscalar) mediator stand-alone cross-section limit of the \noLhigh\ by
up to 15\% (5\%) and it is strongest for mediator masses values around $10\;\GeV$.
Details of the comparison can be found in Appendix~\ref{sec:lowmetresults}.

The signal generation considered in these results includes both the top-quark-pair final states (DM+$\ttbar$) and single-top-quark final
states (DM+$tW$ and DM+$tj$). The limits are expressed in terms of the ratio
of the excluded cross section to the nominal cross section for
a coupling of $g = g_q = g_{\chi} = 1$. With these assumptions, scalar DM models are characterised by
a higher cross section than for pseudoscalar DM models with low mediator masses~\cite{Haisch:2016gry}, while the two
models have very similar cross sections beyond the top quark decay threshold ($m(\phi)$ or $m(a) \sim 2\cdot m_t$).
A DM particle mass of $1\;\GeV$ is considered, although the
results are valid as long as the mass of the mediator is larger than twice the mass of the DM particle. The solid
(dashed) lines show the observed (expected) exclusion limits
for each individual analysis and their statistical
combination. For scalar (pseudoscalar) DM models, the combination extends the
excluded mass range by $50~(25)\;\GeV$ beyond that of the best of the
individual analyses, excluding mediator masses up to $370\;\GeV$.
In addition, the combination improves the expected
cross-section limits by 14\% and 24\%, for low-mass scalar and pseudoscalar
DM mediators, respectively. This directly translates into more stringent
exclusion limits on the couplings.
When only the associated production of DM and two top quarks is considered
in the interpretation of the results, the excluded scalar (pseudoscalar) mediator mass range
obtained from the combination is reduced by
$70\,(20)\,\GeV$ relative to the sensitivity of the combination as reported in Figures~\ref{resscalar}~and~\ref{respseudo}.
As the production of DM in association with a single top quark is most relevant for
higher masses in the scalar mediator models~\cite{Pinna:2017tay},
the impact of this process for masses below $50\;\GeV$ is
negligible. In contrast, for the pseudoscalar mediator models, the ratio of
single-top-quark channel to \ttbar\ channel cross sections is relatively
constant~\cite{Pinna:2017tay}. When considering only DM+$\ttbar$ associated production,
the cross-section upper limit weakens by about 18\% over the
whole mass range.

\begin{figure}[ht!b]
\begin{center}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_05a.pdf}
\caption{}
\label{resscalar}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_05b.pdf}
\caption{}
\label{respseudo}
\end{subfigure}
\caption{Exclusion limits for colour-neutral (a) scalar or (b) pseudoscalar mediator
dark matter models as a function of the mediator mass
$m(\phi)$ or $m(a)$ for a DM mass $m_{\chi} =  1~\GeV$. Associated
production of DM with both single top quarks ($tW$ and $tj$
channels) and top quark pairs is considered. The limits are
calculated at 95\% CL and are expressed in terms of the ratio
of the excluded cross section to the cross section for
a coupling assumption of $g = g_q = g_{\chi} = 1$.  The solid
(dashed) lines show the observed (expected) exclusion limits
for each individual channel and their statistical combination.}
\end{center}
\end{figure}

The negative logarithmic profile likelihood ratios
$-2\,\Delta\!\ln(\Lambda)(\brinv;\theta)$ as a function of \brinv from
the individual analyses and their combination are shown in Figure~\ref{fig:lhscan_rtwo}.

Expected results are obtained using the Asimov dataset technique and calculated from asymptotic formulae~\cite{Cowan:2010js}.
The best-fit values of \brinv for the individual analyses are
compatible, agreeing to within one standard deviation. Their statistical combination
yields a best-fit value of $0.08^{+0.15}_{-0.15}$, consistent with the SM prediction of 0.12\%.
The combined observed 95\% CL upper limit on  \brinv  is $0.38$ while
the expected value is $0.30^{+0.13}_{-0.09}$. The individual analysis results are presented in Table~\ref{tab:DMtt_results},
while the details of the \noL\ combination are reported in Appendix~\ref{sec:lowmetresults}.
The overall uncertainty is dominated by the statistical uncertainty of the
data and, to a lesser extent, by systematic uncertainties associated with the modelling
of the SM processes and jet-related uncertainties.
Higgs boson invisible decays represent a specific case of the DM simplified models
considered in the previous section, where the mass of the scalar mediator is assumed
to be $125\;\GeV$. The two results are compatible with each other, when taking into account
the different order of accuracy used in event generation for the \hinv\ model.

\begin{figure}[ht!b]
\begin{center}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_06a.pdf}
\caption{}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{fig_06b.pdf}
\caption{}
\end{subfigure}
\end{center}
\caption{(a) The expected negative logarithmic profile likelihood
ratios $-2\,\Delta\!\ln(\Lambda)$ as a function of \brinv for
each of the three channels and their statistical combination
and (b) these likelihood ratios for the observed data.
\label{fig:lhscan_rtwo}
}
\end{figure}

\begin{table} [htbp!]
\caption{
\label{tab:DMtt_results}
Summary of results from direct searches for invisible decays of the
125~\GeV\ Higgs boson in the \tth topology using 139~\ifb\ of Run~2
data, and their statistical combination. Shown are the best-fit values
of \brinv, computed as not being bounded below by zero, for consistency with previous results~\cite{HIGG-2015-03}.
Observed and expected upper limits on \brinv at the 95\% CL are computed with the $\mathrm{CL_s}$ method and are new results
with respect to the individual analysis papers quoted in the last table column. The corresponding Asimov datasets for the expected
results are constructed using nuisance parameter values from a fit to
data with $\brinv = 0$, and the quoted uncertainty corresponds to the
68\% confidence interval.
}
\begin{center}
\begin{tabular}{l | c  c  c |l}
\hline
\multirow{2}{*}[1.5mm]{Analysis} & Best fit & Observed & Expected  & \multirow{2}{*}[1.5mm]{Reference} \vspace{-3mm} \\
& \brinv & upper limit & upper limit & \\
\hline
tt0L   		    &      ~~\,$0.48^{+0.27}_{-0.27}$      &     0.95   & $0.52^{+0.23}_{-0.16}$  & \cite{SUSY-2018-12}, this document \\
tt1L   		    &   $-0.04^{+0.35}_{-0.29}$       &     0.74   & $0.80^{+0.40}_{-0.26}$  & \cite{SUSY-2018-07}, this document \\
tt2L    	            &    $-0.08^{+0.20}_{-0.19}$    &     0.36   & $0.40^{+0.18}_{-0.12}$  & \cite{SUSY-2018-08}, this document \\
\hline
\(\ttbar H\) comb. &       ~~\,$0.08^{+0.15}_{-0.15}$     &    0.38   & $0.30^{+0.13}_{-0.09}$ & This document \\
\hline
\end{tabular}
\end{center}
\end{table}

\FloatBarrier

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

In summary, a statistical combination of three analyses using 139~fb$^{-1}$ of $pp$ collisions delivered by the LHC
at a centre-of-mass energy of 13\;\TeV\ and collected by the ATLAS detector is presented.
The three analyses are all designed to select events with two top quarks and invisible particles, and consider all
possible light lepton multiplicities arising from the decays of the two top quarks.

The statistical combination is used
to set 95\% confidence-level constraints on spin-0 simplified
dark matter models. All production modes with top quarks  in the final state (DM+$\ttbar$, DM+$t$) are considered.
For scalar (pseudoscalar) dark matter models, the combination extends the excluded mass range by
$50~(25)\;\GeV$ beyond that of the best of the individual channels, excluding mediator masses up to $370\;\GeV$ with all couplings set to unity. In addition, the combination improves the observed coupling exclusion limit by 24\%, assuming a pseudoscalar mediator mass of $10\;\GeV$.

The specific case where the mediator corresponds to the SM $125~\GeV$ Higgs boson
is also considered when interpreting the results presented in this paper.
An upper limit on the Higgs boson invisible branching ratio of 0.38 ($0.30^{+0.13}_{-0.09}$) is observed (expected) at 95\% confidence level.

\clearpage


\section{Introduction}
Autonomous agents are increasingly developed using machine learning (ML), most often with a reinforcement learning (RL) approach. In the RL setting, ML based agents explore their environment freely, with the aim of performing a task or solving a problem. Many control problems greatly benefit from the RL paradigm, since it alleviates many of their inherent challenges. It lets the agent itself be in charge of data collection, meaning that problems with large state spaces can be sufficiently explored without having to hand-tailor a data set to the problem. This can make the agent more robust, despite the fact that most real-life control environments cannot be fully modelled. RL methods also have the benefit of a goal-oriented approach to problems, meaning that we do not need to specify individual decisions in order to reach the end-goal.

While RL for fitting neural network models is an effective approach for creating autonomous agents, it introduces a problem related to human understanding: the resulting models are not interpretable to humans, opening up questions regarding \emph{what} the agent has learned, and whether it has internalised knowledge human domain experts know to be important. Based on the data processing inequality (see e.g.~\cite{dataprocessinginequality}), we know that an ML model cannot add to the information contained in its input data.  An important goal for explainable AI (XAI) in the context of autonomous agents is therefore to detect what knowledge these agents have conceptualised. We argue that investigating what a model has learned and whether this aligns with human domain knowledge is crucial for human oversight, and ensuring the stability and safety of of autonomous agents.

The game of chess is the most thoroughly explored environment in the history of artificial intelligence (AI).
We argue that chess is also highly suitable for developing and testing both training and explanation methods for control, for the following reasons.
Firstly, the dimensionality of the game is too large for complete exploration of the state space. The RL agent must therefore find an exploration exploitation trade-off in order to succeed in solving its task, as is the case in most real-life control settings.
Secondly, the game is extremely complex, meaning that the objective shifts from finding an optimal solution, to producing efficient heuristics for approximating optimal play.
Thirdly, the entire game with all its complexity can be completely simulated in a fast and efficient manner, and -- with the environment provided in this work -- with the desired dimensionality and thus computational cost. This means that it provides a problem that is complex enough for these methods to be relevant, while not being exceedingly computationally expensive, like many other control-problems.
Finally, the game features a large amount of available expert knowledge, meaning that explanations of the agent can be readily evaluated.

Today, AlphaZero \citep{alphazero_original} is the state-of-the-art ML based chess playing agent, demonstrating RL agents to be capable of outperforming human experts, and matching the strength of the best chess engines in existence. Recent work, available in the form of a preprint \citep{alphazero_concepts_2021}, demonstrates that it is possible to probe which \emph{human concepts} AlphaZero has learned through self-play. This methodology is likely to prove itself highly relevant in control contexts, providing us with the opportunity to investigate what domain knowledge black-box autonomous agents have acquired.

A considerable challenge in this regard is that the re-implementation of the kind of training loop required for such investigation is beyond the computational resources available to most research groups.
The classical game of chess is 8x8 dimensional, and requires the evaluation of approximately \({10^{123}}\) positions in order to determine the value of the initial position \citep{statespacechess}.
It is estimated that AlphaZero performed 44 million self-play games for its peak playing strength. With a large network, this is extremely resource intensive. Additionally, the concept detection method used in \cite{alphazero_concepts_2021} is based on outputs from the intermediary layer blocks of their network, meaning that the dimensionality of the data set depends of the size of the network. Given their specifications of 256x3x3 dimensions per layer, with a specified maximum sample size of the data set at \(5 \cdot {10^5}\), this means one needs to find a least-fit regressor on 2304 variables over such a data set for each such layer.

In order to mitigate this challenge and make this promising venue of investigating XAI methods for detecting concepts learned by RL agents in complex environments, we provide
\begin{itemize}
    \item[\ding{51}] a fast, open-source chess environment for easy simulation of any chess board size, starting position and piece types\footnote{Available piece types are limited to those included in standard 8x8 chess.},
    \item[\ding{51}] fully trained RL agents mastering game-play on smaller chess boards,
    \item[\ding{51}] code for probing concepts learned by trained agents, together with results showing the detected concepts.
\end{itemize}

The paper is organised as follows.
In Sec.~\ref{sec:methods}, we briefly outline the training of the RL agents, present the environment, and describe the concept detection method.
In Sec.~\ref{sec:results}, we present results for training agents in two different environments, and the evolution of the learned concepts throughout their training processes.
In Sec.~\ref{sec:analysis}, we analyse the learned concepts, and how they relate to the agents' learning processes.
In Sec.~\ref{sec:discussion}, we discuss limitations and possible developments of the concept detection methods, including potential relevance to other areas of application.
\section{Method}\label{sec:methods}
\subsection{Training the Agents}\label{sec:agents}
The goal is for an agent to master the game of chess via self-play using reinforcement learning. Our agent is a deep neural network with randomly initialised weights, meaning that it contains no knowledge at the beginning of the training loop. Its ability to play chess arises solely from experience during self-play. Through the exploration of different moves, it learns different playing styles, to recognise useful patterns on the board, and to evaluate board positions.

\begin{figure}
    \centering
    \includegraphics[width=0.35\textwidth]{figures/resnet.png}
    \caption{\label{fig:architecture} A block-diagram showing the architecture of the 6x6-ResNet model, as described in the text. The CNN-architecture for the 4x5-model has the same backbone, but without skip-connections and only three initial convolution layers. The 4x5 model has 32 filters per convolutional layer, and the 6x6-model has 64 filters per layer.}
\end{figure}

In this work, we present two agents; one trained on Silverman 4x5 chess, seen in Fig.~\ref{fig:board-4x5},
and one trained on Los Alamos 6x6 chess, seen in Fig.~\ref{fig:board-6x6}.
Both agents are feed forward convolutional neural networks (CNNs), first described by \cite{cnn}. The 6x6 model has added residual connections, akin to the architecture of AlphaZero~\cite{alphazero_original}. This has its origins in a standard ResNet-architecture, introduced by \cite{resnet}. The main architecture for both the models is shown in Fig. \ref{fig:architecture}.

Despite the reduced board size, compared to the usual 8x8 chess, both the 4x5 and the 6x6 chess variants are quite challenging to master. Our 4x5 agent is trained for 100 iterations, each iteration involving 250 games played. The 6x6 agent is trained for 400 iterations, each iteration involving 500 training games played. Both require tuning of the models' and the training loops' hyperparameters.\footnote{The full models and code used to produce them are available at \url{https://github.com/patrik-ha/explainable-minichess}}
By building the full search tree, we have found that optimal play results in a draw for the Silverman 4x5 variants.
Our 4x5-agent reaches this strategy of play, producing draws in almost all simulation games, which we -- together with the stagnating validation loss -- interpret as having identified and settled on a stable and ideal playing style.
In the following section, we provide further details about the training environment.
\subsection{Environment}
We provide a chess environment for reinforcement learning, with user-defined board-size and starting position, and a simple interface for extracting information about the internal state. The environment also allows customising castling-rules, and creates a representation of each position conforming to the input standard used by AlphaZero. This is highly relevant when used in conjunction with an RL training-loop, and also includes creating a one-to-one mapping for moves between the state-space of the environment, and the state space used for representation in the neural net.
The environment is designed for easy extendability. Although it currently only implements standard chess pieces, it can easily be extended with custom move sets and rules. This is particularly relevant for smaller chess variants, where additional artificial rules may be added in order to increase complexity without increasing the dimensionality.

The chess environment is written primarily in Python. We use a standard bitboard representation first proposed by \cite{bitboardorigin}, enhanced with ``magic-move bitboard generation'' introduced in \cite{kannan2007magic}, providing a significant speedup in move generation. These are both standard techniques for building compact and efficient chess simulators.
Python was chosen due to its ease of implementation, and convenience for the wider machine learning community. Although compiled languages, such as C++, are often faster, having access to the backbone of the environment allows direct access to computations done by the main game-play loop, which will be required later, see Sec. \ref{sec:concepts}. Almost all operations performed during game-play are bit-operations, owing mainly to the bitboard representation of the game. These low-level operations are done on NumPy-integers \citep{numpy}, since most bit-operations on these have bindings that are pre-compiled in C. This alleviates most of the performance loss from using Python. Some parts of the environment are compiled using Numba \citep{numba2015}, effectively providing a JIT-compiled version binding of the Python-code in C.

The environment is also coupled with a standard variant of Monte Carlo Tree Search (MCTS), first described in \cite{mcts}. The implementation utilises simple multi-threading, with minimal communication between the threads. In particular, a tree-paralellized version of MCTS is used, described in \cite{parallelmcts}. Additionally, our implementation exploits the fact that since our neural network is only updated after a given number of self-play iterations, subsequent simulation games between these updates are independent. TensorFlow Lite \citep{tensorflow2015-whitepaper} is utilised to provide fast neural network guidance without having to rely on batching predictions for the GPU. This is doubly important since our implementation requires a separate instance of the predictive model per thread.

This environment is used to train the two agents described in Sec.~\ref{sec:agents}.

\begin{figure}
    \centering
    \begin{subfigure}{0.20\textwidth}
        \includegraphics[width=\textwidth]{figures/4x5.png}
        \caption{\label{fig:board-4x5}}
    \end{subfigure}
    \hspace{3pt}
        \begin{subfigure}{0.25\textwidth}
        \includegraphics[width=\textwidth]{figures/6x6.png}
        \caption{\label{fig:board-6x6}}
    \end{subfigure}
    \caption{The two chess-variant implementations used in this work:
    (\protect\subref{fig:board-4x5}) Silverman 4x5 and
    (\protect\subref{fig:board-6x6}) Los Alamos 6x6.
    }
    \label{fig:boards}
\end{figure}
\subsection{Concept Detection}\label{sec:concepts}
We use the concept detection method used in \cite{alphazero_concepts_2021}, based on the idea introduced in \cite{tcav_2018}. Briefly described, concepts are detected by learning a set of logistic probes for a data set representing the concept of interest, using the activation outputs for each intermediary layer in the neural network. In the binary case, for a given concept \(C\), an intermediary layer of size \(m\) in the neural network, a given list of inputs \({I_0},{I_1}...{I_N}\), the corresponding activation outputs \({O_0},{O_1}...{O_N}\), and binary labels \({P_i} \in \left\{ {0,1} \right\}\) for each activation output, the aim is to find the best-fitting logistic regressor with weights \(\mathbf{w}\) and bias \(\mathbf{b}\) so that
\begin{equation}
    \left\| {\sigma \left( {\mathbf{w} \cdot {O_i} + \mathbf{b}} \right) - {P_i}} \right\|_2^2
\end{equation}
is minimized for all \({P_i}\) and \({O_i}\). Here, \(\sigma \) is the sigmoid function. To ensure that the concepts found are actually compactly represented, an $L1$-penalty weighted by \(\lambda \) is added to the weights, which in turn gives the final minimization objective:
\begin{equation}
    \left\| {\sigma \left( {\mathbf{w} \cdot {O_i} + \mathbf{b}} \right) - {P_i}} \right\|_2^2 + \lambda {\left\| \mathbf{w} \right\|_1} + \lambda \left| \mathbf{b} \right| \,.
\end{equation}

The presence of a concept in different layers of the neural network is calculated by splitting the data set representing the concept into a training set and a validation set, learning the logistic regressor on the training set and calculating its binary accuracy on the validation set.
In this case, for the logistic probe \(L\left(  \cdot  \right)\), the binary accuracy, corrected for random guessing, is
\begin{equation}
    \frac{2}{N} \left( {\sum_i^N {{H\left( {L\left( {{O_i}} \right) - 0.5} \right) - {P_i}} } } \right) - 1 \,,
\end{equation}
where $H$ is the Heaviside step function.

\begin{table}[hb]
    \centering
    \captionsetup{width=.9\linewidth}
    \caption{\label{table:concepts}Concepts probed in the chess agents.}
    \begin{tabular}{ll}
    \toprule
        Name & Description \\
        \midrule
        \texttt{has\_mate\_threat} & Checkmate is available \\
        \texttt{in\_check} & Is in check \\
        \texttt{material\_advantage} & Has more pieces than opponent \\
        \texttt{threat\_opp\_queen} & Opponent's queen can be captured\\
        \bottomrule
    \end{tabular}

\end{table}

We probe the two agents for four concepts in total, listed in Table \ref{table:concepts}. The concept probing data sets are created through large amounts of self-play between all of our model checkpoints. Then, a subset containing 10\% of the positions is randomly sampled from these games, and the positions labelled according to the concept they represent. This processed is repeated until enough positions are gathered to form a balanced data set\footnote{We perform a sanity check of our concept probing approach applying it to a data set with random concept labels, for which it predicts approximately $0$ in all layers.}.
This differs from the approach in \cite{alphazero_concepts_2021}, where concept data sets are obtained from a database of expert-level games.
We choose our approach for two reasons: Primarily, there is no database of recorded games for smaller chess variants. Secondly, we wish to be able to generate balanced data sets of arbitrary sizes, which is not possible when sampling from a fixed-size data set.
To avoid any bias in the concept data set in favour of the models being tested, we sample from a large variety of models, in addition to adding noise to the move selection process. This noise also ensures that each model is able to produce a large amount of potential games.

Each of our concept data sets consists of \(2.5 \cdot {10^5}\) positive and \(2.5 \cdot {10^5}\) negative samples. The validation ratio is \(0.2\), and the $L1$ weighting is \(\lambda  = 0.01\).

\section{Results}\label{sec:results}
The probed concepts for the 4x5 and 6x6 agents are shown in Figs.~\ref{fig:concepts4x5} and~\ref{fig:concepts6x6}. Here, we see concepts evolving over the course of training.
In Fig.~\ref{fig:concepts4x5-material}, we see that the 4x5-agent develops a representation for detecting whether it has a material advantage, meaning the ability to count the material value of the opponent's pieces compared to its own.
Soon after in the training process, proceeding approximately 30 iterations, the 4x5-agent learns to represent whether it is threatening the opponent's queen, see Figs.~\ref{fig:concepts4x5-queen}.
Later, after approximately 50 iterations, the agent learns to represent whether it has a potential mating attack, see Figs.~\ref{fig:concepts4x5-mate-threat}.
Surprisingly, we observe that the agent represents the state of being in check only weakly throughout the entire training process, see Fig.~\ref{fig:concepts4x5-in-check}, despite learning to play optimally.
Observe that all the probed concepts flatten out after about 100 training iterations, regardless of whether training continues. This indicates that this agent is highly unlikely to allocate more resources to represent these concepts provided more experience.

For the 6x6-agent, we see many of the same trends as for the 4x5 agent, compare Figs.~\ref{fig:concepts6x6-material}, \ref{fig:concepts6x6-mate-threat} and \ref{fig:concepts6x6-queen}.
A main difference is that \texttt{material\_advantage} as well as \texttt{has\_mate\_threat} are detectable already at the beginning of the training loop. This is likely due to the residual connections combined with the increased model size of the 6x6-agent, providing it with the capacity to represent the concepts, even in an untrained state. For models with skip connections, it is in general more likely that information from the data is maintained throughout the model, albeit unintentionally.
We also see that this agent learns to detect whether it is in check, see Fig.~\ref{fig:concepts6x6-in-check}, in contrast to the 4x5-agent. Both agents weakly represent this concept during the early stages of training, but only the 6x6-agent continues to develop this concept throughout the entire training process.

\begin{figure}
    \centering
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/4x5/material_advantage.png}
        \caption{}
        \label{fig:concepts4x5-material}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/4x5/has_mate_threat.png}
        \caption{}
        \label{fig:concepts4x5-mate-threat}
    \end{subfigure}
        \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/4x5/threat_opp_queen.png}
        \caption{}
        \label{fig:concepts4x5-queen}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/4x5/in_check.png}
        \caption{}
        \label{fig:concepts4x5-in-check}
    \end{subfigure}
    \caption{Concepts modelled by the 4x5 chess board agent, showing development of the model's ability to detect
    (\protect\subref{fig:concepts4x5-material}) whether the player to move has a significant material advantage,
    (\protect\subref{fig:concepts4x5-mate-threat}) whether the opponent could mate it in a single move,
    (\protect\subref{fig:concepts4x5-queen}) whether the opponent's queen is threatened,
    and
    (\protect\subref{fig:concepts4x5-in-check}) whether it is in check,
    developing in the course of the training checkpoints.
    }
    \label{fig:concepts4x5}
\end{figure}

\begin{figure}
    \centering
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/6x6/material_advantage.png}
        \caption{}
        \label{fig:concepts6x6-material}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/6x6/has_mate_threat.png}
        \caption{}
        \label{fig:concepts6x6-mate-threat}
    \end{subfigure}
        \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/6x6/threat_opp_queen.png}
        \caption{}
        \label{fig:concepts6x6-queen}
    \end{subfigure}
    \begin{subfigure}{0.24\textwidth}
        \includegraphics[width=\textwidth]{results/6x6/in_check.png}
        \caption{}
        \label{fig:concepts6x6-in-check}
    \end{subfigure}
    \caption{Concepts modelled by the 6x6 chess board agent, showing development of the model's ability to detect
    (\protect\subref{fig:concepts6x6-material}) whether the player to move has a significant material advantage,
    (\protect\subref{fig:concepts6x6-mate-threat}) whether the opponent could mate it in a single move,
    (\protect\subref{fig:concepts6x6-queen}) whether the opponent's queen is threatened,
    and
    (\protect\subref{fig:concepts6x6-in-check}) whether it is in check,
    developing in the course of the training checkpoints.
    }
    \label{fig:concepts6x6}
\end{figure}
\section{Analysis}\label{sec:analysis}
In an RL context, most of the early learning takes place because the agent accidentally succeeds, i.e.\ delivers checkmate. As a consequence, concepts connected to actions that are statistically more likely to win the game for a given agent, appear during the earliest stages of training.
For the 4x5-agent, we observe that it quickly creates and bolsters an internal representation of material imbalance, see Fig.~\ref{fig:concepts4x5-mate-threat}. The reason is most likely that this a simple yet highly predictive proxy for which player will win. For very weak agents, doing mostly random action selection strategies, it is highly beneficial to have more pieces than the opponent, as it implies controlling more squares than the opponent, which in turn gives a higher chance of winning. This is also reflected in the appearance and development of the ability to detect threats to the opponent's queen. On average, the queen controls the largest number of squares per piece, and therefore serves as a simple proxy for ``the ability to deliver checkmate'' when playing in a non-informed way. The same tendency was observed in \cite{alphazero_concepts_2021}.

We also observe the emergence of a canonical concept for the 4x5-agent, representing whether the player has a guaranteed win, shown in Fig.~\ref{fig:concepts4x5-mate-threat}. However, for the same agent, we do not see a strong emergence of the ability to detect whether the player to move is in check, despite the model seemingly understanding the concepts of threats and attacks (see Figs.~\ref{fig:concepts4x5-mate-threat} and~\ref{fig:concepts4x5-queen}). We expected that being in check would be an important concept to represent, as being in check significantly limits what moves are available for play. For the model to be able to predict legal moves in this situation, it would therefore be reasonable to assume that it needs to detect if it is in check. The reason for the weak representation of this concept is most likely that the ratio of legal to illegal moves is significantly smaller on this board than on 6x6, or even 8x8-boards. This would also explain why this concept is more strongly represented in our 6x6-agent model (see Fig.~\ref{fig:concepts6x6-in-check}), and in \cite{alphazero_concepts_2021}.

In the case of the 6x6-agent, we observe a localised emergence of concepts. In Figs.~\ref{fig:concepts6x6-queen} and~\ref{fig:concepts6x6-in-check}, either concept emerges primarily in the first and last layers. This is due to the architecture of the 6x6-agent's model, and agrees with the interpretation presented in \cite{linearclassifierprobes} as well as the observations in \cite{alphazero_concepts_2021}. For abstract concepts, like \texttt{threat\_opp\_queen}, it makes sense that these are represented early in the model. On the other hand, it makes sense that concepts that have immediate consequences for which moves are available, like \texttt{in\_check}, are represented close to the relevant output-heads.

We also observe that there is a difference between the two agents regarding \emph{when} during training each concept emerges. We also see that the concepts represented in the smaller agent reach a plateau. This is in contrast to the larger agent, which continues developing concepts even after finding a winning playing style. This is as expected, since a larger dimensionality implies that more examples and thus training time are needed to experience a sufficient number of situations in which the different concepts are relevant.

\FloatBarrier
\section{Discussion}\label{sec:discussion}
We have demonstrated that the concept probing method is capable of detecting whether deep neural network models contain a compact, linear representation of a given concept.
This allows probing neural network based agents for domain knowledge without having to include it in the training loop.
This method is therefore highly relevant for investigating, evaluating and explaining neural network based autonomous agents, as it reveals whether these have internalised domain knowledge as expected. In the following, we discuss challenges and limitations of the method.

While this was not a problem in our case, generating representative concept data sets might prove difficult. Ideally, the concept must be defined programmatically, so that concept examples can be distilled from some large sample-repository. This might be challenging, especially for large and complex problems. Also, it requires an efficient simulation-environment. In cases where the environment is a simulation of a real-world application, the simulation must also be able to generate concept samples that are representative of what the agent would encounter in its ``true environment''. Otherwise, concept detection might be hindered by subtleties of synthetic data.

Exploration is a fundamental aspect of the RL paradigm, opening up the possibility that the model has observed states that, from a domain expert's viewpoint, are unlikely to occur. This is important because it makes it more likely that the model is allowed to generalise its concepts, as opposed to if it was trained, e.g., in a supervised fashion on a curated data set.

In our implementation, we limit ourselves to probing for \emph{linearly} represented concepts, although there are no guarantees that the network represents its knowledge in a linearly separable manner. However, as observed in \cite{linearclassifierprobes}, it is expected that most neural networks with sufficient depth have will end up representing emergent concepts as simply as possible, i.e.\ linearly.

Finally, we wish to highlight is the difference between being able to detect concepts, and knowing how the represented concepts are utilised within the given model. In our case, we see that both models develop some representation of \texttt{material\_advantage}. For a human, this particular concept is usually used as a heuristic for the result of the game, meaning that the player holding fewer pieces is more likely to lose. However, we cannot guarantee that our model has made the same connection, although this aligns most readily with the relevant domain knowledge.
\section{Conclusion}
In this work, we have demonstrated that concept probing can be used to gain insight into what concepts a neural network model learns. We have also demonstrated the feasibility of using these methods to reason whether the model's internal representations align with existing domain knowledge. This is highly relevant for all control contexts in which an autonomous agent has learned to navigate in an environment about which humans possess domain knowledge. This implies state-of-the-art RL methods can be used to solve control tasks, without losing the ability to investigate how the learned agent understands its environment. This is particularly useful for problems where modelling the environment is challenging. Given that the agent performs its task well, this can also be used to provide new understanding of the environment in question.

In the context of neural network based agents, the concept detection method is both architecture and environment agnostic. As long as it is possible to define a concept function, the method can be applied to any layered model structure. This implies easy applicability in control contexts, without having to apply constraints on the environment or training loop. While we demonstrate the applicability of the concept detection methods to agents trained in our environment, an important next step will be applying the method to classical control problems, comparing the extracted concepts to domain knowledge.

\begin{ack}
We thank Vilde B.\ Gj{\ae}rum for valuable discussions and proof reading.
\end{ack}