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

This paper presents a search for a new vector boson $Z'$ in the four-muon (4$\mu$) final state with data recorded in proton-proton ($pp$) collisions at $\sqrt s$ = 13~\TeV{} by the ATLAS detector~\cite{PERF-2007-01} at the Large Hadron Collider (LHC), corresponding to an integrated luminosity of 139~fb$^{-1}$.
The new gauge boson $Z'$, which only interacts with the second and third generation leptons, is predicted by $L_{\mu}-L_{\tau}$ models~\cite{PhysRevD.44.2118} which extend the Standard Model (SM) with an additional $U(1)_{L_{\mu}-L_{\tau}}$ symmetry.
In recent years searching for such a new gauge boson $Z'$  has attracted attention in both theoretical and experimental communities~\cite{Elahi_2016,Elahi_2017,PhysRevLett.114.151801}
because it could address the observed muon magnetic dipole moment ($g$-$2$) anomaly~\cite{Fermilab:Muong2,Keshavarzi_2019,Venanzoni:2016hP,PhysRevD.73.072003}
through the \Zp loop corrections without contradicting other existing data.
The new $Z'$ is also ideally suited to address the lepton flavour anomalies
measured with the ratio of B-meson decays to muons and electrons,
$R_K = \mathcal{B}(B^+\rightarrow K^+\mu^+\mu^-) / \mathcal{B}(B^+\rightarrow K^+ e^+e^-)$~\cite{Aaij2022,Aaij_2013,Aaij_2014,PhysRevLett.125.011802, Blake_2017}.
Such anomalies suggest possible lepton flavor universality violation.
The $L_\mu - L_\tau$ models could provide an explanation for the observed results, since by construction, the $Z'$ arising from the models only couples to muons and taus, and does not couple to electrons.
These models also aim to probe outstanding questions related to dark matter and neutrino mass~\cite{Altmannshofer:2016jzy,Ma_2002,Langacker_2009}.

The \Zp kinematics, mass, and interactions (with the same coupling strength to the second and third lepton families), are described by the Lagrangian below:
\begin{equation*}
L_{Z'} = - {1\over 4} F_{\alpha\beta}F^{\alpha\beta} + {1\over 2} M^2_{Z'} Z'^\alpha Z'_\alpha
- g_{Z'} Z'_\alpha(\bar\ell_2\gamma^\alpha\ell_2 + \bar\mu\gamma^\alpha\mu
-\bar\ell_3\gamma^\alpha\ell_3 - \bar\tau\gamma^\alpha\tau ),
\end{equation*}
where $F_{\alpha\beta} = \partial_\alpha Z'_\beta - \partial_\beta Z'_\alpha$ is the $Z'$ field strength tensor;  $\ell_i = (\nu_i, e_i)^T$ ($i=2,3$, denoting the second and the third generation left-handed lepton doublets); $\mu$ and $\tau$ represent the right-handed muon and tau singlets; and $g_{Z'}$ (which will be referred to as \coupzp in the rest of this paper) is the interaction coupling constant.
The parameter space of ($m_{Z'}, \coupzp$) has not been strongly constrained in experiments since the $Z'$ does not directly couple to the electron nor to any quarks, hence it could not be directly produced from $e^+e^-$ and $pp$ colliding beams.

In $pp$ collisions at the LHC, the $Z'$ could be produced from final state radiation of $\mu$ or $\tau$ pairs of the Drell-Yan (DY) process as shown in \Fig{\ref{fig:feynman_new:a}} with a $4\mu$ final state, which provides the cleanest signature to search for the $Z'$.
For relatively low $Z'$ mass, the most promising experimental signature would be an excess of 4$\mu$ events with the invariant mass of one $\mu^+\mu^-$ pair peaking around the $Z'$ mass.
The major background comes from the SM 4$\mu$ production processes shown in \Fig{\labelcref{fig:feynman_new:b,fig:feynman_new:c,fig:feynman_new:d}}.
The \Zp could also be produced in \Wboson production through the DY process, $pp\rightarrow W \rightarrow Z' \mu \nu \rightarrow 3\mu + \nu$.
The experimental signature would have a final state of 3$\mu$ plus large missing transverse energy.
This final state is not included in this analysis.

\begin{figure}[!ht]
\vspace{-0.6cm}
\centering
\subfloat[]{
\includegraphics[width=0.24\textwidth]{fig_01a.pdf}
\label{fig:feynman_new:a}
}
\subfloat[]{
\includegraphics[width=0.24\textwidth]{fig_01b.pdf}
\label{fig:feynman_new:b}
}
\subfloat[]{
\includegraphics[width=0.22\textwidth]{fig_01c.pdf}
\label{fig:feynman_new:c}
}
\subfloat[]{
\includegraphics[width=0.24\textwidth]{fig_01d.pdf}
\label{fig:feynman_new:d}
}
\caption{Feynman diagrams of $Z'$ production through radiation in a Drell-Yan process (a), and of the corresponding SM background processes (b - d) with a 4$\mu$ final state.}
\label{fig:feynman_new}
\end{figure}

Both the ATLAS and CMS Collaborations have measured the cross-sections of the SM $Z \rightarrow 4\mu$ process~\cite{CMS:2012bw,Aad:2014wra,Sirunyan:2017zjc}.
The measurement by ATLAS was used by theorists to set limits in the parameter space of the $L_\mu - L_\tau$ model~\cite{Altmannshofer:2016jzy}.
The CMS Collaboration has directly searched for the $Z'$ boson in the mass region between 5 to 70~\GeV{} with the 4$\mu$ final state using 77.3~\ifb~of data~\cite{CMS-EXO-18-008}, and set upper limits on the $Z'$ to muon coupling strength, \coupzp, of 0.004 -- 0.3 at 95\% confidence level, depending on the $Z'$ mass.

The organization of this paper is as follows.
The ATLAS detector is described in the next section.
The dataset and Monte Carlo samples used in this analysis are detailed in \Sect{\ref{sec:data-MC}}.
Event reconstruction and selection, followed by background estimation from data, are described in \Sect{\ref{sec:simu}} and \Sect{\ref{sec:background}}.
Event classification using a deep learning approach is presented in \Sect{\ref{sec:dnn}}.
Systematic uncertainties and the statistical approach to interpret data to obtain the
results
are reported in \Sect{\ref{sec:syst}} and \Sect{\ref{sec:result}}.
The conclusion is given in the final section.

\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)\).
Angular distance is measured in units of \(\Delta R \equiv \sqrt{(\Delta\eta)^{2} + (\Delta\phi)^{2}}\).}

The ATLAS detector~\cite{PERF-2007-01} at the LHC covers nearly the entire solid angle around the collision point.\footnote{\AtlasCoordFootnote}
It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadron calorimeters,
and a muon spectrometer incorporating three large superconducting air-core toroidal magnets.

The inner-detector system (ID) is immersed in a \SI{2}{\tesla} axial magnetic field
and provides charged-particle tracking in the range \(|\eta| < 2.5\).
The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track,
the first hit normally being in the insertable B-layer (IBL) installed before Run~2~\cite{ATLAS-TDR-19,PIX-2018-001}.
It is followed by the silicon microstrip tracker (SCT), which usually provides eight measurements per track.
These silicon detectors are complemented by the transition radiation tracker (TRT),
which enables radially extended track reconstruction up to \(|\eta| = 2.0\).
The TRT also provides electron identification information
based on the fraction of hits (typically 30 in total) above a higher energy-deposit threshold corresponding to transition radiation.

The calorimeter system covers the pseudorapidity range \(|\eta| < 4.9\).
Within the region \(|\eta|< 3.2\), electromagnetic calorimetry is provided by barrel and
endcap high-granularity lead/liquid-argon (LAr) calorimeters,
with an additional thin LAr presampler covering \(|\eta| < 1.8\)
to correct for energy loss in material upstream of the calorimeters.
Hadron calorimetry is provided by the steel/scintillator-tile calorimeter,
segmented into three barrel structures within \(|\eta| < 1.7\), and two copper/LAr hadron endcap calorimeters.
The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules
optimised for electromagnetic and hadronic energy measurements respectively.

The muon spectrometer (MS) comprises separate trigger and
high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroidal magnets.
The field integral of the toroids ranges between \num{2.0} and \SI{6.0}{\tesla\metre}
across most of the detector.
Three layers of precision chambers, each consisting of layers of monitored drift tubes, covers the region \(|\eta| < 2.7\),
complemented by cathode-strip chambers in the forward region, where the background is highest.
The muon trigger system covers the range \(|\eta| < 2.4\) with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.

Interesting events are selected by the first-level trigger system implemented in custom hardware,
followed by selections made by algorithms implemented in software in the high-level trigger~\cite{TRIG-2016-01}.
The first-level trigger accepts events from the \SI{40}{\MHz} bunch crossings at a rate below \SI{100}{\kHz},
which the high-level trigger further reduces in order to record events to disk at about \SI{1}{\kHz}.

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{Dataset and Monte Carlo Simulations}
\label{sec:data-MC}

The data used for this analysis were recorded using single-muon and multi-muon triggers, corresponding to an integrated luminosity of 139~\ifb~after the application of data quality requirements~\cite{Aad_2020}.
The transverse momentum (\pT) thresholds for the single-muon trigger vary from 20 to 26 GeV, for the di-muon trigger from 10 to 14 GeV, and for the tri-muon trigger from 4 to 6 GeV, depending on the data-taking periods~\cite{TRIG-2018-01}.
The overall trigger efficiency for the $4\mu$ events selected in this analysis is higher than 98\%{}.

\subsection{Simulation of \texorpdfstring{$Z'$}{Z'} Production}
\label{sec:simulation_sig}

The $L_\mu - L_\tau$ model is used for Monte Carlo (MC) \Zp signal sample production, where the \Zp couples to the left-handed (LH) muon or tau leptons and their corresponding neutrinos, and to the right-handed (RH) muon or tau leptons.
In the model, the \Zp couplings to the first lepton families (electron and its neutrino) and all quarks are set to zero.
The couplings to second and third generation leptons are assumed to be identical, therefore the branching fractions of the \Zp decay to a pair of muons and a pair of muon neutrinos are set to $\frac{1}{3}$ and $\frac{1}{6}$, respectively.
The signal from tau decays in the 4$\mu$ final state is found to be negligible in this analysis and is not included as signal.

The signal events are generated with \MGNLO[2.7.3]~\cite{Alwall:2014hca} at leading-order (LO) accuracy in QCD  by using the \textit{Universal FeynRules Output} (UFO) format~\cite{Degrande:2011ua,Alloul:2013bka}.
The interactions~\cite{delAguila:1972795} mediated by a resonance $Z'$ which couples to the second and third generation leptons are used for four-muon signal event generation.
Based on theoretical calculations~\cite{Grazzini_2015, h2020resummed} for related processes,
the appropriate NNLO/LO correction factor K of 1.3 is used to correct the MC LO signal cross-sections.
In the signal production simulations, the \NNPDF[2.3nlo] set~\cite{Ball:2012cx} is used as parton distribution function (PDF) for the $pp$ collisions.
\PHOTOSpp[3.61]~\cite{Davidson:2010ew} is used as the photon emissions from electroweak vertices and charged leptons via a QED process.

The MC simulated events are generated for a range of masses and coupling parameters of the $L_\mu-L_\tau$ model.
The $Z'$ mass ranges from 5 to 81 GeV, and the value of coupling constant \coupzp ranges from 0.008 to 0.316, as summarized in Table~\ref{tab:signalSamples}.
The value of the \coupzp at each \Zp mass was chosen to allow a sensitive search
for a very small \Zp signal in the full Run 2 dataset.

With the chosen \coupzp, the natural width $\Gamma$ of the $Z'$
(smaller than the experimental resolution)
and the cross-section, $\sigma$
($pp\rightarrow Z' \mu^+\mu^- \rightarrow \mu^+\mu^-\mu^+\mu^-$), are calculated as listed in \Tab{\ref{tab:signalSamples}}.

\begin{table}[h!]
\caption{Summary of the chosen $Z'$ hypotheses and corresponding coupling, width, and cross-section (calculated at LO accuracy in QCD) at each mass point.}
\centering
\resizebox{0.9\textwidth}{!}{
\begin{tabular}{cccc|cccc}
\hline
$m_{Z'}$ {[}GeV{]} & \coupzp & $\Gamma$ {[}GeV{]}        & $\sigma$ {[}fb{]} & $m_{Z'}$ {[}GeV{]} & \coupzp & $\Gamma$ {[}GeV{]} & $\sigma$ {[}fb{]} \\ \hline
5                  & 0.0080  & \SI{2.45E-5}{}    \SI{}{} & 9.96              & 42                 & 0.0900  & \SI{2.71E-2}{}     & 13.38             \\
7                  & 0.0085  & \SI{3.99E-5}{}            & 7.06              & 45                 & 0.1000  & \SI{3.58E-2}{}     & 11.72             \\
9                  & 0.0090  & \SI{5.78E-5}{}            & 5.60              & 48                 & 0.1100  & \SI{4.62E-2}{}     & 9.96              \\
11                 & 0.0095  & \SI{7.89E-5}{}            & 4.65              & 51                 & 0.1200  & \SI{5.84E-2}{}     & 8.24              \\
13                 & 0.0100  & \SI{1.03E-4}{}            & 3.95              & 54                 & 0.1600  & \SI{1.10E-1}{}     & 10.07             \\
15                 & 0.0120  & \SI{1.72E-4}{}            & 4.45              & 57                 & 0.2000  & \SI{1.81E-1}{}     & 10.73             \\
17                 & 0.0140  & \SI{2.65E-4}{}            & 4.80              & 60                 & 0.2665  & \SI{3.39E-1}{}     & 12.92             \\
19                 & 0.0160  & \SI{3.87E-4}{}            & 5.00              & 63                 & 0.2680  & \SI{3.60E-1}{}     & 8.84              \\
23                 & 0.0240  & \SI{1.05E-4}{}            & 7.30              & 66                 & 0.2780  & \SI{4.06E-1}{}     & 6.50              \\
27                 & 0.0320  & \SI{2.20E-3}{}            & 8.50              & 69                 & 0.2890  & \SI{4.59E-1}{}     & 4.89              \\
31                 & 0.0400  & \SI{3.95E-3}{}            & 8.72              & 72                 & 0.3000  & \SI{5.15E-1}{}     & 3.80              \\
35                 & 0.0600  & \SI{1.00E-2}{}            & 12.82             & 75                 & 0.3000  & \SI{5.37E-1}{}     & 2.88              \\
39                 & 0.0800  & \SI{1.99E-2}{}            & 14.77             & 78                 & 0.3080  & \SI{5.89E-1}{}     & 2.40              \\
&         &                           &                   & 81                 & 0.3160  & \SI{6.44E-1}{}     & 2.08              \\
\hline
\end{tabular}
}
\label{tab:signalSamples}
\end{table}

The $Z'$ signal samples were simulated with the ATLAS fast simulation framework (Atlfast-II)~\cite{SOFT-2010-01}
to produce predictions that can be directly compared with the data.

\subsection{Simulation of Background Events}
\label{sec:simulation_bkg}

Dominant SM backgrounds in this analysis come from the SM $Z \rightarrow 4\mu$ processes where the four leptons have an invariant mass close to that of the $Z$ boson.
In the higher mass region the $ZZ^*$ production
contributes a sizable number of prompt $4\mu$ events.
In addition, there are very small contributions
from the Higgs boson, $t\bar{t}V$ ($V=W,Z$), and tri-boson ($VVV$) production processes.
These events are estimated with MC simulations.
The non-prompt muon background events, mostly coming from \Zjet, \ttbar and single-top-quark production processes, are estimated from data in this analysis, as described in \Sect{\ref{sec:background}}.
The \Zjet and \ttbar MC samples are also produced for background studies.

Samples of diboson final states ($q\bar q\rightarrow VV^{(*)}$), including the processes shown in Figures \ref{fig:feynman_new:b} and \ref{fig:feynman_new:c}, were simulated with the \SHERPA[2.2.2]~\cite{Bothmann:2019yzt} generator, including off-shell effects and Higgs boson contributions, where appropriate.
Fully leptonic final states and semileptonic final states, where one boson decays leptonically and the other hadronically, were generated using matrix elements at next-to-leading-order (NLO) accuracy in QCD for up to one additional parton and at LO accuracy for up to three additional parton emissions.
Samples for the loop-induced processes $gg \to VV^{(*)}$, shown in Figure~\ref{fig:feynman_new:d}, were generated using LO-accurate matrix elements for up to one additional parton emission for both the cases of fully leptonic and semileptonic final states.
The matrix element calculations were matched and merged with the \SHERPA parton shower based on Catani--Seymour dipole factorisation~\cite{Gleisberg:2008fv,Schumann:2007mg} using the MEPS@NLO prescription~\cite{Hoeche:2011fd,Hoeche:2012yf,Catani:2001cc,Hoeche:2009rj}.
The virtual QCD corrections were provided by the \OPENLOOPS library~\cite{Buccioni:2019sur,Cascioli:2011va,Denner:2016kdg}.
The \NNPDF[3.0nnlo] set of PDFs was used~\cite{Ball:2014uwa}, along with the
dedicated set of parton-shower parameters (tune) developed by the
\SHERPA authors.

The production of \ttbar events was modelled using the \POWHEGBOX[v2]~\cite{Frixione:2007nw,Nason:2004rx,Frixione:2007vw,Alioli:2010xd} generator at NLO in QCD with the \NNPDF[3.0nlo] PDF set and the \hdamp parameter\footnote{The \hdamp parameter is a resummation damping factor and one of the parameters that controls the matching of \POWHEG matrix elements to the parton shower and thus effectively regulates the high-\pt radiation against which the \ttbar system recoils.} set to 1.5 \mtop~\cite{ATL-PHYS-PUB-2016-020}.
The events were interfaced to \PYTHIA[8.230]~\cite{Sjostrand:2014zea} to model the parton shower, hadronisation, and underlying event, with parameters set according to the A14 tune~\cite{ATL-PHYS-PUB-2014-021} and using the \NNPDF[2.3lo] set of PDFs.
The associated production of top quarks with $W$ bosons ($tW$) was modelled by the \POWHEGBOX[v2]~\cite{Re:2010bp,Nason:2004rx,Frixione:2007vw,Alioli:2010xd} generator at NLO in QCD using the five-flavour scheme and the \NNPDF[3.0nlo] set of PDFs.
The diagram removal scheme~\cite{Frixione:2008yi} was used to remove interference and overlap with \ttbar production.
The events were interfaced to \PYTHIA[8.230] using the A14 tune and the \NNPDF[2.3lo] set of PDFs.

The production of $V+$jets was simulated with the \SHERPA[2.2.1] generator using NLO matrix elements for up to two partons, and LO matrix elements for up to four partons, calculated with the Comix and \OPENLOOPS libraries.
They were matched with the \SHERPA parton shower using the \MEPSatNLO prescription with the set of tuned parameters developed by the \SHERPA authors.
The \NNPDF[3.0nnlo] set of PDFs was used and the samples were normalised to a next-to-next-to-leading-order (NNLO) prediction~\cite{Anastasiou:2003ds}.

The production of \ttV\ events was modelled using the \MGNLO[2.3.3] generator at NLO with the \NNPDF[3.0nlo] PDF.
The events were interfaced to \PYTHIA[8.210]~ using the A14 tune and the \NNPDF[2.3lo] PDF set.
The production of tri-boson ($VVV$) events was simulated with the \SHERPA[2.2.2] generator.
Matrix elements accurate to LO in QCD for up to one additional parton emission were matched and merged with the \SHERPA parton shower based on Catani--Seymour dipole factorisation using the \MEPSatNLO prescription.
Samples were generated using the \NNPDF[3.0nnlo] PDF set, along with the dedicated set of tuned parton-shower parameters developed by the \SHERPA authors.

The generated background MC samples were produced through the full ATLAS detector simulation
based on \textsc{Geant4}~\cite{Agostinelli:2002hh}.
The effect of multiple interactions in the same and neighbouring bunch crossings (\pileup) was modelled by overlaying the simulated hard-scattering event with inelastic $pp$ events generated with \PYTHIA[8.186]~\cite{Sjostrand:2007gs} using the \NNPDF[2.3lo] set of
PDF and the A3 set of tuned parameters~\cite{ATL-PHYS-PUB-2016-017}.
Simulated events were reweighted to match the pile-up conditions in the data.
All simulated events were processed using the same reconstruction algorithms and triggering requirements as used in data.

\section{Event Reconstruction and Pre-selection}
\label{sec:simu}

Proton-proton interaction vertices are reconstructed in events with at least two tracks, each with $\pT > 0.5~\GeV$.
The primary hard-scatter vertex is defined as the one with the largest value of the sum of squared track transverse momenta.

Muons are identified by matching tracks reconstructed in the MS to tracks reconstructed in the ID (referred to as combined muons).
To increase the muon reconstruction efficiency non-combined muon identification algorithms are also used in the analysis,
including using the MS stand-alone tracks in the region $2.5 <|\eta|<2.7$,
and matching the ID tracks with calorimeter hit information within $|\eta| < $0.1,
as well as using the ID tracks associated with at least one local track segment in the MS.
In the $4\mu$ event selection at most one of the selected muons can be a non-combined muon.
Each muon is then required to satisfy the ‘loose’ identification criteria \cite{MUON-2018-03-updated}.
Muons are required to be isolated using a particle-flow algorithm \cite{2017_pflow}
and associated with the primary hard-scatter vertex by satisfying
$|{d_0 \over \sigma_{d_0} }| < $ 3
and $|z_0\times \sin\theta| <$ 0.5~mm,
where $d_0$ is the transverse impact parameter
calculated with respect to the measured beam-line position, $\sigma_{d_0}$ its uncertainty, and $z_0$ is the longitudinal distance between the point at which $d_0$ is measured and the primary vertex.
The minimum muon \pT threshold is 3 \GeV.

In addition to muons, electrons, jets and missing transverse momentum (\MET) are also used to select control samples for background estimation in this analysis.
The reconstructions of these objects are described below.

Electrons are identified with a likelihood discriminator built from the  shower shapes of electron energy loss in the calorimeter, track-cluster matching, and some of the track quality distributions.
Each electron is required to satisfy the ‘medium’ likelihood identification criteria~\cite{Aad_2019}, as well as similar vertex and isolation requirements as muons.
The reconstructed electrons are required to have $p_{\textup{T}} >$ 7~GeV and $|\eta| < $ 2.5, excluding the transition region between the barrel and endcap calorimeters, $1.37 < |\eta| < 1.52$.

Jets are reconstructed with the anti-$k_t$ algorithm \cite{ATL-PHYS-PUB-2018-012, Fastjet} with a radius parameter of $R = 0.4$.
The jet clustering input objects are based on particle-flow \cite{PERF-2015-09} in the ID and the calorimeter.
Jets are required to have $p_{\textup{T}} > 20 ~\GeV$ and $|\eta|<2.5$.
Jets containing B hadrons, referred to as \bjets, are identified with a multivariate discriminant~\cite{ATL-PHYS-PUB-2016-012}.
To reduce the effect of pile-up an additional quality requirement based on the `Jet-Vertex-Tagger' algorithm (JVT) \cite{ATLAS-CONF-2014-018} is applied in jet identification.

\MET is determined as the magnitude of the negative vectorial sum of the transverse momenta of the selected and calibrated physics objects (including muons, electrons, photons, and jets including hadronically decaying tau-leptons) and the ID tracks must come from the main vertex and not associated with any physics object (soft term)~\cite{PERF-2016-07}.

Events containing at least four muons with kinematics consistent with $Z(Z^*/\gamma^*)\rightarrow 4\mu$ production are then selected as follows.
The four leading \pT-ordered muons are required to pass the \pT thresholds of 20, 15, 8, and 3~GeV, respectively.
If a muon is selected as a non-combined muon, its \pT must be greater than 15 GeV.
Any di-muon pair in the event must have an invariant mass $m_{\mu\mu}$ greater than 4 GeV and an angular separation $\Delta R$ larger than 0.2.
To search for the $Z'\rightarrow \mu^+\mu^-$ signature, two muon pairs are selected based on their invariant mass values.
The first pair (referred to as $Z_1$) is selected from all the possible $\mu^+\mu^-$ pairs to have the smallest mass difference between the $Z_1$ mass and the $Z$ mass, $|m_Z - m_{Z_1}|$.
The second $\mu^+\mu^-$ pair is selected from the remaining muons that has the highest invariant mass (referred to as $Z_2$).
The correct signal di-muon pairing fraction varies with the \Zp mass, where the selected di-muon pair that forms $m_{Z1}$ or $m_{Z2}$ originates from the \Zp.
For example, for $m_{Z'}$ = 5, 42, 63, 72, and 81 GeV, the correct di-muon pairing fractions are about 78\%, 50\%, 88\%, 82\%, and 90\%, respectively.
Finally, the selected four muons must have an invariant mass in a range of 80 to 180 GeV, but excluding the Higgs boson resonance mass region of 110 to 130 GeV.

The $Z'$ signal efficiency at various stages of the event selection is shown in Table~\ref{tab:cutgflow} for five representative mass points.
The event selection efficiencies vary significantly depending on the $Z'$ mass.
At generator level,
an MC filter is applied, which requires at least four muons with $p_T > 2$~GeV and $|\eta| < $ 3.0.
The MC filter efficiencies of these representative $Z'$ signal samples are listed in the table as well.

\begin{table}[!h]
\caption{
The $Z'$ signal event selection efficiencies compared to the events passing the previous cut level for several representative mass points.
The overall signal efficiencies are the products of the 4$\mu$ MC filter and the combined
event selection efficiencies.
}
\centering
\begin{tabular}{l|ccccc}
\toprule
$m_{Z'}$ [GeV]                                    & 5      & 42      & 63     & 72    & 81 \\
\midrule
\midrule
MC filter efficiency                   & 32.8\% & 57.7\%  & 61.0\% & 65.3\% & 70.0\% \\
\midrule
Number of identified muons $\geq$ 4               & 47.3\% & 74.1\%  & 70.8\% & 72.4\% & 75.4\% \\
$p_{\textup{T}}^i (i=1,4) > 20, ~15, ~8, ~3$~GeV  & 60.0\% & 82.6\%, & 90.3\% & 93.6\% & 98.2\% \\
$\Delta R(\mu_i,\mu_j)>0.2$ \& vertex requirement & 87.2\% & 95.4\%  & 96.2\% & 96.6\% & 97.2\% \\
Isolation                                         & 54.2\% & 76.9\%  & 79.2\% & 84.1\% & 87.5\% \\
$m_{4\mu}$ within [80, 110] or [130, 180] GeV     & 91.9\% & 88.8\%  & 58.9\% & 33.5\% & 16.8\% \\
\midrule
Combined event selection efficiency               & 12.3\% & 39.9\%  & 28.7\% & 18.4\% & 10.6\%\\

\midrule
Overall 4$\mu$ signal efficiency                  & 4.1\%  & 23.0\%  & 17.5\% & 11.9\% & 7.4\% \\
\bottomrule
\end{tabular}
\label{tab:cutgflow}
\end{table}

The $Z'$ production signature is searched for in the $Z_1$ or $Z_2$ mass spectrum depending on the $Z'$ mass.
The relatively high-mass $Z'$ signals mostly appear as a peak in the $Z_1$ spectrum while the relatively low-mass signals mostly appear as a peak in the $Z_2$ spectrum.
Representative examples of the predicted signal over background, after further selection with a deep learning approach which will be described in \Sect{\ref{sec:dnn}},
are shown in \Fig{\ref{fig:bkg-only-Z1-Z2}}.
In the analysis the $Z_1$ and $Z_2$ mass spectra are scanned to search for a \Zp with mass greater or smaller than 42~\GeV, respectively.
The boundary value of 42~GeV is chosen based on the studies to optimize the search sensitivity.

The numbers of 4$\mu$ events in data and the estimated background yields are given in \Tab{\ref{tab:background-estimate}}.
More details about  the estimation of the reducible backgrounds containing non-prompt muons can be found in \Sect{\ref{sec:background}}.
The total uncertainties of simulated backgrounds are also listed in the table.
The evaluations of systematic uncertainties will be described in \Sect{\ref{sec:syst}}.

\begin{table}[htb]
\caption{The selected 4$\mu$ events in data and the estimated backgrounds and their  combined statistical and systematic uncertainties.}
\centering{
\begin{tabularx}{\textwidth}{c|c|cccc}
\toprule
Data               & Total      & $qq\rightarrow ZZ^*$ & $gg\rightarrow ZZ^*$ & $ttV+VVV+H$ & Reducible background \\
& background & from MC              & from MC              & from MC     & from data            \\
\midrule
1131               &
$1148 \pm 70$      &
$1065^{+70}_{-69}$ &
$15.6 \pm 2.5$     &
$6.2 \pm 2.9$      &
$61.1 ^{+8.3}_{-9.1}$                                                                                \\
\bottomrule
\end{tabularx}
}
\label{tab:background-estimate}
\end{table}

\section{Estimation of Reducible Background from Data }
\label{sec:background}

A data-driven technique is used to estimate the reducible background from \Zjet and \ttbar production.
Events from these processes may contain two prompt leptons from $W$ or $Z$ boson decays, together with additional activity such as heavy-flavor jets or misidentified components of jets yielding reconstructed muons, collectively referred to as non-prompt muons.
Such backgrounds are estimated from data using a control sample of $\mu^+\mu^- \mu_j \mu_j$ events, selected with the standard signal requirements except that non-prompt muons, $\mu_j$, are selected in place of two of the signal muons. Non-prompt muons are defined as muon candidates using the standard selection but fail the requirements on isolation and impact parameter. However, they are required to pass a much looser impact parameter significance requirement, $|{d_0 \over \sigma_{d_0}} |< 10$.
The background in the signal sample is estimated by scaling each event in the selected $\mu^+\mu^- \mu_j \mu_j$ control sample by $f_1 \times f_2$, where the $f_i (i=1, 2)$ is referred to as a "fake-factor" for each of the two non-prompt muons and is computed as a function of their $p_{\textup{T}}$ and $|\eta|$.

The fake-factor $f$ is derived from two independent non-prompt-muon enriched data control samples dominated by \Zjet or \ttbar events.
The \Zjet (\ttbar) sample contains non-prompt muons dominantly from light jets (heavy-flavor jets).
The \Zjet sample is tagged by two isolated prompt leptons from \Zboson decays ($Z\rightarrow \mu^+\mu^-, ~ e^+e^-$).
The \pT-threshold for the two prompt leptons are 27 and 25 GeV and their invariant mass must be within 10 GeV of the \Zboson mass.

The \MET of the event must be less than 25 GeV.
The \ttbar sample is tagged by two prompt isolated high- \pT $e^\pm \mu^\mp$ leptons with \pT thresholds of 28 and 25 GeV and associated with at least one $b$-jet.
The event is required to have $\MET >$ 50 GeV.
If there is an additional isolated high-\pT lepton in the event, the transverse mass calculated with the \MET and the third lepton is required to be less than 60 GeV to reduce the contribution from the $t\bar t W$ events.
Events with four muons that satisfy the signal selection criteria are removed from both \Zjet and \ttbar control samples.
Additional reconstructed muon objects, passing either the signal-muon or the non-prompt-muon selection criteria are counted in the control samples.
The fake-factor $f$ is calculated as the ratio of the probability for a non-prompt muon to satisfy the signal muon selection criteria to the probability for a non-prompt muon to fail the signal muon selection.
The factor $f$ is derived as a function of $|\eta|$ and \pT of the non-prompt muons,
which varies from 0.11 to 0.76 (0.01 to 0.27) in the \Zjet (\ttbar) control sample.
Then, the fake-factors derived from two data control samples are combined using the $m_{Z1}$ spectra obtained from simulated $4\mu$ background events from \Zjet and \ttbar processes.
A simultaneous fit of the $m_{Z1}$ spectra of the MC events to data selected in the $\mu^+\mu^- \mu_j \mu_j$ control sample is performed to determine the fractions of the \Zjet and \ttbar events in each bin of the mass spectrum.

The overall systematic uncertainty of $f$ is about 14\%.
It is determined with alternative non-prompt muon selections, such as changing the muon isolation criteria in data control samples (7.6\%), and the uncertainties of the \Zjet and \ttbar event fractions when combining the fake-factors derived from two control samples (8.6\%).
The statistical uncertainties of the control samples are also accounted as part of the systematic uncertainties (6.1\%).

Additional reducible background contributions come from the $WZ$ process. These events contain three prompt muons from $W$ and $Z$ decays and one non-prompt muon. This background is estimated by scaling a $3\mu +\mu_j$ control sample by the $f$ derived from the \Zjet sample.
The total estimated number of reducible background events,
$61.1^{+8.3}_{-9.1}$,
is listed in Table~\ref{tab:background-estimate}.

\section{Event Classification with Deep Learning Approach}
\label{sec:dnn}

The selected 4$\mu$ events are classified with a "Deep Learning" approach to further separate the $Z'$ signal from the SM background.
A parametrized deep neural network (pDNN) architecture~\cite{Baldi_2016} is used in the analysis.
This algorithm allows the training of a single classifier for multiple signal mass hypotheses in the search range.
In practice, the kinematic inputs together with the $Z'$ mass parameters of signal and background are used for training.
The MC generated $Z'$  masses (listed in Table~\ref{tab:signalSamples}) are used as the multiple signal mass parameters,
while the distribution of the mass parameter for background events is randomly drawn from the same distribution as for the signal class.
The algorithm was implemented in the Keras \cite{chollet2015keras} framework with the TensorFlow \cite{tensorflow2015-whitepaper} backend.
Two classifiers are trained for low (high) $Z'$ mass searches using MC mass parameters lower (higher) than 40 GeV.
During the training of the classifier, the training samples were composed of simulated signal and background 4$\mu$ events using the pre-selection described in Section~\ref{sec:simu}.
A set of kinematic distributions were used for pDNN training input features.
They are the \pT and $\eta$ of each muon, the \pT of the $Z_1$ and $Z_2$, the mass difference of the $Z_1$ and $Z_2$, $\Delta R$ and $\Delta\eta$ of each muon pair that forms the $Z_1$ and $Z_2$, and the \pT and mass of the $4\mu$ system.
Examples of the input variable distributions, both for representative signals ($m_{Z'}$ = 15, and 51 GeV) and background, are shown in \Fig{\ref{fig:Z1Z2}} to compare the predicted and observed
transverse momenta ($p_{\textup{T}}^{Z_1}$ and $p_{\textup{T}}^{Z_2}$) distributions
and the mass difference of the $Z_1$ and $Z_2$.
In addition to the major background from the SM $Z(Z^*)\rightarrow 4\mu$ production, other backgrounds, including 4$\mu$ events containing non-prompt muons estimated from data, and from $ttV$, $VVV$, and Higgs boson production processes, are included in the plots.

\begin{figure}[!htbp]
\centering

\subfloat[]{
\includegraphics[width=0.33\textwidth]{fig_02a.pdf}
}
\subfloat[]{
\includegraphics[width=0.33\textwidth]{fig_02b.pdf}
}
\subfloat[]{
\includegraphics[width=0.33\textwidth]{fig_02c.pdf}
}
\caption{Distributions of $p_{\text{T}}^{Z_1}$ (a), and $p_{\text{T}}^{Z_2}$ (b), and
the mass difference of the $Z_1$ and $Z_2$ candidates (c).
Small background contributions are denoted as "other backgrounds", including 4$\mu$ events containing non-prompt muons estimated from data and from $ttV$, $VVV$, and Higgs boson production processes.
}
\label{fig:Z1Z2}
\end{figure}

The training results are a set of weights of the pDNN model which are applied to real data or MC samples to obtain the pDNN output discriminating variable.
The $m_{Z1}$ or $m_{Z2}$ are used as the mass parameter when applying the model to data.

\begin{figure}[!hb]
\centering
\subfloat[low mass: $m_{Z'} \leq 42~\GeV$]{
\includegraphics[width=0.48\textwidth]{fig_03a.pdf}
\label{fig:pnn_score:low}
}
\subfloat[high mass: $m_{Z'} \geq 42~\GeV$]{
\includegraphics[width=0.48\textwidth]{fig_03b.pdf}
\label{fig:pnn_score:high}
}
\caption{
The pDNN output discriminant variable distributions for low mass (a) and high mass (b) with a signal sample at 35 \GeV ~and 51 \GeV, respectively.
}
\label{fig:pnn_score}
\end{figure}

To improve the power of the pDNN classifier the model is optimized through an automatic process developed with the package \textsc{Tune}~\cite{liaw2018tune} and \textsc{HEBO Search} algorithm~\cite{alex2020hebo}.
The neural-network structure was chosen to have two fully connected hidden layers, each with 256 (32) nodes for the high (low) mass search in the analysis.
Other training hyper-parameters, such as the \textit{learning rate} with \textit{decay} and \textit{class weight}, for both high and low mass searches were also selected.
Two pDNN output discriminant classifiers for high
and low mass
searches were obtained and are shown in Figure~\ref{fig:pnn_score} for data and simulated signal and background events.

To maximize the search sensitivity, a scan of the pDNN output scores was performed to find the optimal cut values for each \Zp mass hypothesis.
These cut values vary from 0.42 to 0.74 (0.12 to 0.16) for the low (high) mass region.
These cuts keep high signal efficiencies between 98\% and 95\%  (90\% and 50\%), while the corresponding background reductions range from 10\% to 50\%  (50\% to 96\%) for the low (high) mass region.

After the 4$\mu$ event selection with the pDNN classifier, the final discriminant to search for the $Z'$ resonance signature is the $Z_1$ (for $m_{Z'} \geq 42$~GeV) or $Z_2$ (for $m_{Z'} \leq 42$~GeV) mass spectrum as shown in \Fig{\ref{fig:bkg-only-Z1-Z2}}.
Data are compared to the estimated background together with two representative signals with masses of 15 and 51 GeV, shown in \Fig{\labelcref{fig:bkg-only-Z1-Z2:a,fig:bkg-only-Z1-Z2:b}}, respectively.
The values of the gauge coupling strengths ($g$) for the two mass points are chosen for the purpose of illustration.

One should note that the dimuon masses ($m_{Z1},~ m_{Z2}$) are not used  as the pDNN training variables, therefore cutting on the pDNN output scores does not sculpt the dimuon invariant mass distributions of the backgrounds (see \Fig{\ref{fig:bkg-only-Z1-Z2}}).
To interpret the observed mass spectra systematic uncertainties in the predictions are estimated in Section~\ref{sec:syst}.

\begin{figure}[!htb]
\centering
\subfloat[]{
\includegraphics[width=0.48\textwidth]{fig_04a.pdf}
\label{fig:bkg-only-Z1-Z2:a}
}
\subfloat[]{
\includegraphics[width=0.48\textwidth]{fig_04b.pdf}
\label{fig:bkg-only-Z1-Z2:b}
}
\caption{
Mass spectra of $m_{Z2}$ (left) and $m_{Z1}$ (right) for the pDNN-selected events with a signal sample at 15 \GeV ~and 51 \GeV, respectively.
}
\label{fig:bkg-only-Z1-Z2}
\end{figure}

\section{Systematic Uncertainties}
\label{sec:syst}

Systematic uncertainties in the simulated event yields and shapes, for both signal and background processes, may arise from the calibration of the physics objects and from the theoretical modelling used in the predictions.

The major experimental uncertainties come from the muon reconstruction, identification, and isolation requirement efficiencies.
These efficiencies are corrected based on studies performed in data control regions.
The energy and momentum scales and resolutions of the simulated objects are calibrated to reproduce data from $Z\rightarrow \mu^+\mu^-$ and $J/\psi \rightarrow \mu^+\mu^-$ decays~\cite{MUON-2018-03-updated}.
The uncertainties on the 4$\mu$ detection efficiency are determined
by varying the nominal calibrations in the MC samples by one standard deviation, including muon momentum resolutions and scales, and the trigger, reconstruction, identification, and isolation requirement efficiencies.
The overall relative experimental uncertainties in the 4$\mu$ event selection efficiency is about 3.9\%, dominated by the uncertainty of the isolation efficiency (2.9\%) and the low $p_{\textup{T}}$ calibration uncertainty (2.0\%).
The signal event selection efficiency uncertainties vary from 8.3\% to 3.9\%, depending on the \Zp mass.
In addition, the uncertainty in the combined 2015-2018 integrated luminosity is 1.7\%~\cite{ATLAS-CONF-2019-021}, obtained using the LUCID-2 detector~\cite{Avoni_2018} for the primary luminosity measurements.

Sources of theoretical uncertainties come from the choice of QCD scales (renormalization $\mu_R$ and factorization $\mu_F$), strong coupling constant $\alpha_S$, and PDFs, as well as the parton shower models.
These uncertainties affect the signal and background event selection efficiencies, normalization, and the shape of their kinematic distributions.
The scales are varied independently from 0.5 to 2.0 times the nominal values and the largest deviation is chosen as the systematic uncertainty.
The PDF uncertainty is estimated by comparing events generated with different PDF sets, as well as the uncertainties from the nominal PDF set itself.
The maximal variation (envelope) is accounted for as the systematic uncertainty, following the PDF4LHC~\cite{Butterworth:2015oua} recommendations.
The $\alpha_{S}$ uncertainty is estimated by comparing events generated with different $\alpha_{S}$ values using the nominal PDF set.
The parton shower uncertainty is estimated by comparing events with different parton shower parameters in the Sherpa MC samples.
For the $Z(Z^*)\rightarrow 4\mu$ process, the relative uncertainties of event yields for scale, PDF, $\alpha_s$, and parton shower,  are 4.6\%, ~1.8\%, ~1.0\%, and 1.9\%, respectively.
The $t\bar t V$,
$VVV$,
and Higgs boson,
processes contribute
5.2, 0.50, and 0.53 background events, respectively.
The total relative theoretical uncertainty for these background is estimated about 16\%.
The interference effect between the signal and background is estimated (using MC samples generated by \MGNLO)
by evaluating the cross-section ratio, ${\Delta \sigma} / {\sigma_{Z'}}$,
where $\Delta \sigma$ is
the difference of the cross-section from the inclusive MC sample and the sum of the cross-sections of the signal and background samples in the $Z'$ detection phase space.
The effect varies from 1.8 to 7.5\% for the $Z'$ mass from 5 - 81 GeV,
which is accounted for as additional signal yield systematic uncertainties.
In addition, the MC $Z'$ signal filter acceptance uncertainty is estimated to be about 2\%, which
is calculated by varying the QCD scales, the PDF-sets, and the strong coupling constant using MC events at the generator level for different $Z'$ mass points.

All the uncertainties (from both experiment and theory) on the final discriminant di-muon mass spectra are included as nuisance parameters in the signal extraction fitting process described in Section~\ref{subsec:final_result}.

\section{Data Interpretation and Results}
\label{sec:result}

The statistical analysis is performed by comparing the data to
the sum of the background prediction and the signal
to search for the $Z'$ signature and information about the $pp\rightarrow \mu^+\mu^-Z'\rightarrow 4\mu$ signal production cross-section and the associated coupling strength.
In case of no significant data excess over the background prediction, upper limits on the signal production cross-section times the decay branching fraction for different \Zp masses are set at 95\% confidence level (CL).

\subsection{Statistical Procedure and Results}
\label{subsec:final_result}

The statistical procedure used to interpret the data is described in \Refn{\cite{HIGG-2012-17}}.
The likelihood function is constructed from the product of Poisson probabilities:
\begin{equation*}
L (\text{data} |~\mu, \theta)= \prod_{i=1}^N \text{Poisson}(\text{data}_i |~\mu \cdot s_i(\theta) + b_i(\theta))\times {\text G}(\tilde\theta |\theta)
\end{equation*}
where $i$ is the bin-index of the fitted variable distribution;
$\mu$ is the signal strength, and $\theta$ denotes the nuisance parameters, which represent the uncertainties of the measurements;
${\text G}(\tilde\theta |\theta)$ is the Gaussian scaling function
of the nuisance parameters constructed as deviations from the nominal model of the systematic uncertainties, where $\tilde\theta$ provides a maximum likelihood estimate for $\theta$.
The parameter of interest in the statistical analysis is the global signal strength factor $\mu$, which acts as a scale factor on the total number of events predicted by the signal model.
This factor is defined such that $\mu = 0$ corresponds to the background-only hypothesis and $\mu > 0 $ corresponds to a \Zp signal in addition to the background.
Hypothesised values of $\mu$ are tested based on the profile likelihood ratio~\cite{Cowan:2010js}, which compares data with background-only ($b$) and signal+background ($s+b$) models using the following test statistic:
\begin{equation*}
q_\mu =
\begin{cases}
-2 \text{ln} \frac{L( \text{data} |~\mu, \hat{\hat{\theta}}_\mu)}{L(\text{data} |~\hat\mu, \hat\theta)} & \hat{\mu} \leq \mu \\
0                                                                                & \hat{\mu} > \mu    \\
\end{cases}
\end{equation*}
where $\hat\mu ~\text{and} ~\hat\theta$ are the values
for $\mu$ and $\theta$ when maximizing the $L$ with all parameters floating, which are named as the unconditional maximum-likelihood (ML) estimators;
$\hat{\hat{\theta}}_\mu$ is the conditional ML estimator of $\theta$ for a fixed value of $\mu$.
This test statistic extracts the information on the signal strength from a full likelihood fit to the data.
The likelihood function includes all the parameters that describe the systematic uncertainties and their correlations.

In case of no significant data excess over the background, exclusion limits are set, based on the CLs prescription~\cite{Read:2002hq}, which calculates the ratio of the probabilities based on the $s+b$ background-only models:

\begin{equation*}
\text{CL}_s(\mu) = \frac{\text{CL}_{s+b}}{\text{CL}_{b}}
= \frac{\text{P}(q_\mu> q_\mu^{obs})}{\text{P}(q_0 > q_0^{obs})}
\end{equation*}
A value of $\mu$ is regarded as excluded at 95\% confident level when CLs is less than 5\%.

The significance of an excess in the data is first quantified
with the local $p_0$, the probability that the background
can produce a fluctuation greater than or equal
to the excess observed in data.
The equivalent formulation in terms of number of standard deviations is
referred to as the local significance, which is computed from test statistic $q_0$:
\begin{equation*}
p_0 = \text{P}(q_0 < q_0^{obs}).
\end{equation*}
The global probability
for the most significant excess to be observed anywhere in a given search region is estimated with the method described in Ref.~\cite{Gross_2010}.

The statistical tests are performed by first searching for data excesses by means of a $p_0$ scan, and in case no such significant data excess is found, by setting limits on the \Zp production cross-section times the branching fraction.
In both steps the values and reconstructed mass distributions of hypothesized \Zp mass $m_{Z'}$ are used.
The steps between the mass points used in these test procedures are set to match the mass resolutions determined by reconstruction of the di-muon invariant mass.
A linear interpolation of the expected event yields and shapes between the MC generated signal samples is used in the fitting process.
The asymptotic approximation~\cite{Cowan:2010js} upon which the results are based has been validated against the method (detailed in \Refn{\cite{HIGG-2012-17}}) of using pseudo-experiments for several mass points.

In both steps of the statistical tests, data are fit to the $m_{Z1}$ and $m_{Z2}$ spectrum with background ($b$) only and signal+background ($s+b$) hypotheses.
In the fitting process each mass spectrum is divided in the signal region (SR) and the background control region (CR).
For each $Z'$ mass point the SR is defined in a mass window of $m_{Z'} \pm 3\sigma_{m_{\mu\mu}}$ of the di-muon mass spectrum.
The sidebands outside of the SR are defined as the CR.
The di-muon mass resolution $\sigma_{m_{\mu\mu}}$ is determined by the fully simulated $Z'$ mass distribution, which combines the \Zp natural width and the detector resolution, ranging from 0.10 to 1.75 GeV.
The mass resolution is mostly dominated by the detector resolution.
Finer binning is used in the SR to enhance the sensitivity.
The background CR is used to constrain the overall normalization for the background in the signal region.
The shape of the major background from $Z(Z^*)\rightarrow 4\mu$ is fixed with prior uncertainties included in the fitting process, but the normalization (or strength) floats in the fit.
Other background normalizations and shapes are fixed with prior uncertainties included in the fitting.

\subsection{The $p_0$ Scan Results}
\label{subsec:p-scan}

The $p_0$-values corresponding to the background-only hypothesis are scanned in the mass range of this analysis.
A binned profile-likelihood fit \cite{Cowan:2010js} is performed simultaneously across the $Z'$ signal-region and the background control region using the predicted and observed mass spectrum as inputs.
Data are fit to the $m_{Z1}$ and $m_{Z2}$ distributions
for $m_{\Zp} \geq 42~\GeV$ and $m_{\Zp} \leq 42~\GeV{}$, respectively,
with a "sliding" mass window as the defined SR changes for different \Zp mass points.
The chosen bin-size inside the SR is around 0.3~$\sigma_{m_{\mu\mu}}$, for each mass point. The total number of bins in the CR is 20.
The fit mass range of $m_{Z1}$ ($m_{Z2}$) is [30, 85] GeV ([0, 45] GeV).

The $p_0$-values at different \Zp mass hypothesis points are
computed and transformed into Gaussian standard deviations to indicate the significance as
shown in \Fig{\ref{fig:p_value_scan}}.
The smallest $p_0$-value is at 39.6 GeV, corresponding to a local 2.65$\sigma$ deviation from the background-only hypothesis, while the global deviation~\cite{Gross_2010} is found to be 0.52$\sigma$, indicating that no significant data excess over the expected background is observed.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{fig_05.pdf}
\caption{The $p_0$-value scan across the \Zp mass signal regions.}
\label{fig:p_value_scan}
\end{figure}

\subsection{Upper Limits}
\label{subsec:limits}

The upper limits on the production cross-section times branching fraction of the $pp\rightarrow \mu^+\mu^-Z'\rightarrow 4\mu$ process are calculated using a similar fitting procedure described in Section~\ref{subsec:p-scan}.
Confidence intervals are computed based on the profile-likelihood-ratio test statistics~\cite{Cowan:2010js}.
The observed and expected upper limits at 95\% CL on the cross-section times branching fraction, $\sigma (pp\rightarrow Z'\mu\mu\rightarrow 4\mu)$, are shown in \Fig{\ref{fig:xs_limit}}.
Assuming the same coupling strength \coupzp of the $Z'$ boson to the second and third lepton families and to the left- and right-handed fermions, the upper limits on the coupling parameter \coupzp are extracted from the limits of the $Z'$ production cross-section times branching fraction using the $L_\mu-L_\tau$ model, where the branching fraction of $\mathcal{B}$($Z'\rightarrow \mu^+\mu^-) = {1\over 3}$, which is determined by counting all the possible \Zp decay modes in this model.
At each generated $Z'$ mass point, a limit on the coupling strength \coupzp has been obtained from the cross-section limit.
The observed and expected upper limits on the coupling parameter \coupzp are shown in \Fig{\ref{fig:coup_limit}}.
The limits on the coupling \coupzp are in the range of 0.003 (for $m_{Z'}$ = 5~GeV) to 0.2 (for $m_{Z'}$ = 81 GeV) depending on the \Zp mass ranging from 5 to 81 GeV.
This ensures that the ratio of the $Z'$ natural width and mass, $\Gamma(Z') / m_{Z'}$, is well below 1\% in this mass range.
\begin{figure}[!htbp]
\centering
\subfloat[]{
\includegraphics[width=0.49\textwidth]{fig_06a.pdf}
\label{fig:xs_limit}
}
\subfloat[]{
\includegraphics[width=0.49\textwidth]{fig_06b.pdf}
\label{fig:coup_limit}
}
\caption{95\% CL upper limits (expected and observed) on the cross-sections times branching fraction (a) and coupling parameter (b). The discontinuity at 42 GeV represents the border of the low/high mass classifiers.
}
\label{fig:limit}
\end{figure}

Motivated by theoretical interpretations in \Refn{~\cite{Altmannshofer:2016jzy}},
a 2-dimensional exclusion contour at 95\% CL in the parameter-space of $(m_{Z'}, \coupzp)$ of the $L_\mu-L_\tau$ model from this analysis is produced and shown in \Fig{\ref{fig:limit_comp}}.
The parameter space exclusion regions calculated by theorists using data from
the Neutrino Trident experiment~\cite{PhysRevLett.66.3117}
and the $B_s$ mixing measurements by a global analysis performed in \Refn{~\cite{Altmannshofer:2016jzy}}
are also shown in \Fig{\ref{fig:limit_comp}}.
This had left a large gap in the parameter space not yet excluded, that was allowed to explain the LHCb $b\rightarrow s \mu^+\mu^-$ anomalies~\cite{Aaij2022,Aaij_2013}.
This gap is now largely excluded by this analysis.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{fig_07.pdf}
\caption{
The coupling parameter $g$ limits from this search as a function of the $Z'$ mass
compared to the limits~\cite{Altmannshofer:2016jzy} from the Neutrino Trident (red) and the $B_s$ mixing (green) experimental results.
}
\label{fig:limit_comp}
\end{figure}

\FloatBarrier

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

A search for a new vector gauge boson $Z'$ predicted by the $L_\mu-L_\tau$ models has been performed with a 4$\mu$ final state in the invariant mass range of [80, 180]~\GeV{}, excluding the Higgs boson mass window [110, 130]~\GeV{}, using 139~\ifb{} of $\sqrt{s} = $ 13~\TeV{} proton-proton collision data collected with the ATLAS detector.
No significant excess of events over the expected SM background is observed.
Therefore, upper limits are set on the \Zp production cross-section times the decay branching fraction
of the $pp \rightarrow Z'\mu^+\mu^- \rightarrow \mu^+\mu^-\mu^+\mu^-$ process,
varying from 0.31 to 4.3 fb at 95\% CL,
in a $Z'$ mass range of [5, 81]~\GeV{},
from which the coupling strength \coupzp of the \Zp to muons above 0.003 to 0.2 (depending on the \Zp mass) are excluded in the same mass range.
This search shows significant sensitivity improvements over previous indirect and direct searches of the \Zp with 4$\mu$ final state.
An interesting parameter space of the $L_\mu-L_\tau$ model prediction
that was not excluded before and could explain the $b\rightarrow s \mu^+\mu^-$ anomalies is now largely excluded by this search.

\clearpage

\clearpage
\printbibliography

\clearpage

\end{document}


\section{Introduction}

Graph editing problems are fundamental problems in graph theory. They define a set of basic operations on a graph $G$ and ask for the minimum number of these operations necessary in order to turn $G$ into a graph of a desired target graph class $\mathcal G$~\cite{nss-ccsemp-01,ly-nphpn-80,y-nenp-78,k-aog-96}.
For instance, in the Cluster Editing problem~\cite{sst-cgmp-04} the operations are insertions or deletions of individual edges and the target graph class are cluster graphs, i.e., unions of vertex-disjoint cliques.
In graph drawing, a particularly interesting graph class are planar graphs,
for which several related graph editing problems have been studied, e.g., how many vertex deletions are needed to turn an arbitrary graph into a planar one~\cite{ms-opgvd-12} or how many vertex splits are needed to obtain a planar graph~\cite{jackson1984splitting,Faria_2001}.
In this paper, we are interested in the latter operation: vertex splitting.
A \emph{vertex split} creates two copies of a vertex $v$, distributes its edges among these two copies and then deletes $v$ from~$G$.

Further, we are translating the graph editing problem into a more geometric or topological drawing editing problem.
This means that we apply the splitting operations not to the vertices of an abstract graph, but to the vertices of a planar graph drawing, or more generally to a planar embedded (or \emph{plane}) graph.
In a plane graph, each vertex has an induced cyclic order of incident edges, which needs to be respected by any vertex split in the sense that we must split its cyclic order into two contiguous intervals, one for each of the two copies.
From a different perspective, the two faces that serve as the separators of these two edge intervals are actually merged into a single face by the vertex split.

Finally, we consider outerplanar graphs as the target graph class.
Thus, we want to apply a minimum number of vertex splits to a plane graph $G$, which merge a minimum number of faces in order to obtain an outerplanar embedded graph $G'$, where all vertices are incident to a single face, called the \emph{outer face}. We denote this minimum number of splits as the \emph{outerplane splitting number} $\ensuremath{\operatorname{osn}}(G)$ of $G$ (see~\cref{fig:introfig}).
Outerplanar graphs are a prominent graph class in graph drawing (see, e.g., \cite{Biedl11,Frati22,LazardLL19,LenhartL96prox}) as well as in graph theory and graph algorithms more generally (e.g., \cite{BodlaenderF02,f-saicrt-96,mz-emaog-99}).
For instance, outerplanar graphs admit planar circular layouts or 1-page book embeddings~\cite{bk-btg-79}.
Additionally, outerplanar graphs often serve as a simpler subclass of planar graphs with good algorithmic properties. For instance, they have treewidth 2 and
their generalizations to $k$-outerplanar graphs still have bounded treewidth~\cite{Bodlaender98,Biedl_2015}, which allows for polynomial-time algorithms for \NP-complete problems that are tractable for such bounded-treewidth graphs. This, in turn, can be used to obtain a PTAS for these problems on planar graphs~\cite{Baker94}.

We are now ready to define our main computational problem as follows.

\begin{problem}[\textsc{Outerplane Splitting Number}\xspace]
	Given a plane biconnected graph $G=(V,E)$ and an integer $k$, can we transform $G$ into an outerplane graph $G'$ by applying at most $k$ vertex splits to $G$?
\end{problem}

\begin{figure}[t]
\centering
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[page=3]{intro.pdf}
    \caption{}
    \label{fig:subfig3}
\end{subfigure}
\hfil
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[page=4]{intro.pdf}
    \caption{}
    \label{fig:subfig3}
\end{subfigure}
\caption{(a)~An instance of \textsc{Outerplane Splitting Number}\xspace, where the colored vertices will be split; (b) resulting outerplane graph after the minimum 3 splits.}
\label{fig:introfig}
\end{figure}

\paragraph{Contributions.}
In this paper, we introduce the above  problem \textsc{Outerplane Splitting Number. We start by showing} the key property for our subsequent results, namely that (minimum) sets of vertex splits to turn a  plane biconnected graph $G$ into an outerplane one correspond to (minimum) connected face covers in $G$ (\cref{sec:fvig}), which in turn are equivalent to (minimum) feedback vertex sets in the dual graph of $G$.
Using this tool we then show that for general  plane biconnected graphs \textsc{Outerplane Splitting Number}\xspace is \NP-complete (\cref{sec:hardness}), whereas for maximal planar graphs we can solve it in polynomial time (\cref{sec:fvs}).
Finally, we provide upper and lower bounds on the outerplane splitting number for maximal planar graphs (\cref{sec:bounds}).

\paragraph{Related Work.}
Splitting numbers have been studied mostly for abstract (non-planar) graphs with the goal of turning them into planar graphs.
The \textsc{Planar Splitting Number} problem is \NP-complete in general~\cite{Faria_2001}, but exact splitting numbers are known for complete and complete bipartite graphs~\cite{HartsfieldJR85,jackson1984splitting}, as well as for the 4-cube~\cite{FariaFN98}. For two-layer drawings of general bipartite graphs, the problem is still \NP-complete, but \FPT~\cite{splittingbipfpt} when parametrized by the number of split vertices.
It has also been studied for other surfaces such as the torus~\cite{Hartsfield86} and the projective plane~\cite{Hartsfield87}.
Another related concept is the split thickness of a graph $G$ (or its folded covering number~\cite{ku-twcg-16}), which is the smallest $k$ such that $G$ can be transformed into a planar graph by applying at most $k$ splits per vertex.
Recognizing graphs with split thickness $2$ is \NP-hard, but there is a constant-factor approximation algorithm and a fixed-parameter algorithm for graphs of bounded treewidth~\cite{Eppstein_2017}.
Recently, the complexity of the embedded splitting number problem of transforming non-planar graph drawings into plane ones has been investigated~\cite{splittingnoellenburg}.
Beyond the theoretical investigations of splitting numbers and planarity, there are also applied work in graph drawing making use of vertex splitting to untangle edges~\cite{Wu_2020_TVCG} or to improve layout quality for community exploration~\cite{admpt-hgvwcaea-22,HenryBF08}.

Regarding vertex splitting for achieving graph properties other than planarity, Trotter and Harary~\cite{th-dmig-79} studied vertex splitting to turn a graph into an interval graph.
Paik et al.~\cite{PaikRS98} considered vertex splitting to remove long paths in directed acyclic graphs and
Abu-Khzam et al.~\cite{Abu-KhzamBFS21} studied heuristics using vertex splitting for a cluster editing problem

\begin{figure}[t]
\centering
\begin{subfigure}[b]{.24\textwidth}
    \centering
    \includegraphics[page=1]{split_intro}
    \caption{}
    \label{fig:split_intro_1}
\end{subfigure}
\hfil
\begin{subfigure}[b]{.24\textwidth}
    \centering
    \includegraphics[page=2]{split_intro}
    \caption{}
    \label{fig:split_intro_2}
\end{subfigure}
\begin{subfigure}[b]{.24\textwidth}
    \centering
    \includegraphics[page=3]{split_intro}
    \caption{}
    \label{fig:split_intro_3}
\end{subfigure}
\begin{subfigure}[b]{.24\textwidth}
    \centering
    \includegraphics[page=4]{split_intro}
    \caption{}
    \label{fig:split_intro_4}
\end{subfigure}
\caption{(a)~Two touching faces $f_1, f_2$ with a common vertex $v$ on their boundary. (b)~Result of the split of $v$ with respect to $f_1, f_2$ joining them into a new face $f$. (c-d)~Merging 4 faces $f_1, \ldots, f_4$ covering a single vertex $v$ with 3 splits.}
	\label{fig:split_intro}
\end{figure}

\paragraph{Preliminaries.} The key concept of our approach is to merge a set of faces of a given plane graph $G=(V,E)$ with vertex set $V=V(G)$ and edge set $E=E(G)$ into one big face which is incident to all vertices of $G$. Hence, the result is outerplanar.
The idea is that if two faces $f_1$ and $f_2$ share a vertex $v$ on their boundary (we say $f_1$ and $f_2$ \emph{touch}, see \cref{fig:split_intro_1}), then we can split $v$ into two new vertices $v_1, v_2$. In this way, we are able to create a narrow gap, which merges $f_1, f_2$ into a bigger face $f$ (see \cref{fig:split_intro_2}). With this in mind, we formally define an \emph{embedding-preserving split} of a vertex $v$ w.r.t. two incident faces $f_1$ and $f_2$. We construct a new plane graph $G'=(V', E')$ with $V' = V \setminus \{v\} \cup \{v_1,v_2\}$. Consider the two neighbors of $v$ both incident to $f_1$ and
let $w_1$ be the second neighbor in clockwise order.
Similarly, let $w_i$ be the second vertex adjacent to $v$ and incident to $f_2$. We call $w_d$ the vertex preceding $w_1$ in the cyclic ordering or the neighbors, with $d$ being the degree of $v$, see~\cref{fig:split_intro_1}. Note that while $w_1=w_{i-1}$ and $w_i=w_d$ is possible, $w_d\ne w_1$ and $w_{i-1}\ne w_i$.
For the set of edges, we now set $E'= E \setminus \{ (v,w_1), \ldots, (v,w_d) \} \cup \{ (v_2, w_1), \ldots, (v_2, w_{i-1}) \} \cup \{ (v_1, w_i), \ldots, (v_1, w_{d}) \}$ and assume that they inherit their embedding from $G$.
From now on we refer to this operation simply as a \emph{split} or when $f_1,f_2$ are clear from the context, we may refer to \emph{merging} the two faces at $v$.
The vertices $v_1,v_2$ introduced in place of $v$ are called \emph{copies} of $v$.
If a copy $v_i$ of a vertex $v$ is split again, then any copy of $v_i$ is also called a copy of the original vertex $v$.

We can now reformulate the task of using as few splits as possible.
Our objective is to find a set of faces $S$ that satisfies two conditions.
(1) Every vertex in $G$ has to be on the boundary of at least one face $f \in S$, that is, the faces in $S$ \emph{cover} all vertices in $V$.\footnote{Testing whether such $S$ with $|S| \le k$ exists, is the \NP-complete  problem \textsc{Face Cover}~\cite{bm-ccvfpg-88}.}
And (2) for every two faces $f,f' \in S$ there exists a set of faces $\{f_1,\dots,f_k\}\subseteq S$
such that $f=f_1, \ldots, f_k = f'$,
and $f_i$ touches $f_{i+1}$ for $1\leq i < k$.
In other words, $S$ is connected in terms of touching faces. We now introduce the main tool in our constructions that formalizes this concept.

\section{Face-Vertex Incidence Graph}\label{sec:fvig}
Let $G=(V,E)$ be a plane biconnected graph and $F$ its set of faces. The \emph{face-vertex incidence graph\xspace} is defined as $H = (V \cup F, E_H)$ and contains the edges $E_H = \{ (v,f) \in V \times F : v \text{ is on the boundary of } f \}$.
Graph $H$ is by construction bipartite and we
assume that it is plane by placing each vertex $f \in F$ into its corresponding face in~$G$.

\begin{definition}\label{def:cfc}
Let $G$ be a plane biconnected graph, let $F$ be the set of faces of $G$, and let $H$ be its face-vertex incidence graph\xspace.
A \emph{face cover} of $G$ is a set $S \subseteq F$ of faces such that every vertex $v\in V$ is incident to at least one face in $S$.
A face cover $S$ of
$G$ is a \emph{connected} face cover if the induced subgraph $H[S \cup V]$ of $S \cup V$ in $H$ is connected.
\end{definition}
\noindent
We point out that the problem of finding a connected face cover\xspace is not equivalent to the Connected Face Hitting Set Problem \cite{SchweitzerS10}, where a connected set of vertices incident to every face is computed.
\noindent We continue with two lemmas that are concerned with
merging multiple faces at the same vertex (\cref{fig:split_intro_3}).

\begin{lemma}\label{lem:single-vertex-split}
Let $G$ be a plane biconnected graph and $S \subseteq F$ a subset of the faces $F$ of $G$ that all have the vertex $v \in V$ on their boundary. Then $|S|-1$ splits are sufficient to merge the faces of $S$ into one.
\end{lemma}

\begin{proof}
Let $f_1, \ldots, f_k$ with $k = |S|$ be the faces of $S$ in the clockwise order as they appear around $v$ ($f_1$ chosen arbitrarily).
We iteratively merge $f_1$ with $f_i$ for $2 \leq i \leq k$, which requires in total $|S|-1$ splits (see \cref{fig:split_intro_3} and \cref{fig:split_intro_4}). \qed
\end{proof}

\begin{lemma}\label{lem:cover-to-split}
Let $G$ be a plane biconnected graph and let $S$ be a connected face cover\xspace of $G$. Then $|S|-1$ splits are sufficient to merge the faces of $S$ into one.
\end{lemma}

\begin{proof}
Let $H' = H[S \cup V]$ and
compute a spanning tree $T$ in $H'$. For every vertex $v \in V(T) \cap V(G)$, we apply \cref{lem:single-vertex-split} with the face set $F'(v) = \{ f \in S \cap V(T) \; | \; (v,f)  \in E(T) \}$.
We root the tree at an arbitrary face $f' \in S$, which provides a hierarchy on the vertices and faces in $T$.  Every vertex $v \in V(T) \cap V(G)$ requires by \cref{lem:single-vertex-split} $|F'(v)|-1$ splits. Note that that for all leaf vertices in $T$, $|F'(v)|=1$, i.e., they will not be split. Each split is charged to the children of $v$ in $T$. Since $H$ is bipartite, so is $T$. It follows that every face $f \in S \setminus \{ f' \}$ is charged exactly once by its parent, thus $|S|-1$ splits suffice. \qed
\end{proof}

\begin{lemma}\label{lem:split-to-cover}
Let $G$ be a plane biconnected graph and $\sigma$ a sequence of $k$
splits to make $G$ outerplane. Then $G$ has a connected face cover\xspace of size $k+1$.
\end{lemma}

\begin{proof}
    Since by definition applying $\sigma$ to $G$ creates a single big face that is incident to all vertices in $V(G)$ by iteratively merging pairs of original faces defining a set $S\subseteq F$, it is clear that $S$ is a face cover of $G$ and since the result of the vertex splits and face merges creates a single face, set $S$ must also be connected. \qed
\end{proof}

As a consequence of \cref{lem:cover-to-split,lem:split-to-cover} we obtain that \textsc{Outerplane Splitting Number}\xspace and computing a minimum connected face cover\xspace are equivalent.

\begin{theorem}\label{thm:split-cfc}
    Let $G$ be a plane biconnected graph. Then $G$ has outerplane splitting number $k$ if and only if it has a connected face cover of size $k+1$.
\end{theorem}

\section{\NP-completeness}\label{sec:hardness}
In this section, we prove that finding a connected face cover\xspace of size $k$ (and thus \textsc{Outerplane Splitting Number}\xspace) is \NP-complete. The idea is to take the dual of a planar biconnected \textsc{Vertex Cover}\xspace instance and subdivide every edge once (we call this an \emph{all-1-subdivision}). Note that the all-1-subdivision of a graph $G$ corresponds to its vertex-edge incidence graph and the all-1-subdivision of the dual of $G$ corresponds to the face-edge incidence graph of $G$. A connected face cover\xspace then corresponds to a vertex cover in the original graph, and vice versa. The following property greatly simplifies the arguments regarding \cref{def:cfc}.

\begin{prp}\label{prp:subdiv}

	Let $G'$ be an all-1-subdivision of a biconnected planar graph $G$ and $S$ a set of faces that cover $V(G')$. Then $S$ is a connected face cover\xspace of $G'$.
\end{prp}

\begin{proof}
Let $H$ be the all-1-subdivision of the dual of $G$, and assume to the contrary that the induced subgraph $H' = H[S \cup V(G)]$ is not connected. Then there exists an edge $(u,v) \in E(G)$ such that $u$ and $v$ are in different connected components in $H'$.
Let $w$ be the subdivision vertex of $(u,v)$ in $G'$.
As a subdivision vertex, $w$  is incident to only two faces, one of which, say $f$, must be contained in $S$.
But $f$ is also incident to $u$ and $v$ and hence $u$ and $v$ are in the same component of $H'$ via face $f$, a contradiction. Hence $H'$ is connected and $S$ is a connected face cover\xspace of $G'$. \qed
\end{proof}

\begin{figure}[t]
	\centering
	\begin{subfigure}[b]{.32\textwidth}
		\centering
		\includegraphics[page=2]{VC_reduction.pdf}
		\caption{primal (black) and\\dual graph (gray)}
		\label{fig:subfig99}
	\end{subfigure}
	\hfil
	\begin{subfigure}[b]{.32\textwidth}
		\centering
		\includegraphics[page=3]{VC_reduction.pdf}
		\caption{subdivision of\\the dual graph}
		\label{fig:subfig99}
	\end{subfigure}
	\hfil
	\begin{subfigure}[b]{.32\textwidth}
		\centering
		\includegraphics[page=4]{VC_reduction.pdf}
		\caption{vertex cover (red) in\\the primal graph}
		\label{fig:subfig99}
	\end{subfigure}
	\caption
	Link between the primal graph $G$, its vertex cover, the dual $D$ and its subdivision $D^*$.}
	\label{fig:vertex_cover}
\end{figure}

The proof of the next theorem is very similar to the reduction of Bienstock and Monma to show \NP-completeness of \textsc{Face Cover}~\cite{bm-ccvfpg-88}; due to differences in the problem definitions, such as the connectivity of the face cover and whether the input graph is plane or not, we provide the full reduction for the sake of completeness.

\begin{theorem}
Deciding whether a plane biconnected graph $G$ has a connected face cover\xspace of size at most $k$ is \NP-complete.
\end{theorem}

\begin{proof}
Clearly the problem is in $\NP$. To prove hardness, we first introduce some notation. Let $G$ be a plane biconnected graph and $D$ the corresponding dual graph.
Furthermore, let $D^*$ be
the all-1-subdivision of $D$.
We prove now that a connected face cover\xspace $S^*$ of size $k$ in $D^*$ is in a one-to-one correspondence with a vertex cover $S$ of size $k$ in $G$ (see \cref{fig:vertex_cover}). More specifically, we show that the dual vertices of the faces of $S^*$ that form a connected face cover\xspace in $D^*$, are a vertex cover for $G$ and vice versa. The reduction is from the \NP-complete \textsc{Vertex Cover}\xspace problem in biconnected
planar graphs in which all vertices have degree 3 (cubic graphs)~\cite{Mohar01}.

\paragraph{\textsc{Connected Face Cover}\xspace $\Rightarrow$ \textsc{Vertex Cover}\xspace:}
Let $G$ be such a biconnected plane \textsc{Vertex Cover}\xspace instance. Assume we have a connected face cover\xspace $S^*$ with $|S^*|=k$ for $D^*$. Note that the faces of $D^*$ correspond to the vertices in $G$. We claim that the faces $S^*$, when mapped to the corresponding vertices $S \subseteq V(G)$ are a vertex cover for $G$. Assume otherwise, that is, there exists an edge $e^* \in E(G)$ that has no endpoints in $S$. However, $e^*$ has a dual edge $e \in E(D)$ and therefore a subdivision vertex $v_e \in V(D^*)$. Hence, there is a face $f \in S^*$ that has $v_e$ on its boundary by definition of connected face cover\xspace. And when mapped to $D$, $f$ has $e$ on its boundary, which implies that the primal edge $e^*$ has at least one endpoint in $S^*$; a contradiction.

\paragraph{\textsc{Vertex Cover}\xspace $\Rightarrow$ \textsc{Connected Face Cover}\xspace:}
To prove that a vertex cover $S$ induces a connected face cover\xspace $S^*$ in $D^*$, we have to prove that $S^*$ covers all vertices and the induced subgraph in the face-vertex incidence graph\xspace $H$ is connected. We proceed as in the other direction. $S$ covers all edges in $E(G)$, thus every edge $e \in E(D)$ is bounded by at least one face of $S^*$. Hence, every subdivision vertex in $V(D^*)$ is covered by a face of $S^*$. Furthermore, every vertex in $D^*$ is adjacent to a subdivision vertex, thus, also covered by a face in $S^*$.
Since $S^*$ is covering all vertices, we obtain from \cref{prp:subdiv} that $S^*$ is a connected face cover\xspace.
\qed
\end{proof}

\section{Feedback Vertex Set Approach}\label{sec:fvs}

A \emph{feedback vertex set} $S^\circ\subset V(G)$ of a graph $G$ is a vertex subset such that the induced subgraph $G[V(G)\setminus S^\circ]$ is acyclic.
We show here that finding a connected face cover\xspace $S$ of size $k$ for a plane biconnected graph $G$ is equivalent to finding a feedback vertex set $S^\circ \subset V(D)$ of size $k$ in the dual graph $D$ of $G$.
The \emph{weak dual}, i.e., the dual without a vertex for the outer face, of an outerplanar graph is a forest. Thus we must find the smallest number of splits in $G$ which transform $D$ into a forest.
In other words, we must must break all the cycles in $D$, and hence all of the vertices in the feedback vertex set $S^\circ$ of $D$ correspond to the faces of $G$ that should be merged together (see~\cref{fig:algsteps}).

\begin{figure}[t]
\centering
\begin{subfigure}[b]{.32\textwidth}
    \centering
    \includegraphics[page=2]{polyalg_steps_alt.pdf}
    \caption{primal (gray) and\\dual graph (black)}
    \label{fig:subfig3}
\end{subfigure}
\hfil
\begin{subfigure}[b]{.32\textwidth}
    \centering
    \includegraphics[page=3]{polyalg_steps_alt.pdf}
    \caption{spanning tree (black) of\\the connected face cover\xspace}
    \label{fig:subfig3}
\end{subfigure}
\hfil
\begin{subfigure}[b]{.32\textwidth}
    \centering
    \includegraphics[page=4]{polyalg_steps_alt.pdf}
    \caption{outerplanar graph after\\splitting along the tree}
    \label{fig:subfig3}
\end{subfigure}
	\caption
	The connected face cover\xspace\ (blue) is a feedback vertex set (red) in the dual.
	}
	\label{fig:algsteps}
\end{figure}

\begin{prp}\label{prp:fvsconnec}
	Let $H$ be the face-vertex incidence graph\xspace of a plane biconnected graph $G$ and let $S^\circ$ be a feedback vertex set in the dual $D$ of $G$.
	Then $S^\circ$ induces a connected face cover\xspace $S$ in $G$.
\end{prp}

\begin{proof}
We need to show that $S^\circ$ is a face cover and that it is connected.
First, assume there is a vertex $v \in V(G)$ of degree $\deg(v) = d$ that is not incident to a vertex in $S^\circ$, i.e., a face of $G$.
Since $G$ is biconnected, $v$ is incident to $d$ faces $f_1, \dots, f_d$, none of which is contained in $S^\circ$.
But then $D[V(D) \setminus S^\circ]$ has a cycle $(f_1, \dots, f_d)$, a contradiction.

Next, we define $\overline{S^\circ} = V(D)\setminus S^\circ$ as the complement of the feedback vertex set $S^\circ$ in $D$.
Assume that $H[V\cup S^\circ]$ has at least two separate connected components $C_1, C_2$.
Then there must exist a closed curve in the plane separating $C_1$ from $C_2$, which avoids the faces in $S^\circ$ and instead passes through a sequence $(f_1, \dots, f_\ell)$ of faces in $\overline{S^\circ}$, where each pair $(f_i, f_{i+1})$ for $i \in \{1, \dots, \ell-1\}$ as well as $(f_\ell, f_1)$ are adjacent in the dual $D$.
Again this implies that there is a cycle in $D[V(D) \setminus S^\circ]$, a contradiction.
Thus $S^\circ$ is a connected face cover\xspace. \qed
\end{proof}

\begin{theorem}\label{thm:split-fvs}
    A plane biconnected graph $G$ has outerplane splitting number $k$ if and only if its dual $D$ has a minimum feedback vertex set of size $k+1$.

\end{theorem}

\begin{proof}
    Let $S^\circ$ be a minimum feedback vertex set of the dual $D$ of $G$ with cardinality $|S^\circ| = k+1$ and let $H$ be the face-vertex incidence graph\xspace of $G$.
    We know from \cref{prp:fvsconnec} that $H' = H[V(G) \cup S^\circ]$ is connected and hence $S^\circ$ induces a connected face cover\xspace $S$ with $|S| = k+1$. Then by \cref{lem:cover-to-split} $G$ has $\ensuremath{\operatorname{osn}}(G) \le k$.

    Let conversely $\sigma$ be a sequence of $k$ vertex splits that turn $G$ into an outerplane graph $G'$ and let $F$ be the set of faces of $G$.
    By \cref{lem:split-to-cover} we obtain a connected face cover\xspace $S$ of size $k+1$ consisting of all faces that are merged by $\sigma$.
    The complement $\overline{S} = F\setminus S$ consists of all faces of $G$ that are not merged by the splits in $\sigma$ and thus are the remaining (inner) faces of the outerplane graph $G'$.
    Since $G'$ is outerplane and biconnected, $\overline{S}$ is the vertex set of the weak dual of $G'$, which must be a tree.
    Hence $S$ is a feedback vertex set in $D$ of size $k+1$ and the minimum feedback vertex set in $D$ has size at most $k+1$. \qed
\end{proof}

Since all faces in a maximal planar graph are triangles, the maximum vertex degree of its dual is 3.
Thus, we can apply the polynomial-time algorithm of Ueno et al.~\cite{UenoKG88} to this dual, which computes the minimum feedback vertex set in graphs of maximum degree 3 by reducing the instance to polynomial-solvable matroid parity problem instance, and obtain

\begin{corollary}
	We can solve \textsc{Outerplane Splitting Number}\xspace\ for maximal planar graphs in polynomial time.
\end{corollary}

Many other existing results for feedback vertex set  extend to \textsc{Outerplane Splitting Number}\xspace, e.g., it has a kernel of size $13k$~\cite{BonamyK16} and admits a \PTAS~\cite{DemaineH05}.

\section{Lower and Upper Bounds}\label{sec:bounds}

In this section we provide some upper and lower bounds on the outerplane splitting number in certain maximal planar graphs.

\subsection{Upper Bounds}

Based on the equivalence of \cref{thm:split-fvs} we obtain upper bounds on the outerplane splitting number from suitable upper bounds on the feedback vertex set problem, which has been studied for many graph classes, among them cubic graphs~\cite{BondyHS87}.
Liu and Zhao~\cite{LiuZ96} showed that cubic graphs $G=(V,E)$ of girth at least four (resp., three) have a minimum feedback vertex set of size at most $\frac{|V|}{3}$ (resp., $\frac{3|V|}{8}$).
Kelly and Liu~\cite{kl-msfvspgglf-17} showed that connected planar subcubic graphs of girth at least five have a minimum feedback vertex set of size at most $\frac{2|V|+2}{7}$.
Recall that the girth of a graph is the length of its shortest cycle.

\begin{proposition}\label{cor:upper}
	The outerplane splitting number of a maximal planar graph $G=(V,E)$ of minimum degree (i) 3, (ii) 4, and (iii) 5, respectively, is at most (i) $\frac{3 |V| - 10}{4}$, (ii) $\frac{2|V|-7}{3}$, and (iii) $\frac{4|V|-13}{7}$, respectively.
\end{proposition}

\begin{proof}
    Maximal planar graphs with $n = |V|$ vertices have $2n-4$ faces.
    So the corresponding dual graphs have $2n-4$ vertices.
    Moreover, since the degree of a vertex in $G$ corresponds to the length of a facial cycle in the dual, graphs with minimum vertex degree 3, 4, or 5 have duals with girth 3, 4, or 5, respectively.
    So if the minimum degree in $G$ is 3, we obtain an upper bound on the feedback vertex set of $(3n - 6)/4$; if the minimum degree is 4, the bound is $(2n-4)/3$; and if the minimum degree is 5, the bound is $(4n-6)/7$.
    The claim then follows from \cref{thm:split-fvs}. \qed
\end{proof}

\subsection{Lower Bounds}

We first provide a generic lower bound for the outerplane splitting number of maximal planar graphs.
Let $G$ be an $n$-vertex maximal planar graph with $2n-4$ faces.
Each face is a triangle incident to three vertices.
In a minimum-size connected face cover $S^*$, the first face covers three vertices.
Due to the connectivity requirement, all other faces can add at most two newly covered vertices.
Hence we need at least $\frac{n-1}{2}$ faces in any connected face cover.
By \cref{thm:split-cfc} this implies that $\ensuremath{\operatorname{osn}}(G) \ge \frac{n-3}{2}$.

\begin{proposition}
Any maximal planar graph $G$ has outerplane splitting number at least $\frac{|V(G)| - 3}{2}$.
\end{proposition}

Next, towards a better bound, we define a family of maximal planar graphs $T_d = (V_d, E_d)$ of girth 3 for $d \ge 0$ that have outerplane splitting number at least $\frac{2 |V_d| - 8}{3}$.
The family are the complete planar 3-trees of depth $d$, which are defined recursively as follows.
The graph $T_0$ is the 4-clique $K_4$.
To obtain $T_d$ from $T_{d-1}$ for $d \ge 1$ we subdivide each inner triangular face of $T_{d-1}$ into three triangles by inserting a new vertex and connecting it to the three vertices on the boundary of the face.

\begin{proposition}
  The complete planar 3-tree $T_d$ of depth $d$ has outerplane splitting number at least $\frac{2 |V_d| - 8}{3}$.
\end{proposition}

\begin{proof}
Each $T_d$ is a maximal planar graph with $n_d = 3 + \sum_{i=0}^{d} 3^i = \frac{3^{d+1} + 5}{2}$ vertices.
All $3^d$ leaf-level vertices added into the triangular faces of $T_{d-1}$ in the last step of the construction have degree 3 and are incident to three exclusive faces, i.e., there is no face that covers more than one of these leaf-level vertices.
This immediately implies that any face cover of $T_d$, connected or not, has size at least $3^d$.
From $n_d = \frac{3^{d+1} + 5}{2}$ we obtain $d = \log_3 \frac{2n_d-5}{3}$ and $3^d = \frac{2n_d-5}{3}$.
\cref{thm:split-cfc} then implies that $\ensuremath{\operatorname{osn}}(T_d) \ge \frac{2n_d-8}{3}$. \qed
\end{proof}

\section{Open Problems}
We have introduced the \textsc{Outerplane Splitting Number}\xspace\ problem and established its complexity for plane biconnected graphs.
The most important open question revolves around the embedding requirement.
Splitting operations can be defined more loosely and allow for any new embedding and neighborhood of the split vertices.
In general, it is also of interest to understand how the problem differs when the input graph does not have an embedding at all, as in the original splitting number problem.
Since \textsc{Outerplane Splitting Number}\xspace\ can be solved in polynomial time for maximal planar graphs but is hard for plane biconnected graphs, there is a complexity gap to be closed when faces of degree more than three are involved.
Vertex splitting in graph drawings has so far been studied to achieve planarity and outerplanarity. A natural extension is to study it for other graph classes or graph properties.

\bibliographystyle{splncs04}