∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>
11institutetext: The Abdus Salam International Centre for Theoretical Physics,
Strada Costiera 11, I - 34151, Trieste, Italy. 22institutetext: Kyungpook
National University, Center for High Energy Physics, Daegu, South Korea.
33institutetext: University of California at Irvine, Physics Department,
Irvine, CA 92697, USA.

# Down type iso-singlet quarks at the HL-LHC and FCC-hh

Arpon Paule1,addr1  Sezen Sekmene2,addr2  Gokhan Unele3,addr3

###### Abstract

We study the discovery potential of down type iso-singlet quarks, $D$,
predicted by the $E_{6}$ GUT model in the ${pp\rightarrow D\bar{D}\rightarrow
ZZd\bar{d}\rightarrow\ell^{+}\ell^{-}\ell^{+}\ell^{-}d\bar{d}}$ channel, at
the HL-LHC and FCC-hh colliders. The analysis is performed using a high level
analysis description language and its runtime interpreter. The study shows
that, using solely this channel, HL-LHC can discover $D$ quarks up to a mass
of 730 GeV whereas FCC-hh up to 2980 GeV with data collected in their complete
run periods.

###### Keywords:

Isosinglet quarks, FCC-hh, HL-LHC, CutLang

##  1 Introduction

The long-awaited discovery of the Higgs boson at the LHC experiments 20121  in
the year 2012 completed the experimental validation of the standard model
(SM). However, there are some well known issues that are not addressed by the
SM, such as the mass hierarchy problem, the unification of the fundamental
interactions, the origin of the baryon asymmetry of the Universe, and a
plausible explanation for dark matter. To address these issues, SM is proposed
to be extended into a more complete theory. In general, candidate extensions
predict the existence of new fundamental particles and interactions. The Large
Hadron Collider (LHC) experiments are conducting a great diversity of searches
for discovering these new particles and interactions. Results of all these
searches so far have been found to be consistent with the SM predictions. The
forthcoming High Luminosity Large Hadron Collider (HL-LHC) and Future Circular
Collider (FCC) machines, with their higher luminosity, energy and better
detector acceptance and efficiency, will increase the sensitivity of these
searches, and expand them to more difficult scenarios, enabling access to
higher particle masses and lower effective cross sections.

One class of candidate extensions to the SM consists of Grand Unified Theories
(GUTs) based on a gauge group larger than that of the SM. The GUT models merge
strong and electroweak interactions in a single gauge group, thereby allowing
a solution to at least two of the above mentioned problems, namely, the
complete unification of the fundamental interactions (except gravity) and the
baryon asymmetry of the observed Universe. Specifically, when unifying gravity
with other interactions both within the contexts of the superstring and
supergravity theories, the exceptional Lie group $E_{6}$ has been shown to be
the gauge symmetry group which can be compactified from 10 (or 11) dimensions
down to the $3+1$ that we observe HEWETT1989193 .

The GUT model using the Exceptional Lie Group $E_{6}$ as the gauge symmetry
group is referred to as the $E_{6}$ model. It predicts the existence of iso-
singlet quarks (in literature, denoted by $D,S,$ and $B$) having charge
$Q=-1/3$. The discovery potential of the ATLAS experiment for the down type
iso-singlet quark $D$ of the first SM family has been previously investigated
in Mehdiyev2007 ; Mehdiyev2008 . The discovery reach for $D$ quarks were
estimated at a phenomenology study before the LHC data taking to be 950 GeV
for 100 fb-1 integrated luminosity using the combination of all $D$ decay
channels Mehdiyev2008 .

Dedicated searches for down-type iso-singlet quarks predicted by the model
described in this paper in the LHC data are currently ongoing in the ATLAS
experiment. In the meanwhile, the closest estimates of sensitivity come from
searches for vector-like quarks (VLQs), which have similar production
mechanisms. However almost all existing VLQ searches are exclusively designed
to target third generation vector-like partners $B$ and $T$ of the bottom and
top quarks. The most stringent limits to date come from an ATLAS combination
of 7 VLQ searches performed with 13 TeV data, looking at different final
states Aaboud:2018xuw ; Aaboud:2017zfn ; Aaboud:2018uek ; Aaboud:2017qpr ;
Aaboud:2018saj ; Aaboud:2018xpj ; Aaboud:2018wxv , which excluded $T$ ($B$)
masses below 1.31 (1.03) TeV for any combination of decays into SM particles
Aaboud:2018pii . Several CMS studies also searched for third generation VLQs.
One search in the single lepton channel excluded $T$ masses less than 1.295
TeV in exclusive decays to $tW$ Sirunyan:2017pks . A different fully hadronic
search excluded $T$ and $B$ quark masses between 0.74-1.37 TeV
Sirunyan:2019sza , while a leptonic search excluded $T$ quarks with masses
below 1.14 to 1.30 TeV and $B$ quarks with masses below 0.91 to 1.24 TeV
Sirunyan:2018omb  for various branching fraction combinations. However those
limits do not directly apply to the down-type isosinglet quark scenario
studied here, as the searches mainly focus on third generation final states
that contain $b$ quarks. For the specific case of light-flavor VLQ, a CMS
search with at least one lepton excluded pair-produced VLQs below masses 845
and 685 GeV for branching ratios $B(W)=1$ and $B(W)=0.5$, $B(Z)=B(H)=0.25$,
respectively Sirunyan:2017lzl .

In this study, we investigate the possibility of observing the pair production
of first generation down type iso-singlet quarks $D$, in the decay channel
$D\rightarrow Zd\rightarrow\ell^{+}\ell^{-}d$ (where $\ell=e,\mu$) using the 4
leptons plus 2 jets final state at the HL-LHC and the proton-proton scenario
for the FCC. Due to its low effective cross section, this process could not be
observed at the current LHC conditions. With their higher luminosity, energy
and detector acceptances, HL-LHC and FCC are expected to significantly improve
sensitivity in this channel. Despite its low effective cross section,
exploring this channel is critical, as it provide the most precise
reconstruction of the $D$ quark mass in case of discovery.

Additionally, this work aims to test the feasibility of a new and practical
analysis writing approach for high energy physics. The search method in this
study is implemented and performed using an analysis description language and
its runtime interpreter CutLang, which allows quick analysis prototyping and
histogramming Unel:2019reo ; Sekmen:2018ehb .

The paper starts by introducing the down-type iso-singlet quark model in
Section 2 followed by a description of the HL-LHC and FCC colliders and
relevant experimental conditions in Section 3, and the analysis description
language and runtime interpreter CutLang in Section 4. Detailed explanation of
the search for $D$ quarks and the search results are presented in Section 5
followed by the conclusions in Section 6.

##  2 Down-type iso-singlet quark model

If the group structure of the SM, $SU_{C}(3)\times SU_{W}(2)\times U_{Y}(1)$,
originates from the breaking of the $E_{6}$ group at the GUT scale, then the
quark sector of the SM is extended by the addition of an iso-singlet quark per
family as:

$\displaystyle\begin{pmatrix}u_{L}\\\ d_{L}\\\ \end{pmatrix}$ , $\displaystyle
u_{R},d_{R},D_{L},D_{R};$ (1) $\displaystyle\begin{pmatrix}c_{L}\\\ s_{L}\\\
\end{pmatrix}$ , $\displaystyle c_{R},s_{R},S_{L},S_{R};$
$\displaystyle\begin{pmatrix}t_{L}\\\ b_{L}\\\ \end{pmatrix}$ , $\displaystyle
t_{R},b_{R},B_{L},B_{R}.$

In the considered model, the $S$ and $B$ quarks are assumed to be heavy and
decoupled from the spectrum, leaving the $D$ quark as the only one accessible
for searches at the present and near future colliders. A second assumption,
following from the general behavior of CKM (Cabibbo, Kobayasi, Maskawa), is
that mixing inside a given family is stronger compared to mixing between
different families. Therefore, we only consider the Lagrangian relevant for
the weak interaction of $d$ and $D$ quarks as given in _akir_1997 :

$\begin{split}L_{D}&=\frac{\sqrt{4\pi\alpha_{em}}}{2\sqrt{2}\sin{\theta_{W}}}\big{[}\bar{u}^{\theta}\gamma_{\alpha}(1-\gamma_{5})d\cos{\phi}\\\
&\quad+\bar{u}^{\theta}\gamma_{\alpha}(1-\gamma_{5})D\sin{\phi}\big{]}W^{\alpha}\\\
&\quad-\frac{\sqrt{4\pi\alpha_{em}}}{4\sin{\theta_{W}}}\left[\frac{\sin{\phi}\cos{\phi}}{\cos{\theta_{W}}}\bar{d}\gamma_{\alpha}(1-\gamma_{5})D\right]Z^{\alpha}\\\
&\quad-\frac{\sqrt{4\pi\alpha_{em}}}{4\cos{\theta_{W}}\sin{\theta_{W}}}\times\\\
&\qquad\quad\big{[}\bar{D}\gamma_{\alpha}(4\sin^{2}{\theta_{W}}-3\sin^{2}{\phi}(1-\gamma_{5}))D\\\
&\qquad\quad+\bar{d}\gamma_{\alpha}(4\sin^{2}{\theta_{W}}-3\sin^{2}{\phi}(1-\gamma_{5}))d\big{]}Z^{\alpha}\\\
&+h.c.\quad,\end{split}$ (2)

where the superscript $\theta$ represents the usual CKM mixings taken to be in
the up sector for simplicity of calculation, $\theta_{W}$ is the weak mixing
angle and $\phi$ is the mixing angle between the $d$ and $D$ quarks, which is
responsible for the decay of the $D$ quark. The limits on $\phi$ can be
obtained from the current precision measurements for the $3\times 3$ CKM
matrix elements, assuming that its $3\times 4$ extension has the sum of the
squares of the elements of a row equal to 1.

The evaluation of the presently measured values and their errors yield
${\lvert\sin{\phi}\rvert}\leq 0.035$ (0.043) allowing a 1(2) sigma variation
on the first row elements PhysRevD.98.030001 . The cross section calculation
results for FCC are essentially insensitive to $\sin(\phi)$, since the studied
pair production proceeds mostly via gluon exchange. However at HL-LHC,
especially for large values of $D$ quark mass, the production is mostly via
the $q\bar{q}$ channel which has a slight $\sin(\phi)$ dependence for the
cross section due to the $t$ channel sub-process propagating via $W$ boson as
shown in 1, sub-figure (d). Since this sub-process contributes with an
opposite sign, reducing the mixing angle effectively increases the $D\bar{D}$
production cross section.

The branching fractions for the three possible $D$ decay modes, ${D\rightarrow
Wu}$, ${D\rightarrow Zd}$ and ${D\rightarrow hd}$ are about $50\%$, $25\%$ and
$25\%$ respectively for masses above $\sim 800$ GeV Sultansoy:2006cw . In this
study, we consider the pair production of $D$ quarks and their subsequent
decay in the $D\rightarrow Zd$ channel to explore the discovery prospects of
two possible future collider scenarios.

##  3 Considered collider scenarios

###  3.1 High-Luminosity LHC

The LHC reached its design value of peak luminosity $10^{34}~{}cm^{-2}s^{-1}$
in June, 2016. The High-Luminosity Large Hadron Collider (HL-LHC) project aims
to improve the performance of the LHC in order to increase the potential for
discoveries after 2027 Apollinari:2120673 ; doi:10.1142/9581 ;
ApollinariG.:2017ojx . To implement this, HL-LHC will have several cutting-
edge technologies, such as, 11–12 T superconducting magnets; very compact with
ultra-precise phase control superconducting cavities for beam rotation; new
technology for beam collimation; and long high-power superconducting links
with zero energy dissipation. HL-LHC is expected to reach the peak luminosity
of $5\times 10^{34}$ $cm^{-2}s^{-1}$, allowing an integrated luminosity of
$250fb^{-1}$ per year. Therefore, it gives an integrated luminosity of
$3000~{}fb^{-1}$ in the operation period of about a dozen years after the
upgrade. This integrated luminosity corresponds to ten times the amount LHC is
expected to collect after 12 years of operation.

To meet the challenges brought by this higher luminosity at the HL-LHC, such
as higher radiation dose, higher particle rate, higher pileup, and higher
event rate, etc, the ATLAS and CMS detectors will undergo an extensive upgrade
(i.e. the “Phase 2” upgrade). The ATLAS inner tracker (ITk) is being
completely rebuilt for Phase 2, as a result of which, the pseudorapidity
coverage will extend up to $|\eta|=4$. Moreover, new front-end electronics and
a new readout system in the calorimeters will allow triggering higher
resolution objects at the lowest trigger level at an increased rate, and lead
to improved reconstruction. In addition, new inner barrel chambers will be
installed in the muon detector system for increased coverage. The CMS detector
will similarly undergo major upgrades which include a replacement of the
silicon strip and pixel components in the tracking detector increasing the
coverage up to $|\eta|=4$. The hadronic calorimeter will be read out by
silicon photomultipliers. The endcap electromagnetic and hadron calorimeters
will be replaced with a new combined sampling calorimeter that will provide
highly-segmented spatial information in both the transverse and longitudinal
directions, as well as high-precision timing information. The muon system will
be extended with new chambers in the forward region, bringing the coverage up
to $|\eta|=2.8$. Additionally, both ATLAS and CMS envisage adding timing
detectors to provide the capability of adding timing information to
reconstruction CERN-LHCC-2015-020 ; Contardo:2020886 .

###  3.2 Future Circular Collider

The Future Circular Collider (FCC) was launched as a world-wide international
collaboration hosted at CERN in response to the 2013 Update of the European
Strategy for Particle Physics (EPPSU) Abada2019 ; Benedikt:2653674 . In the
2020 Update of EPPSU, it has been proposed to investigate the technical and
financial feasibility of FCC European:2720129 . FCC scenarios are studied for
three different types of particle collisions, namely hadron (proton-proton and
heavy ion), electron-positron and proton-electron collisions. The proposed
energy frontier proton-proton collider, FCC-hh, which is considered in this
study, is designed to provide proton–proton collisions with a centre-of-mass
energy of $100~{}TeV$ and an integrated luminosity of $20~{}ab^{-1}$ for 25
years of operation. The FCC-hh collider layout has two high luminosity
interaction points for general purpose detectors. The factor 7 increase in
energy over the present LHC requires a vast modification compared to the
designs of current general purpose LHC detectors. The detectors for 100 TeV
should be able to measure multi-TeV jets, leptons and photons from heavy
resonances with masses up to 50 TeV, while at the same time measuring the
known SM processes with high precision, and still being sensitive to a broad
range of BSM signatures with moderate momentum. In addition, future detectors
will need to operate at $~{}1000$ pileup events per bunch-crossing. The
detector acceptance is targeted to increase up to $|\eta|=4.4$ in order to
improve sensitivity to vector boson fusion processes.

##  4 CutLang analysis description language and runtime interpreter

As mentioned earlier, one goal of this study is to test the feasibility of the
new “analysis description language” approach in analysis writing and running
in phenomenological studies. An analysis description language is a domain-
specific, declarative language designed to express the physics contents of an
analysis in a standard and unambiguous way. In this approach, the description
of the analysis components is decoupled from the software framework that run
the analysis.

This study uses the language ADL Brooijmans:2016vro ; Brooijmans:2018xbu ;
Brooijmans:2020yij , which consists of a plain text file containing blocks
with a keyword-value structure. The blocks make clear the separation of
analysis components such as object definitions, variable definitions, and
event selections while the keywords specify analysis concepts and operations.
The syntax includes mathematical and logical operations, comparison and
optimization operators, reducers, four-vector algebra and common HEP-specific
functions (e.g. $\delta\phi$, $\delta R$, etc.). ADL files can refer to self-
contained functions encapsulating variables with complex algorithms (e.g.
$M_{T2}$, aplanarity, etc.) or non-analytic variables (e.g. efficiency tables,
machine learning discriminators, etc.).

ADL can be used for performing an analysis by any framework capable of
interpreting and running it. Here, we use CutLang Sekmen:2018ehb ;
Unel:2019reo , a runtime interpreter, which is able to operate directly on
events without the need for compilation. CutLang is written in C++ and is
based on ROOT ROOT  classes for Lorentz vector operations and histogramming.
It uses automatically generated dictionaries and grammar rules based on unix
tools Lex and Yacc lexandyacc  . The typical output of an analysis in CutLang
is a file containing surviving events and histograms which can be used for
statistical analysis.

Not having the necessity to write or compile code, combined with the simple,
human-readable nature of ADL syntax makes it a very practical construct for
quickly performing phenomenological analyses such as the one in this study.

##  5 Search for down-type iso-singlet quarks

###  5.1 Signal and background processes

The main tree level Feynman diagrams for the pair production of $D$ quarks at
hadron colliders are presented in Figure 1. The model Lagrangian in equation
(2) was implemented into the tree level event generator, CompHEP BOOS2004250 ;
Pukhov:1999gg . The resulting pair production cross sections at generator
level for HL-LHC and FCC-hh for the $gg$ and $q\bar{q}$ channels and their sum
are shown in Figure 2 as a function of $D$ quark mass. The pair production
cross section is somewhat smaller than the single production, for example for
a $D$ quark of 1 TeV the former is 38.6 fb whereas it is 94.5 fb for the
latter. However as the single production results depend heavily on the mixing
angle and the SM background is especially large due to QCD jets this paper
focuses on pair production.

(a) gluons, s channel

(b) gluons, t channel

(c) up quarks, s channel

(d) up quarks, t channel

Figure 1: Tree level Feynman Diagrams for the process $pp\rightarrow D\bar{D}$
Figure 2: $pp\rightarrow D\bar{D}$, $q\bar{q}\rightarrow D\bar{D}$ and
$gg\rightarrow D\bar{D}$ cross sections vs $D$ quark mass for HL-LHC and FCC-
hh energies, calculated using CompHep. The $d-D$ mixing angle is taken as
$\sin{\phi}=0.035$

The $E_{6}$ GUT model does not predict the masses of the iso-singlet quarks.
Therefore, this study scans some plausible values for the $D$ quark mass (up
to $2500$ GeV) to estimate the experimental reach at both HL-LHC and FCC-hh
machines. The iso-singlet quarks are expected to immediately decay into SM
particles due to their large masses. In this analysis, we have considered the
decay process $D\bar{D}\rightarrow ZZd\bar{d}$, with subsequent leptonic
decays of both $Z$ bosons, $Z\rightarrow\ell^{+}\ell^{-}$.

The main SM background to the signal process is $pp\rightarrow ZZjj$
production, with subsequent leptonic decays of both Z bosons. The SM cross-
section of $pp\rightarrow ZZjj$ is calculated using MadGraph5_aMC@NLO
Alwall2014  considering up to 4 QED and QCD interaction vertices and found to
be $2.918$ pb and $68.04$ pb for HL-LHC and FCC-hh, respectively.

Processes with Higgs decaying to two $Z$ bosons also provide final states
resembling that of the signal, however they are not considered as significant
backgrounds in this study due to relatively low effective cross sections as
well as one of the $Z$ bosons from the Higgs boson decay being virtual. At 14
TeV, the Higgs production cross sections are estimated as 54.6 pb from gluon
fusion, 4.3 pb from VBF, 1.5 pb from $WH$, 0.98 pb from $ZH$ and 0.55 pb from
$bbH$ production channels. To obtain an estimate for $ZZjj$ final states,
these numbers are multiplied by the $h\rightarrow ZZ$ branching fraction and
the hadronic branching fraction of $W$ and $Z$ bosons. Moreover, the gluon
fusion cross section is corrected to account for multi-jet events
Greiner:2124367 . Extrapolating linearly from 8 and 13 TeV results, the $h+2j$
cross section from gluon fusion is estimated as 5.4 pb at NLO level. Folding
in the appropriate branching fractions, the total effective cross section for
the Higgs-related backgrounds becomes $\sim$0.31 pb, which is a small fraction
of the direct $ZZjj$ production cross section. The approximate estimate of
these processes for 100 TeV is $\sim$5 pb, which is similarly small compared
to the SM $ZZjj$ cross section. Moreover, one of the $Z$ bosons originating
from the Higgs decays would be virtual. Therefore the majority of such events
would be rejected by the requirement of two reconstructed $Z$ bosons having an
invariant mass of 91.2 GeV in our analysis.

The $E_{6}$ model signal events with $D$ quarks decaying to SM particles and
SM background events were generated using CompHEP and MadGraph5_aMC@NLO
respectively. The CompHEP setup was adjusted to impose a generator level
requirement of 10 GeV on the transverse momenta of the SM d-quarks originating
from the $D\rightarrow Zd$. The NNPDF 3.1 parton distribution function set
Ball:2017nwa , which is the most up-to-date set available has been used both
for 14 and 100 TeV. Further decays and showering and hadronization processes
were simulated using Pythia6 Sjostrand:2006za . Pythia was set up to only
allow electron and muon decays of the $Z$ bosons. Subsequently, the detector
effects were modelled with the fast detector simulation program Delphes
deFavereau2014  using the configurations RefHL-LHCcard  and RefFCCcard  for
generic HL-LHC and FCC-hh detectors.

###  5.2 Object and event reconstruction and selection

The complete object and event reconstruction and selection algorithm for the
analysis is given in ADL format in Table 1. This is, in fact, the exact ADL
code run in CutLang to produce the results presented in this paper.

The analysis is performed in the $4\ell+2j$ channel, and thus uses leptons and
jets. Both for HL-LHC and FCC-hh cases, leptons considered are electrons and
muons, which are both required to have transverse momentum $p_{T}>20$ GeV and
pseudorapidity $|\eta|<4$. Electrons (muons) are required to have an isolation
of 0.1 (0.2) within a cone of $dR<0.3$. Jets are reconstructed with the
anti-$k_{T}$ algorithm with a radius of $R=0.5$, and are required to have
$p_{T}>50$ GeV (which is higher than generator level requirement) and
$|\eta|<4$. Increased pseudorapidity acceptance at the HL-LHC and FCC-hh
detectors compared to LHC will provide an increased sensitivity for the
analysis. Events are required to have at least 4 leptons and at least 2 jets
as defined above.

####  5.2.1 Leptonic Z boson reconstruction

The two $Z$ boson candidates from the $D$ decay are reconstructed from the
selected leptons. For an efficient Z boson reconstruction, we consider the
following criteria:

  1. 1.

Mass of the reconstructed Z boson candidate should be as close as possible to
91.2 GeV,

  2. 2.

the Z boson candidate should be flavour and charge neutral (i.e, reconstructed
from a $e^{+}e^{-}$ or a $\mu^{+}\mu^{-}$ pair)

In this analysis, we are focused on final states with Z bosons with moderate
momentum, which decay to non-collimated leptons that can be independently
reconstructed. However, especially at the FCC-hh energies, higher mass $D$
quarks yield a Z boson $p_{T}$ spectra with a higher component of boosted Z
bosons that would decay to collimated lepton pairs. Such collimated lepton
pairs would partially fail to be identified as two individual leptons due to
the lepton isolation requirement and be counted as a single lepton, resulting
in the event failing the 4 lepton criteria. A more effective treatment of the
boosted final states would require Z boson reconstruction via explicit tagging
of the boosted Z boson via collimated lepton jets. These boosted channels can
be added when collimated lepton jet tagging performance or simulation for the
FCC-hh conditions become available, and they would increase the analysis
sensitivity.

For the resolved final state, leptons are paired to reconstruct both $Z$
bosons simultaneously in the $\chi^{2}$ expression below, which both selects
the dilepton combinations with masses as close as possible to the measured $Z$
mass of 91.2 GeV and ensures the same flavor requirement on dileptons in a
candidate:

$\displaystyle\chi^{2}_{ZZ}$ $\displaystyle\equiv$
$\displaystyle(m_{Z1}-91.2)^{2}+(m_{Z2}-91.2)^{2}$ $\displaystyle+$
$\displaystyle\left(999\times
PdgID\left[Z_{1}\right]\right)^{2}+\left(999\times
PdgID\left[Z_{2}\right]\right)^{2}.$

More information on technical implementation of Z reconstruction and the
$\chi^{2}_{ZZ}$ in ADL and CutLang is given in Appendix A. The reconstructed
$Z$ candidates are additionally required to have a total electric charge of 0.
Mass distributions of both $Z$ candidates reconstructed from $e^{+}e^{-}$ and
$\mu^{+}\mu^{-}$ pairs are shown in Figure 3 for different $D$ quark masses
for HL-LHC and FCC-hh.

Figure 3: Invariant mass distribution for both reconstructed $Z$ boson
candidates for HL-LHC (top) and FCC-hh (bottom) conditions. Candidates are
reconstructed from both $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$ pairs.

####  5.2.2 $D$ quark reconstruction

Each $D$ quark candidate ($D_{1}$ and $D_{2}$) is reconstructed from a $Z$
boson candidate and a jet. Once again, the reconstruction is based on a
$\chi^{2}$ optimization which takes into account the following conditions:

  1. 1.

$D$ quark mass is presumed unknown. However, masses of the two reconstructed
$D$ quark candidates should be as close as possible to each other. We express
this condition as:

$\chi^{2}_{m_{D}}\equiv((m_{D_{1}}-m_{D_{2}})/m_{D})^{2}\;,$ (4)

where $m_{D_{1}}$ and $m_{D_{2}}$ are the invariant masses of the two $D$
quark candidates and $m_{D}=(m_{D_{1}}+m_{D_{2}})/2$.

  2. 2.

Transverse momentum of the jets directly originating from the $D$ quark decay
is expected to be high. To ensure selecting jets with high momentum, we use
the Heavyside step function with a weight factor:

$\displaystyle\chi^{2}_{p_{T,j}}$ $\displaystyle\equiv$ $\displaystyle
H(p_{T,j}^{cut}-p_{T,j_{1}})\times((p_{T,j}^{cut}/p_{T,j_{1}})-1.0)$ (5)
$\displaystyle+$ $\displaystyle
H(p_{T,j}^{cut}-p_{T,j_{2}})\times((p_{T,j}^{cut}/p_{T,j_{2}})-1.0),$

where $p_{T,j_{1}}$ and $p_{T,j_{2}}$ are the transverse momenta of the jets
and $p_{T,j}^{cut}$ is the selection threshold to be applied to the jet
transverse momenta. To determine the optimal value for this threshold which
would obtain the best signal-background separation, we show the $p_{T}$
distributions of the candidate jets in Figure 4 for signals with different
$m_{D}$ and the background at HL-LHC and FCC-hh. The jets are selected by
minimizing the condition defined in Eq. 4. Due to its much higher center-of-
mass energy, FCC-hh yields a much harder jet $p_{T}$ spectrum. Based on these
distributions, we select $p_{T,j}^{cut}=300$ and $500$ GeV as thresholds for
HL-LHC and FCC-hh respectively.

Figure 4:  Transverse momentum distribution for both jets used in $D$ quark
reconstruction for HL-LHC (top) and FCC-hh (bottom) conditions. The jets are
selected by minimizing the condition defined in Eq. 4.

  3. 3.

Angular separation between the two $D$ quarks,

$dR_{DD}=\sqrt{(\eta_{D_{1}}-\eta_{D_{2}})^{2}+(\phi_{D_{1}}-\phi_{D_{2}})^{2}}\;,$
(6)

should reflect that the $D$ quarks are centrally produced, with negligible
Lorentz boost. The most characteristic configuration would correspond to $D$
quarks having $|\eta|\simeq 0$ and being back-to-back on the transverse plane,
which gives $\delta\phi\simeq\pi$, where $\delta\phi$ represents the $\phi$
difference of the two particles. As a result, $dR$ is expected to be dominated
by $\delta\phi$ and peak around 3.14. This can be seen in Figure 5, which
shows the $dR$ distributions for signals and the background for HL-LHC and
FCC-hh, after applying a minimization based on Eq 4. Both signals and the
background peak around 3.14, but the backgrounds display a wider distribution.
Based on this information, we define a variable that can be minimized to zero:

$\displaystyle\chi^{2}_{dR_{DD}}$ $\displaystyle\equiv$
$\displaystyle(dR_{DD}/3.14-1.0)^{2}.$ (7)

Figure 5: Distribution of angular distance between the two reconstructed $D$
quark candidates $D_{1}$ and $D_{2}$ for HL-LHC (top) and FCC-hh (bottom)
conditions. The jets are selected by minimizing the condition defined in Eq.
4.

We then combine the three conditions in Eqs 4, 5 and 7 to obtain a $\chi^{2}$
and select the $D$ candidates by running a minimization based on the sum:

$\chi^{2}_{DD}\equiv\chi^{2}_{m_{D}}+\chi^{2}_{p_{T,j}}+\chi^{2}_{dR_{DD}}\simeq
0\quad.$ (8)

Here, we tried different relative weighting of $\chi^{2}_{m_{D}}$,
$\chi^{2}_{p_{T,j}}$ and $\chi^{2}_{dR_{DD}}$, but the above choice gives the
optimal result.

####  5.2.3 Final selection on $\chi_{DD}^{2}$

Figure 6 shows the distribution of $\chi^{2}_{DD}$ values obtained after
minimization for HL-LHC (top) and FCC-hh (bottom) conditions for signals with
different $m_{D}$ and background. As expected, the signals exhibit a
distribution much closer to zero compared to the background. A selection of
$\chi_{DD}^{2}<0.5$ was applied to further reduce the SM contamination. The
threshold value was chosen to ensure a high signal significance.

Figure 6: Distribution of $\chi^{2}_{DD}$ values obtained after minimization
for HL-LHC (top) and FCC-hh (bottom) conditions. Table 1: Analysis description
using the ADL/CutLang syntax. This description can be directly processed with
CutLang over events.

⬇

1## Object definitions

2object goodJet

3 take JET

4 select Pt(JET) > 50

5 select abs(Eta(JET) < 4

6object goodEle

7 take ELE

8 select Pt(ELE) > 20

9 select abs(Eta(ELE)) < 4

10object goodMuo

11 take MUO

12 select Pt(MUO) > 20

13 select abs(Eta(MUO)) < 4

14object goodLep : Union (goodEle,goodMuo)

15## Reconstructed particles

16define Zreco1 = goodLep[-1] goodLep[-1]

17define Zreco2 = goodLep[-3] goodLep[-3]

18define dj1 = goodJet[-2]

19define dj2 = goodJet[-4]

20define Dreco1 = Zreco1 dj1

21define Dreco2 = Zreco2 dj2

22## Event variables

23define PTj1 = Pt(dj1)

24define PTj2 = Pt(dj2)

25define mZ1 = m(Zreco1)

26define mZ2 = m(Zreco2)

27define mD1 = m(Dreco1)

28define mD2 = m(Dreco2)

29define mD = ( mD1 + mD2 ) / 2

30define dRDD = dR( Dreco1 , Dreco2 )

31# define PTjcut = 300 # for HL-LHC

32define PTjcut = 500 # for FCC-hh

33define chi2DDcut = 0.5

34## Chi2 variable definitions

35define chi2ZZ = (mZ1 - 91.2)^2 + (mZ2 - 91.2)^2 + (999*pdgID(Zreco1))^2 +
(999*pdgID(Zreco2))^2

36define chi2mD = ((mD1 - mD2)/mD)^2

37define chi2PTj = Hstep(PTjcut -PTj1)*(PTjcut/PTj1 - 1.0) + Hstep(PTjcut -
PTj2)*(PTjcut/PTj2 - 1.0)

38define chidRDD = (dRDD/3.14 - 1.0)^2

39define chi2DD = chimD + chiPTj + chidRDD

40## Event selection

41region DDselection

42 select ALL

43 select Size(goodEle) >= 0

44 select Size(goodMuo) >= 0

45 select Size(goodLep) >= 4

46 select chi2ZZ ~= 0

47 ## chi^2 optimization for Z reconstruction

48 select q(Zreco1) == 0 # Z is neutral

49 select q(Zreco2) == 0 # Z is neutral

50 histo hmZ1, "Z candidate1 mass (GeV)", 320, 0.0, 3200.0, mZ1

51 histo hmZ2, "Z candidate2 mass (GeV)", 320, 0.0, 3200.0, mZ2 # mZ1 & mZ2
histogram plotting (fig. 3)

52 select Size ( goodJet ) >= 2

53 select chi2DD ~= 0

54 ## chi^2 optimization for D reconstruction

55 histo hchi2DD, "chi2DD ", 200, 0.0, 10.0, chi2DD # chi2DD histogram
plotting (fig. 6)

56 select chi2DD < chi2DDcut

57 histo hmD, "D candidate mass (GeV)", 320, 0.0, 3200.0, mD # mD histogram
plotting (fig. 7 & 8)

Table 2: Percentage selection efficiencies for various signals and background for the HL-LHC selection. Cumulative Selection criteria |  | Selected events (% of Total)
---
Background | Signal
600 GeV | 800 GeV | 1000 GeV
`ALL` | 100 | 100 | 100 | 100
`Size(goodLep) >= 4` | 20.1 | 26.7 | 28.1 | 30.0
`chi2ZZ ~= 0` | 20.1 | 26.7 | 28.1 | 30.0
`{Zreco1}q == 0` | 20.1 | 26.7 | 28.1 | 30.0
`{Zreco2}q == 0` | 20.1 | 26.7 | 28.1 | 30.0
`Size(goodJet) >= 2` | 5.03 | 25.0 | 26.9 | 29.3
`chi2DD ~= 0` | 5.03 | 25.0 | 26.9 | 29.3
`chi2DD < chi2DDcut` | 0.187 | 8.18 | 15.2 | 22.5
Table 3: Percentage selection efficiencies for various signals and background for the FCC-hh selection. Cumulative Selection criteria |  | Selected events (% of Total)
---
Background | Signal
800 GeV | 1600 GeV | 2500 GeV
`ALL` | 100 | 100 | 100 | 100
`Size(goodLep) >= 4` | 28.3 | 42.2 | 48.6 | 51.6
`chi2ZZ ~= 0` | 28.3 | 42.2 | 48.6 | 51.6
`{Zreco1}q == 0` | 28.3 | 42.2 | 48.6 | 51.6
`{Zreco2}q == 0` | 28.3 | 42.2 | 48.6 | 51.6
`Size ( goodJet ) >= 2` | 13.2 | 41.0 | 48.1 | 51.3
`chi2DD ~= 0` | 13.2 | 41.0 | 48.1 | 51.3
`chi2DD < chi2DDcut` | 0.302 | 11.2 | 36.2 | 46.0

###  5.3 Results

The percentage selection efficiencies for signal and background events for the
event selection criteria described above are given in Tables 2 and 3 for HL-
LHC and FCC-hh. Overall signal selection efficiency is seen to increase as $D$
mass increases.

The distribution of the average reconstructed $D$ quark invariant mass
$(m_{D_{1}}+m_{D_{2}})/2$ in the signal and background events that remain
after selection are shown in Figures 7 and 8 for different generated $D$ quark
masses for HL-LHC and FCC-hh, respectively. Signal events are seen to peak
visibly over the falling background distributions. In order to reduce the
statistical fluctuations due to limited amount of statistics, the signal and
background distributions can be modelled with a Gaussian function and a
Crystal Ball function, respectively. The signal and background yields are
obtained from the total events distribution, by fitting it to the sum of these
two functions. The initial fit parameters for the Crystal Ball and Gaussian
functions were determined by performing independent fits to the signal and
background distributions. The resulting fits are also shown in the same
figures.

Figure 7: Distribution of average reconstructed $D$ quark invariant mass
$(m_{D_{1}}+m_{D_{2}})/2$ for HL-LHC conditions for background and signals
with $m_{D}=600$ GeV (top), $800$ GeV (middle) and $1000$ GeV (bottom).
Results of the fit to the sum of a Gaussian and Crystal Ball functions are
also shown.

Figure 8: Distribution of average reconstructed $D$ quark invariant mass
$(m_{D_{1}}+m_{D_{2}})/2$ for FCC-hh conditions for background and signals
with $m_{D}=800$ GeV (top), $1600$ GeV (middle) and $2500$ GeV (bottom).
Results of the fit to the sum of a Gaussian and Crystal Ball functions are
also shown.

The fit results are then used for estimating the final signal and background
yields denoted as $S$ and $B$. These are obtained by integrating the fitted
Gaussian and Crystal Ball functions in a range defined by two standard
deviations mass window around the Gaussian mean. The obtained values for each
$D$ quark mass are then used for calculating the signal significance
$\sigma_{DD}$ defined as:

$\displaystyle\sigma_{DD}$
$\displaystyle\equiv\sqrt{2\times\left[\left(S+B\right)\ln{\left(1+\frac{S}{B}\right)}-S\right]}\;.$
(9)

The yields $S$ and $B$ along with the significance obtained for each simulated
mass point are shown in Tables 4 and 5 for HL-LHC and FCC-hh, respectively.
Signal significance values are also shown in Figure 9, plotted against the $D$
quark mass. A linear function is fitted to the plot to estimate the dependence
of significance on $D$ quark mass. The $D$ quark mass values, for which it
would be possible to make an observation ($3\sigma$) or a discovery
($5\sigma$), are then calculated from the linear function obtained from the
fit, and are shown in Table 6 for HL-LHC and FCC-hh. Finally, the integrated
luminosities required for $3\sigma$ observation and $5\sigma$ discovery at HL-
LHC and FCC-hh are plotted versus $D$ quark mass in Figure 10.

Figure 9: Signal significance as a function of $D$ quark mass for HL-LHC and
FCC-hh.

Figure 10: The integrated luminosity needed for $3\sigma$ observation and $5\sigma$ discovery as a function of $D$ quark mass for HL-LHC (top) and FCC-hh (bottom). Table 4: Signal and background yields and significance for different $D$ quark masses at HL-LHC. D quark mass | S | B | $\sigma$
---|---|---|---
600 GeV | 46 | 19 | 8.2
800 GeV | 15 | 10 | 3.9
1000 GeV | 6 | 7 | 1.9
Table 5: Signal and background yields and significance for different $D$ quark masses at FCC-hh. D quark mass | S | B | $\sigma$
---|---|---|---
800 GeV | 21359 | 2690 | 250
1600 GeV | 2123 | 740 | 59
2500 GeV | 241 | 318 | 12
Table 6: Upper limit on $D$ quark masses for $3\sigma$ observation and $5\sigma$ discovery for HL-LHC and FCC-hh. Experiment |  $3\sigma$ observation |  $5\sigma$ discovery
---|---|---
HL-LHC | 880 GeV | 730 GeV
FCC-hh | 3260 GeV | 2980 GeV

##  6 Conclusions

In this paper, we studied the feasibility of discovering pair-produced down
type iso-singlet quarks $D$ at the High Luminosity LHC and the hadronic
scenario for the Future Circular Collider. The search was designed in the
$4\ell+2j$ channel, targeting the $D\rightarrow
Zd\rightarrow\ell^{+}\ell^{-}d$ decay mode, which is not accessible at the
LHC. Despite its relative low sensitivity, this channel is expected to provide
the most precise reconstruction of the $D$ quark mass. Furthermore, in case of
$D$ quark discovery through a higher sensitivity channel, $ZZ\rightarrow
4\ell$ channel would help to estimate relative branching ratios, thus leading
to a preliminary understanding of the underlying model properties. However,
extracting further information on the model would require observing and
measuring iso-singlet partners of different quark types.

The analysis consisted of a basic event selection followed by a two-step
reconstruction of the $D$ quark masses, where the $Z$ bosons were
reconstructed in the first step. A $\chi^{2}$ optimization was used for
finding the combination giving the best $D$ quark candidates. A further
selection was applied on the $\chi^{2}$ to discriminate signal events from the
background. Finally, a fit was performed on the average $D$ quark invariant
mass distribution to obtain event yields and sensitivity.

The study showed that the $5\sigma$ discovery reach for $D$ quark mass at HL-
LHC is possible, and is around 730 GeV for the full run period, while FCC-hh
can reach up to 2980 GeV, considering only the $4\ell+2j$ decay channel. It
also demonstrated that FCC-hh requires about two orders of magnitude less
integrated luminosity than HL-LHC for discovering $D$ quarks at a given mass.
Therefore searches for $E_{6}$ GUT models using $4\ell+2j$ channel would
benefit from FCC-hh. Sensitivity of FCC-hh could further be enhanced by the
addition of final states with boosted Z bosons decaying to boosted collimated
lepton jets.

As a side note, this study showed an example of how extensively the analysis
description language (ADL) concept and its runtime interpreter implementation,
CutLang, can be used to benefit particle physics analyses. This approach
allows performing the analysis algorithm steps (e.g. object definitions,
object reconstructions, histogramming) in an easy and descriptive way.

###### Acknowledgements.

SS is supported by the National Research Foundation of Korea (NRF), funded by
the Ministry of Science & ICT under contract NRF-2008-00460. This study
started during CERN Summer Student Program 2019, in which AP was a
participant. GU would like to dedicate this paper to his parents Turkan and
Olgun Unel (who both passed away during its preparation), for nurturing him
into becoming the person he is.

##  A Appendix: Explanation of the ADL implementation

We provide details on the ADL/CutLang implementation of the analysis in this
paper, in particular with the aim to clarify the implementation of composite
object reconstruction and optimization.

###  A.1 Object selection

Electron, muon and jet selection based on object properties, e.g. transverse
momentum and pseudorapidity, is expressed as :

    object goodJet
      take JET
      select Pt(JET) > 50
      select abs(Eta(JET)) < 4

    object goodEle
      take ELE
      select Pt(ELE) > 20
      select abs(Eta(ELE)) < 4

    object goodMuo
      take MUO
      select Pt(MUO) > 20
      select abs(Eta(MUO)) < 4

where `JET`, `ELE`, `MUO` are the original objects from the input event files
and `goodJet`, `goodEle` and `goodMuo` are the derived objects.

Selected electrons and muons can be combined to define the unified set of
leptons as

    object goodLep  : Union (goodEle, goodMuo)

###  A.2 Definitions

ADL allows to define aliases for event variables or reconstructed particles
through the usage of the `define` keyword. Shorthand notations for
reconstructed $Z$ bosons and $D$ quarks, optimization criteria and selection
variables are given in this section of the ADL file.

###  A.3 Event selection

Event selections in ADL are described within `region` blocks defined for each
selection region. This analysis has a single search channel, which is
described in the `DDselection` region block, which starts by selecting the
analysis final state of $4\ell+2j$ as:

    region DDselection
      select Size(goodLep)  >= 4
      select Size(goodJet)  >= 2

Subsequent optimization and selection requirements are also given in this
region as described below.

###  A.4 Z and D reconstruction

In this analysis, the particles for reconstructing $Z$ and $D$ must be
combined such that the resulting $Z$ and $D$ would best satisfy the criteria
defined by an optimization rule. Therefore the indices of the particles
combined are not known before the optimization, and are only determined after
the optimization. In CutLang, negative numbers are used for specifying such
indices of particles that would be combined through a $\chi^{2}$ optimization.
Following this approach, $Z$ reconstruction is written as

    define Zreco1 = goodLep[-1] goodLep[-1]
    define Zreco2 = goodLep[-3] goodLep[-3]

Here the lepton indices are to be determined at run time for each event
according to an optimization rule, yet to be defined. The repeated indices
stress that in combining two leptons to reconstruct a $Z$ boson, the order is
unimportant. The `chi2ZZ` variable to be minimized in order to reconstruct the
two $Z$ candidates is defined as

    define chi2ZZ = (mZ1 - 91.2)^2
                  + (mZ2 - 91.2)^2
                  + (999*{Zreco1}pdgID)^2
                  + (999*{Zreco2}pdgID)^2

Given a variable $x$ with an optimal value $v$, the operator `~=` is used to
calculate the particle combination that gives an $x$ value closest to $v$. The
optimization criteria `chi2ZZ` is finally called after the initial event
selection, with the syntax

    select chi2ZZ ~= 0

Note that, CutLang takes the PDG ID of a reconstructed object to be the sum of
the PDG ID of its constituent objects. As the constituents of `Zreco` must be
a lepton-antilepton pair, pdgID of the `Zreco` itself has to be zero. Here,
$999$ is a high enough weight factor to ensure flavour neutrality of the $Z$
boson candidates. A further requirement of $Z$ boson charge to be $0$ is also
applied as:

    select q(Zreco1) == 0
    select q(Zreco2) == 0

Next, $D$ quark candidates are reconstructed using the previously obtained $Z$
bosons and jets. As in the case for $Z$ bosons, the indices of the optimal
jets cannot be known, therefore, are written as negative indices.

    define dj1 = goodJet[-2]
    define dj2 = goodJet[-4]
    define Dreco1 = Zreco1 dj1
    define Dreco2 = Zreco2 dj2

All these expressions are then used for defining terms in the optimization
condition for $D$ reconstruction:

    define chi2mD  = ((mD1 - mD2)/mD)^2

    define  chi2PTj
        = Hstep(PTjcut -PTj1)*(PTjcut/PTj1 - 1.0)
        + Hstep(PTjcut - PTj2)*(PTjcut/PTj2 - 1.0)

    define  chidRDD = (dRDD/3.14 - 1.0)^2

These terms are added to obtain $\chi^{2}_{DD}$

    define chi2DD  = chimD + chiPTj + chidRDD

Finally, a selection criteria is applied on $\chi^{2}_{DD}$ in the
`DDselection` region as:

    select chi2DD < chi2DDcut

###  A.5 Histogramming

CutLang is designed to be a complete tool for event processing and
visualization tasks in an analysis, and therefore allows to define and fill
histograms at runtime. The CutLang syntax to plot 1D histograms of a variable
is given below :

    histo [label], "[title]", [no. of bins],
          [lower limit], [upper limit], [variable]

Following are the definitions of histograms which were eventually plotted in
Figures 3, 6, 7 and 8:

    histo hmZ1, "Z candidate1 mass (GeV)", 320,
                                     0, 3200, mZ1
    histo hmZ2, "Z candidate2 mass (GeV)", 320,
                                     0, 3200, mZ2
    histo hchi2DD, "chi2DD ", 200, 0, 10, chi2DD
    histo hmD, "D candidate mass (GeV)", 320,
                                     0, 3200, mD

## References

  * (1) G. Aad, et al., Physics Letters B 716(1), 1 (2012).  DOI https://doi.org/10.1016/j.physletb.2012.08.020.  URL http://www.sciencedirect.com/science/article/pii/S037026931200857X
  * (2) J.L. Hewett, T.G. Rizzo, Physics Reports 183(5), 193 (1989).  DOI https://doi.org/10.1016/0370-1573(89)90071-9.  URL http://www.sciencedirect.com/science/article/pii/0370157389900719
  * (3) R. Mehdiyev, S. Sultansoy, G. Unel, M. Yilmaz, The European Physical Journal C 49(2), 613 (2007).  DOI 10.1140/epjc/s10052-006-0147-9.  URL https://doi.org/10.1140/epjc/s10052-006-0147-9
  * (4) R. Mehdiyev, A. Siodmok, S. Sultansoy, G. Unel, The European Physical Journal C 54(3), 507 (2008).  DOI 10.1140/epjc/s10052-008-0550-5.  URL https://doi.org/10.1140/epjc/s10052-008-0550-5
  * (5) M. Aaboud, et al., JHEP 07, 089 (2018).  DOI 10.1007/JHEP07(2018)089
  * (6) M. Aaboud, et al., JHEP 10, 141 (2017).  DOI 10.1007/JHEP10(2017)141
  * (7) M. Aaboud, et al., JHEP 08, 048 (2018).  DOI 10.1007/JHEP08(2018)048
  * (8) M. Aaboud, et al., JHEP 08, 052 (2017).  DOI 10.1007/JHEP08(2017)052
  * (9) M. Aaboud, et al., Phys. Rev. D 98(11), 112010 (2018).  DOI 10.1103/PhysRevD.98.112010
  * (10) M. Aaboud, et al., JHEP 12, 039 (2018).  DOI 10.1007/JHEP12(2018)039
  * (11) M. Aaboud, et al., Phys. Rev. D 98(9), 092005 (2018).  DOI 10.1103/PhysRevD.98.092005
  * (12) M. Aaboud, et al., Phys. Rev. Lett. 121(21), 211801 (2018).  DOI 10.1103/PhysRevLett.121.211801
  * (13) A.M. Sirunyan, et al., Phys. Lett. B 779, 82 (2018).  DOI 10.1016/j.physletb.2018.01.077
  * (14) A.M. Sirunyan, et al., Phys. Rev. D 100(7), 072001 (2019).  DOI 10.1103/PhysRevD.100.072001
  * (15) A.M. Sirunyan, et al., JHEP 08, 177 (2018).  DOI 10.1007/JHEP08(2018)177
  * (16) A.M. Sirunyan, et al., Phys. Rev. D 97, 072008 (2018).  DOI 10.1103/PhysRevD.97.072008
  * (17) G. Unel, S. Sekmen, A.M. Toon.  CutLang: a cut-based HEP analysis description language and runtime interpreter (2019).  URL https://arxiv.org/abs/1909.10621
  * (18) S. Sekmen, G. Ünel, Comput. Phys. Commun. 233, 215 (2018).  DOI 10.1016/j.cpc.2018.06.023
  * (19) O. Çakir, M. Yilmaz, Europhysics Letters (EPL) 38(1), 13 (1997).  DOI 10.1209/epl/i1997-00528-3.  URL https://doi.org/10.1209%2Fepl%2Fi1997-00528-3
  * (20) M. Tanabashi, et al., Phys. Rev. D 98, 030001 (2018).  DOI 10.1103/PhysRevD.98.030001.  URL https://link.aps.org/doi/10.1103/PhysRevD.98.030001
  * (21) S. Sultansoy, G. Unel, Phys. Lett. B 669, 39 (2008).  DOI 10.1016/j.physletb.2008.09.022
  * (22) G. Apollinari, O. Brüning, T. Nakamoto, L. Rossi, CERN Yellow Report (arXiv:1705.08830. 5), 1 (2017).  DOI 10.5170/CERN-2015-005.1.  URL https://cds.cern.ch/record/2120673.  21 pages, chapter in High-Luminosity Large Hadron Collider (HL-LHC) : Preliminary Design Report
  * (23) O. Brüning, L. Rossi, _The High Luminosity Large Hadron Collider_ (WORLD SCIENTIFIC, 2015).  DOI 10.1142/9581.  URL https://www.worldscientific.com/doi/abs/10.1142/9581
  * (24) G. Apollinari, I. Béjar Alonso, O. Brüning, P. Fessia, M. Lamont, L. Rossi, L. Tavian.  High-Luminosity Large Hadron Collider (HL-LHC): Technical Design Report V. 0.1 (2017).  DOI 10.23731/CYRM-2017-004
  * (25) ATLAS Phase-II Upgrade Scoping Document.  Tech. Rep. CERN-LHCC-2015-020. LHCC-G-166, CERN, Geneva (2015).  URL https://cds.cern.ch/record/2055248
  * (26) D. Contardo, M. Klute, J. Mans, L. Silvestris, J. Butler, Technical Proposal for the Phase-II Upgrade of the CMS Detector.  Tech. Rep. CERN-LHCC-2015-010. LHCC-P-008. CMS-TDR-15-02, Geneva (2015).  URL https://cds.cern.ch/record/2020886
  * (27) A. Abada, et al., The European Physical Journal Special Topics 228(4), 755 (2019).  DOI 10.1140/epjst/e2019-900087-0.  URL https://doi.org/10.1140/epjst/e2019-900087-0
  * (28) M. Benedikt, et al., Future Circular Collider - European Strategy Update Documents.  Tech. Rep. CERN-ACC-2019-0005, CERN, Geneva (2019).  URL https://cds.cern.ch/record/2653674
  * (29) E.S. Group, _2020 Update of the European Strategy for Particle Physics_ (CERN Council, Geneva, 2020).  DOI 10.17181/ESU2020.  URL https://cds.cern.ch/record/2720129
  * (30) G. Brooijmans, et al.  Les Houches 2015: Physics at TeV colliders - new physics working group report (2016).  URL https://arxiv.org/abs/1605.02684
  * (31) G. Brooijmans, et al.  Les Houches 2017: Physics at TeV Colliders New Physics Working Group Report (2018).  URL https://arxiv.org/abs/1803.10379
  * (32) G. Brooijmans, et al.  Les Houches 2019 Physics at TeV Colliders: New Physics Working Group Report (2020).  URL https://arxiv.org/abs/2002.12220
  * (33) R. Brun, F. Rademakers, Nucl. Instrum. Meth. A 389, 81 (1997).  DOI 10.1016/S0168-9002(97)00048-X.  See also http://root.cern.ch/
  * (34) The lex & yacc page.  http://dinosaur.compilertools.net
  * (35) E. Boos, et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 534(1), 250 (2004).  DOI https://doi.org/10.1016/j.nima.2004.07.096.  URL http://www.sciencedirect.com/science/article/pii/S0168900204015402
  * (36) A. Pukhov, et al.  CompHEP: A Package for evaluation of Feynman diagrams and integration over multiparticle phase space (1999).  URL https://arxiv.org/abs/hep-ph/9908288v2
  * (37) J. Alwall, et al., Journal of High Energy Physics 2014(7), 79 (2014).  DOI 10.1007/JHEP07(2014)079.  URL https://doi.org/10.1007/JHEP07(2014)079
  * (38) N. Greiner, S. Höche, G. Luisoni, M. Schönherr, J.C. Winter, V. Yundin, PoS RADCOR2015(arXiv:1601.03722), 009. 11 p (2016).  DOI 10.22323/1.235.0009.  URL https://cds.cern.ch/record/2124367.  Comments: 11 pages, 6 figures, Contribution to appear in the proceedings of 12th International Symposium on Radiative Corrections (Radcor 2015) and LoopFest XIV
  * (39) R.D. Ball, et al., Eur. Phys. J. C 77(10), 663 (2017).  DOI 10.1140/epjc/s10052-017-5199-5
  * (40) T. Sjostrand, S. Mrenna, P.Z. Skands, JHEP 05, 026 (2006).  DOI 10.1088/1126-6708/2006/05/026
  * (41) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, M. Selvaggi, T.D.. collaboration, Journal of High Energy Physics 2014(2), 57 (2014).  DOI 10.1007/JHEP02(2014)057.  URL https://doi.org/10.1007/JHEP02(2014)057
  * (42) M. Selvaggi.  Beta card for hl-lhc and he-lhc studies.  https://github.com/delphes/delphes/blob/master/cards/delphes_card_HLLHC.tcl (2017)
  * (43) M. Selvaggi.  Official delphes card prepared by fcc-hh collaboration.  https://github.com/delphes/delphes/blob/master/cards/FCC/FCChh.tcl (2017)


# A Learnheuristic Approach to A Constrained Multi-Objective Portfolio
Optimisation Problem

Sonia Bullah<EMAIL_ADDRESS>Computer Science and Applied
Mathematics, University of the WitwatersrandJohannesburgSouth Africa and
Terence L. van Zyl  Institute for Intelligent Systems, University of
JohannesburgJohannesburgSouth Africa<EMAIL_ADDRESS>
(2023)

###### Abstract.

Multi-objective portfolio optimisation is a critical problem researched across
various fields of study as it achieves the objective of maximising the
expected return while minimising the risk of a given portfolio at the same
time. However, many studies fail to include realistic constraints in the
model, which limits practical trading strategies. This study introduces
realistic constraints, such as transaction and holding costs, into an
optimisation model. Due to the non-convex nature of this problem,
metaheuristic algorithms, such as NSGA-II, R-NSGA-II, NSGA-III and U-NSGA-III,
will play a vital role in solving the problem. Furthermore, a learnheuristic
approach is taken as surrogate models enhance the metaheuristics employed.
These algorithms are then compared to the baseline metaheuristic algorithms,
which solve a constrained, multi-objective optimisation problem without using
learnheuristics. The results of this study show that, despite taking
significantly longer to run to completion, the learnheuristic algorithms
outperform the baseline algorithms in terms of hypervolume and rate of
convergence. Furthermore, the backtesting results indicate that utilising
learnheuristics to generate weights for asset allocation leads to a lower risk
percentage, higher expected return and higher Sharpe ratio than backtesting
without using learnheuristics. This leads us to conclude that using
learnheuristics to solve a constrained, multi-objective portfolio optimisation
problem produces superior and preferable results than solving the problem
without using learnheuristics.

multi-objective portfolio optimisation, learnheuristics, metaheuristics,
genetic algorithms, surrogate-assisted algorithms, machine learning, deep
learning, constraints

††copyright: acmcopyright††journalyear: 2023

##  1\. Introduction

This study explores a multi-objective portfolio optimisation problem by
including numerous factors that produce above-average returns while attempting
to minimise the investor’s overall risk (Paskaramoorthy et al., 2020).

This is conducted by introducing realistic constraints into the optimisation
model alongside machine learning, deep learning (Laher et al., 2021) and
metaheuristics to solve non-convex problems. The metaheuristic algorithms
explored are Evolutionary Algorithms (EA), namely NSGA-II, R-NSGA-II, NSGA-III
and U-NSGA-III. In contrast, the machine learning models used to enhance these
metaheuristic algorithms will be that of the surrogate-assisted variants of
each metaheuristic algorithm used (Perumal and van Zyl, 2020; Stander et al.,
2020).

In the study by Ma _et al._ (Ma et al., 2021), the researchers compare machine
learning and deep learning models with portfolio optimisation models to
determine the best combination of investment strategies. Ultimately, it is
concluded that the Random Forest Model + Mean-Variance with Forecasting (MVF)
Model produces the most efficient and profitable results for practical
investments. However, the only issue is that a high turnover results in a huge
challenge for these models to overcome. A high turnover implies that the
transaction costs or transaction fees (both terms can be used interchangeably)
will be quite high. Since these costs are very much present in realistic
trading scenarios, they should not be overlooked when analysing the portfolio
optimisation problem. Hence, this study will introduce them and other
realistic factors as constraints in the portfolio optimisation problem. Adding
these additional constraints can be done by employing a hybrid approach, which
includes merging metaheuristics with machine learning and deep learning
techniques. This approach can be called ‘learnheuristics’ (Doering et al.,
2019).

Overall, this study will focus on how realistic constraints, such as
transaction costs and holding costs, can be introduced into the portfolio
optimisation problem in such a way as to account for their presence in
realistic trading scenarios as well as how learnheuristics can be used to
assist in solving a constrained, multi-objective and non-convex portfolio
optimisation problem. By introducing a unique, hybrid approach, this study
will solve the current issue of unrealistic trading scenarios in portfolio
optimisation problems while producing superior results compared to those
produced by alternative methods and benefiting the investor simultaneously.

##  2\. Background

###  2.1. Portfolio Optimisation

Portfolio optimisation forms the basis of this study and is a fundamental
problem faced by many investors (Qu et al., 2017). The portfolio optimisation
process comprises selecting the optimal portfolio from all possible portfolios
to maximise investor wealth while minimising risk. A popular way of doing this
includes using the Sharpe ratio, which can be directly optimised (Zhang et
al., 2020).

One of the major contributions to the portfolio selection problem is
Markowitz’s theory of portfolio selection (Markowitz, [n. d.]). Markowitz
derives a model that assumes that investors aim to maximise their returns for
a particular level of risk, or the converse. Many other studies have extended
Markowitz’s theory to include constraints. Kizys _et al._ (Kizys et al., 2019)
acknowledge that the complexity of this model increases significantly when
further constraints are added. However, metaheuristics can be used to overcome
such obstacles.

Doering _et al._ (Doering et al., 2019) extend Streichert _et al._ (Streichert
et al., 2004) by formulating an unconstrained portfolio optimisation problem
as per Markowitz’s mean-variance theory in such a way that minimises risk as
follows:

$Min\sum_{i=1}^{N}\sum_{j=1}^{N}\omega_{i}\omega_{j}\sigma_{ij}$

and maximises profit or returns in the following manner:

$Max\sum_{i=1}^{N}\omega_{i}\mu_{i}$

subject to:

$\sum_{i=1}^{N}\omega_{i}=1,\;\;\;\;where\;0\leq\omega_{i}\leq 1,\forall
i=1,2,...N,$

and N represents the available assets, $\mu_{i}$ represents the expected
return of asset $i$, $\sigma_{ij}$ represents the covariance between two
assets $i$ and $j$, and $\omega_{i}$ represents the decision variables of the
portfolio. This study will adjust such models to include realistic constraints
such as transaction costs.

Furthermore, portfolio optimisation can be classified as a multi-objective
optimisation problem (van Zyl et al., 2021) with the target of finding
solutions that can illustrate the trade-off present for different objectives
(Qu et al., 2017). A multi-objective portfolio optimisation problem requires
two objectives to be met (Streichert et al., 2004).

Machine learning and deep learning models can be combined with traditional
optimisation models (Ma et al., 2021) to find the best combination for
practical investments. These ideas will be explored in the sections below.

###  2.2. Machine Learning Models and Deep Learning Models for Multi-Objective
Portfolio Optimisation

Zhang _et al._ (Zhang et al., 2020) focuses primarily on using deep learning
models to evaluate portfolio optimisation problems by first providing an
objective function in which the Sharpe ratio is defined and maximised using
gradient ascent.

Ma _et al._ (Ma et al., 2021) make use of popular machine learning models,
such as Support Vector Regression (SVR) and Random Forest (RF), as well as
deep learning models, such as Deep Multilayer Perceptron (DMLP), LSTM Neural
Network and Convolutional Neural Network (CNN). Furthermore, they pair each
model with a traditional portfolio optimisation model, namely the Mean-
Variance and Omega Models.

In addition to using machine learning and deep learning models, many studies
also introduce hybrid algorithms. This can be seen in Kizys _et al._ (Kizys et
al., 2019), in which a new algorithm, _ARPO_ , is used to consider inevitable,
real-world constraints that many researchers refrain from analysing.
Similarly, van Zyl _et al._ (van Zyl et al., 2021) implement _ParDen_ to
include machine learning models that are either generative or discriminative
and use them as a surrogate to reduce costs that arise from backtesting
procedures.

Optimisation problems depend on time-intensive functions due to the problem’s
computationally expensive nature (Martin, 2021). Implementing such functions
influences the computational overhead of the algorithm, which is used to
discover new solutions every time the function is called. Therefore, a
_surrogate_ model can be incorporated to accelerate the convergence and
approximate the time-intensive functions of the models and algorithms
currently being used to solve these optimisation problems.

Surrogate-assisted optimisation satisfies the objective of producing efficient
approximate solution evaluations (ASEs), which are then used to minimise
expensive solution evaluations (ESEs), all while providing a minimal
approximation error. Overall, the convergence of the optimisation algorithm
being utilised will improve (Martin, 2021).

###  2.3. Transaction Costs And Holding Costs

Transaction costs arise from turnover; the higher the turnover, the more
transaction costs will increase. This suggests that the performance of a model
will differ with and without factoring in the relevant transaction costs (Ma
et al., 2021). Ma _et al._ (Ma et al., 2021) extends (van Zyl et al., 2021; Li
and Hoi, 2014; Doering et al., 2019; Kizys et al., 2019) and (Abate et al.,
2022), in which the concept of additional costs is discussed briefly and where
transaction costs are recognised as a major constraint in the realistic
portfolio optimisation problem.

Likewise, holding costs are also important to consider in practical investment
and trading scenarios. The simplest holding cost model accounts for charges
incurred for borrowing assets when taking a short position (Boyd et al.,
2017). However, these models can become significantly more complex under
certain conditions. Similar to transaction costs, holding costs generally
consist of non-negative values. However, negative values are not unlikely.

Trading costs and holding costs can be combined to formulate a multi-objective
optimisation problem with the inclusion of a leverage constraint as done by
van Zyl _et al._ (van Zyl et al., 2021). This is shown by:

$\text{Maximise:}\quad w^{\prime}\mu-\frac{\gamma}{2}w^{\prime}\sum
w-\gamma_{t}\phi^{trade}(w-w_{0})-\gamma_{h}\phi^{hold}(w)$ $\text{Subject
to:}\quad\lVert w\rVert_{1}\leq L^{max}$

where $w_{0}$ denotes the initial portfolio, $w-w_{0}$ represents the trades
needed to secure $w$, $\phi^{trade}$ is the trading cost function with the
trade-off parameter $\gamma_{t}$, $\phi^{hold}(w)$ is the trading cost
function with the trade-off parameter $\gamma_{h}$ and, finally, $L^{max}$
denotes the maximum leverage.

It can be noted that a fully-invested portfolio, in which a long position is
maintained, has a leverage of 1.

###  2.4. Introducing Portfolio Constraints

Most portfolio optimisation problems consider some form of risk management,
and many research papers fall short of considering stochastic optimisation
problems, which include using constraints in optimisation models (Doering et
al., 2019).

Abate _et al._ (Abate et al., 2022) consider portfolio constraints, as
considering mean-variance optimisation alone leads to unfair asset
allocations. They also implement a list of constrained strategies. These
strategies include classical constraints, flexible portfolio constraints,
norm-based constraints, variance-based constraints, tracking error volatility
constraints and beta constraints. They note that tighter constraints can
produce diversified portfolios and that choosing the proper constraints plays
a vital role in the evaluation process. Investors should typically choose
constraints that align with their target or goal.

Kizys _et al._ (Kizys et al., 2019) extend the work of Doering _et al._
(Doering et al., 2019) and Abate _et al._ (Abate et al., 2022), in which both
articles consider constraints that are unavoidable for investors to ignore in
realistic investment scenarios.

The transaction costs play a significant role when choosing constraints.
Constrained portfolio optimisation models are favourable to investors and can
be seen as an effective alternative to classic investment approaches (Doering
et al., 2019).

###  2.5. Metaheuristics for Solving Optimisation Problems

Metaheuristic algorithms can provide solutions to portfolio optimisation
problems that are non-convex. In general, the single-period optimisation
problem (SOP) can be efficiently solved given that the problem is convex (Boyd
et al., 2017). However, once this is adjusted to become a multi-objective
optimisation problem (MOP) with non-convex constraints, the problem becomes
NP-hard and is more complicated to solve. Metaheuristics play a vital role in
this case as they can be used to solve NP-hard optimisation problems and
provide exceptional solutions (Doering et al., 2019), while considering the
realistic, non-convex constraints introduced into the problem.

Doering _et al._ (Doering et al., 2019) classify metaheuristics into two
groups: nature-inspired and non-nature metaheuristics. This study will
primarily focus on a sub-class of nature-inspired metaheuristics, evolutionary
algorithms. In particular, it will focus on the Genetic Algorithm (GA).
However, this study focuses on multi-objective problems, so multi-objective
metaheuristics must be incorporated. Therefore, the various forms of genetic
algorithms, namely, NSGA-II, R-NSGA-II, NSGA-III and U-NSGA-III (van Zyl et
al., 2021), will be explored.

NSGA-II is a genetic algorithm based on non-domination and is popular for
solving multi-objective optimisation problems (Seshadri, 2006). First, the
algorithm initialises the population and considers any constraints that may be
specified. Next, the population goes through a non-dominated sort. Each
individual is then given a rank or fitness value based on the front to which
they belong. Furthermore, a crowding distance is assigned to each individual
in the population, quantifying how close the relevant individual is to their
neighbours. Next, a binary tournament selection process is used to select
parents from the population, and offspring are produced from the selected
population. The relevant population and offspring are again sorted using a
non-dominated sort procedure, and only the best individuals are chosen.

The R-NSGA-II algorithm is an extension of the NSGA-II algorithm with an
improved survival selection process (van Zyl et al., 2021). In this algorithm,
a smaller Euclidean distance between parents and the reference point is
favoured or preferred (Filatovas et al., 2017). The algorithm begins the same
way as NSGA-II, resulting in a single, merged population. The merged
population is then segmented into various non-dominated levels using non-
dominated sorting. Then each front of the population is clustered to ensure
that the solutions are diverse. The whole population size is then diminished
to match that of the parent population 111In the case of not all members being
selected, the ones with the smaller Euclidean distances between the member
itself and the reference point will be conserved. (Filatovas et al., 2017).
The difference between this algorithm and NSGA-II is that solutions are based
on the rank or fitness value while segmenting the front (van Zyl et al.,
2021).

Unlike NSGA-II, which uses selection based on rank or fitness and crowding
distance, NSGA-III uses reference directions to ensure diversity amongst the
solutions (Seada and Deb, 2014). The NSGA-III algorithm begins with a random
population and disperse reference points on a hyper-plane and aims to find
only one population member per reference point. The population is segregated
into various non-dominated levels, and then non-dominated sorting is
performed, just as the NSGA-II algorithm does. After this, an offspring
population is produced from which a combined population is formed. Next,
points are iteratively selected until a point is reached at which all
solutions from a complete front cannot be selected. The final front will not
be able to be fully selected. However, only a few solutions from this front
need to be selected in this case. Then, each population member and member of
the final front is linked to a reference point obtained from the minimum
perpendicular distance of each member. Next, a niching strategy is implemented
and followed to obtain members of the final front corresponding to the
population’s marginalised reference points. Population members associated with
these points are preferred. The algorithm attempts to find one population
member corresponding to each reference point close to the Pareto-optimal front
(Seada and Deb, 2014).

The population will be segmented into non-dominated fronts using a non-
dominated sorting process for a multi-objective optimisation problem. U-NSGA-
III uses a niched tournament selection process as this has a superior
performance compared to random selection (van Zyl et al., 2021). This places
weight on non-dominated solutions over dominated solutions and solutions
nearer to the reference directions than those farther away. Thereafter, the
algorithm follows the same process as that of the NSGA-II algorithm (Seada and
Deb, 2014).

###  2.6. The Learnheuristic Approach

_‘Learnheuristics’_ is a hybrid approach that combines metaheuristics with
machine learning or statistical learning models (Doering et al., 2019). This
study will use a learnheuristic approach to pair each of the four multi-
objective genetic algorithms with a surrogate model. The surrogate-assisted
algorithms will then be used to evaluate the portfolio optimisation problem.

Calvet _et al._ (Calvet et al., 2017) explain that the hybridisation of
machine learning and metaheuristics can be split into two classes, namely one
in which machine learning is used to enhance metaheuristics and the other
class being one where metaheuristics are used to improve the machine learning
models or techniques used.

Using machine learning to enhance metaheuristics can be split into
specifically-located and global-level hybridisations. Specifically-located
hybridisations take on approaches such as fine-tuning parameters, quick
initialisation, inexpensive evaluation techniques, population management,
introducing the knowledge in operators and, lastly, using machine learning
techniques as local searches. On the other hand, global hybridisations focus
on reducing the search space, algorithm selection, using hyperheuristics
(which are search methods used for choosing heuristics that will provide
solutions to computational search problems), cooperative strategies and,
finally, using new types of metaheuristics.
Conversely, using metaheuristics to improve machine learning techniques spans
topics such as classification, regression, clustering and rule mining.

Taking on a learnheuristic approach implies that either one or both of the
classes mentioned above be implemented. Table 1 indicates that while the
related work either acknowledges or implements methods that consider some form
of realistic constraints, only the survey presented by Doering _et al._
(Doering et al., 2019) focuses on metaheuristics to solve multi-objective
optimisation problems. Therefore, this study will combine techniques outlined
in the related work and attempt to enhance metaheuristics by using machine
learning techniques and producing improved results.

Table 1. A summary of the work related to this study. Author |  Optimisation |  Constraints |  Algorithms |  Advantages
---|---|---|---|---
Ma _et al._ |  Single- & Multi-objective |  With constraints |  Machine learning & Deep learning |  Provides solutions to NP-hard optimisation problems; produces high-quality solutions.
Zhang _et al._ |  Single-objective |  With constraints |  Deep learning |  Improved performance; method is rational and practical.
Doering _et al._ |  Single-objective |  With and without constraints |  Metaheuristic |  Yields a high turnover.

##  3\. Methodology

To test whether or not a learnheuristic approach to a multi-objective
portfolio optimisation problem will produce higher returns while minimising
the risk of an investment, we generate a hybrid learnheuristic algorithm and
compare the results to those of the baseline algorithms. In the past, research
surrounding multi-objective portfolio optimisation, such as the research
conducted by Ma _et al._ (Ma et al., 2021), has been explored without
realistic constraints, such as transaction costs. Therefore, while attempting
to maximise returns and minimise risk concurrently, our study will factor in
these constraints in a learnheuristic approach similar to the process outlined
in the study by Calvet _et al._ (Calvet et al., 2017).

The implementation of this study can be separated into five main
interdependent sections. The sections will all be described in detail below
and are summarised in Fig. 1.

Figure 1. A summary of the methodology followed in this study.

###  3.1. Data Cleaning And Exploration

The S&P 500 data was pre-processed to ensure correct datatypes and eliminate
redundant attributes. Duplicates, missing values and null values were also
dropped. We used only the top ten stocks regarding expected returns as our
portfolio to conduct the rest of the study.

###  3.2. Forecasting

To forecast the data, some form of noise needs to be added. This noise
represents an error that will increase or decrease predictions. This makes it
difficult to make accurate predictions about the data. We use a mean of 0 and
a standard deviation of 0.1 to acquire random samples from a normal
distribution. This, in turn, outputs an array containing the noise, which is
now added to the data, resulting in forecasted data values. Lastly, any null
values need to be dropped before the portfolio optimisation.

###  3.3. Genetic Algorithms

In this study, the four genetic algorithms are needed to satisfy two main
objectives, namely, optimising the portfolio and then finding the optimal
constraint value simultaneously and, secondly, producing the best set of
weights that can be used to perform backtesting. All genetic algorithms are
implemented using the _Pymoo_ (Blank and Deb, 2022a) library, and the
portfolio optimisation method followed is dependent on _PyPortfolioOpt_ (Blank
and Deb, 2020).

First, an objective function is adjusted to perform portfolio optimisation on
the data. This is done by obtaining the forecasted dataset’s expected returns
and sample covariance. The expected returns are computed using the Capital
Asset Pricing Model

$R_{i}=R_{f}+\beta_{i}(E(R_{m})-R_{f}),$

where $R_{f}$ represents the risk-free rate, $\beta_{i}$ is the beta (or
sensitivity) of asset _i_ and $(E(R_{m})-R_{f})$ the market risk premium. In
addition, the sample covariance is obtained using the Ledoit-Wolf covariance
shrinkage estimator

$\hat{\Sigma}=\alpha F+(1-\alpha)S,$

where $\alpha$ represents the shrinkage constant, $F$ the a structured
estimator and $S$ a sample covariance matrix.

This, in turn, allows us to find optimal portfolios and produce an efficient
frontier once the results are correctly plotted. Furthermore, it should be
noted that this method will solve an unconstrained portfolio. However, the
focus of this study is on constrained portfolios. As such, we add a constraint
to the optimisation problem.

Once an efficient frontier has been produced, we proceed by adding a
constraint to the efficient frontier. Instead of adding an arbitrary
constraint with a random value, we search for the best possible parameter.
Once this has been found, the Sharpe ratio is maximised, and the overall
portfolio performance is returned.

The genetic algorithms are executed by randomly generating an initial
population and assigning a fitness score to each individual. Next, individuals
go through the selection and reproduction stages before crossover and mutation
occur, in which new genetic information is added to the child population.
Finally, generational replacement occurs, and then the algorithm terminates
based on a stopping condition.

Once each genetic algorithm returns the array of portfolio performance results
obtained, we adjust the array by removing the Sharpe Ratio values from each
sub-array as we do not want a three-dimensional plot. This leaves us with sub-
arrays of size two, which contain the optimal portfolio’s return and
volatility (or risk). These results are then plotted.

This process is run ten times on a population of 12 over 30 generations. We
compare results by plotting all four genetic algorithms on the same axis. To
compare the success of each algorithm, we note the average time taken for each
algorithm to run to completion and the average hypervolume produced by the
algorithm. Firstly, we begin a timer each time an algorithm starts and end it
once the algorithm completes its last generation. The algorithms’ average time
(in seconds) is returned at the end. After each algorithm runs, we use the
hypervolume function from the _Pymoo_ (Blank and Deb, 2022a) library. A
hypervolume is a performance indicator which uses a reference point to
calculate the area influenced by the solution set concerning a specific
reference point (Blank and Deb, 2022a). It should be noted that higher
hypervolumes suggest better-performing algorithms. Once we have the relevant
hypervolumes for each generation of each algorithm, we calculate the average
hypervolume from all ten runs and compare the algorithms to each other to
determine how well they perform over the total number of generations.

Within the objective function, we also add an objective to the efficient
frontier, which aims to adjust the gamma value. This parameter for L2
regularisation can be increased when requiring further non-negligible weights.
The genetic algorithms search for the best possible value of gamma, which, in
turn, produces the optimal set of weights that are used to perform
backtesting. This process is elaborated on in the next section. Like the first
case, once the set of weights for each stock is found, the portfolio is
optimised, and the maximum Sharpe ratio is found. The final weights of each
stock of the optimised portfolio are then added to an array, which is utilised
in the backtesting procedure.

###  3.4. Learnheuristics: Surrogate-Assisted Genetic Algorithms

Once our genetic algorithms run, we adjust each function to include a
surrogate model. We call the GPSAF implementation from the _Pysamoo_ (Blank
and Deb, 2022b) framework to do this. Furthermore, we manually fine-tune the
_alpha_ , _beta_ , _n_max_doe_ , _n_max_infills_ , _n_offsprings_ and _seed_
parameters to produce an optimal convergence rate for each algorithm. We
continue adjusting these parameters until the surrogate-assisted genetic
algorithms are ‘learning’ throughout the iterations. We eventually produce an
improved solution compared to the non-surrogate-assisted genetic algorithms.
Again, these learnheuristic algorithms are run ten times on a population size
of 12 over 30 generations, with the average time taken and average hypervolume
being recorded at the end. We compare results by plotting each genetic
algorithm alongside the surrogate-assisted version of the algorithm on the
same axis and take note of the difference in convergence rate, average
hypervolume and the average time taken for the algorithms to run to
completion.

###  3.5. Backtesting

We run simulations or backtests to help optimise our strategy. To do this, we
call the _Backtest_ function from the _portfolio-backtest_ 222The
documentation can be found on: https://pypi.org/project/portfolio-backtest/
library. The weights obtained by the genetic algorithm are substituted into
the function and the target return, target variance, start period and end
period. We first run a backtest using the function’s built-in method of
determining the weights for each asset and then use each of the four
surrogate-assisted genetic algorithms to search for the optimal set of weights
that can then be passed into the backtesting function. Lastly, we record the
data, namely the annual volatility, expected annual return and Sharpe ratio
outputted from each backtest and compare the results.

All algorithms and implementations were run on a Quad-Core Intel Core i5
running macOS 12.3. The complete code, datasets and results can be found on
the GitHub repository.

##  4\. Results and Discussion

In this section, we present the results of the four multi-objective genetic
algorithms used to optimise the portfolio and then consider the limitations
and future work.

###  4.1. Results

Based on the results presented in Fig. 2(a) below, it is evident that the
NSGA-III algorithm runs to completion in the shortest amount of time, followed
by the U-NSGA-III algorithm and then the R-NSGA-II algorithm, while the NSGA-
II algorithm takes the longest time to run. However, given that the NSGA-II
algorithm produces the highest hypervolume followed by the R-NSGA-II, U-NSGA-
III, and, lastly, NSGA-III algorithms, it suggests that the best performing
algorithm is NSGA-II.

Algorithm Generations Population HV Time (sec) NSGA-II 30 12 0.7225 1.0115
R-NSGA-II 30 12 0.5790 0.7941 NSGA-III 30 12 0.4986 0.7557 U-NSGA-III 30 12
0.5126 0.7681

(a) The average hypervolume and average time taken (in seconds) for each
algorithm to run over 30 generations with a population size of 12. (a) A
comparison of all four genetic algorithms with a population size of 12 over 30
generations.

Figure 2. A comparison of the different genetic algorithms.

Based on Fig. 2(a), it can be seen that the NSGA-II algorithm is the first to
converge and produces a smoother line plot. The remaining three algorithms
contain more spikes in their respective line plots, however, the general
outcome suggests that the hypervolume increases as the number of generations
increases.

Next, we present the results of each of the four learnheuristic algorithms (or
surrogate-assisted genetic algorithms). Fig. 2(b) indicates that the
surrogate-assisted R-NSGA-III algorithm runs to completion in the shortest
amount of time, followed by the U-NSGA-III algorithm and then the NSGA-III
algorithm, while the NSGA-II algorithm takes the longest time to run. However,
similar to the non-surrogate variant, the NSGA-II algorithm yields the highest
hypervolume, thus making it the most efficient algorithm.

Algorithm Generations Population HV Time (sec) SA_NSGA-II 30 12 0.7543 18.0166
SA_R-NSGA-II 30 12 0.7181 15.2820 SA_NSGA-III 30 12 0.6403 16.6078 SA_U-NSGA-
III 30 12 0.7093 16.4171

(b) The average hypervolume and average time taken (in seconds) for each
surrogate-assisted metaheuristic algorithm to run over 30 generations with a
population size of 12.

Algorithm | Volatility (%) | Return (%) | SR | Time (sec)
---|---|---|---|---
No GA | 22.4 | 26.6 | 1.10 | -
SA_NSGA-II | 18.2 | 55.6 | 2.94 | 40.7051
SA_R-NSGA-II | 18.2 | 55.6 | 2.94 | 43.0664
SA_NSGA-III | 18.2 | 55.6 | 2.94 | 29.5306
SA_U-NSGA-III | 18.2 | 55.6 | 2.94 | 28.7683

(c) The annual volatility, expected annual return and Sharpe ratio (SR) of a
tangency portfolio generated from backtests as well as the time taken (in
seconds) for each algorithm to run.

Figure 3. A comparison of the different surrogate-assisted metaheuristic
algorithms.

Fig. 4(a), Fig. 4(b), Fig. 4(c) and Fig. 4(d) below each indicate that the
learnheuristic algorithms produce higher hypervolumes than the baseline
algorithms. By comparing Fig. 2(a) and Fig. 2(b), we observe that while the
learnheuristic algorithms take significantly longer to run to completion, they
produce much higher hypervolumes and converge faster while also creating
smoother plots. As Blank _et al._ (Blank and Deb, 2020) mentioned, the
learnheuristic algorithms accelerate the convergence of the multi-objective
metaheuristic algorithms used to solve the optimisation problem.

(a) A comparison of the NSGA-II and surrogate-assisted NSGA-II algorithms.

(b) A comparison of the R-NSGA-II and surrogate-assisted R-NSGA-II algorithms.

(c) A comparison of the NSGA-III and surrogate-assisted NSGA-III algorithms.

(d) A comparison of the U-NSGA-III and surrogate-assisted U-NSGA-III
algorithms.

Figure 4. A comparison of the different algorithms.

Lastly, we run simulations or backtests on our data and record the annual
volatility, expected annual return and Sharpe ratio generated for a tangency
portfolio of ten assets. Based on Fig. 2(c), we see that the tangency
portfolio for each of the four learnheuristic algorithms outperforms the
backtest run without a learnheuristic algorithm. The portfolios that run by
utilising the surrogate-assisted genetic algorithms to find the optimal
weights for asset allocation yield higher returns for a lower risk percentage
and produce greater Sharpe ratios. In these cases, the expected annual returns
and Sharpe ratios are more than double compared to the backtest run without
using a genetic algorithm to find the optimal asset weights. In the study by
Zhang _et al._ (Zhang et al., 2020), the Sharpe ratio is maximised within the
objective function before deep learning models evaluate the portfolio
optimisation problem. Similarly, in this study, we solve the optimisation
problem by maximising the Sharpe ratio in the objective function, which is
then called by each metaheuristic algorithm.

Although each of the four learnheuristic algorithms produces optimal results,
Fig. 2(c) indicates that the surrogate-assisted U-NSGA-III algorithm runs to
completion in the shortest time, followed by NSGA-III, NSGA-II and, finally,
R-NSGA-II. Therefore, the surrogate-assisted U-NSGA-III algorithm is the most
efficient learnheuristic for backtesting purposes. Overall, it is evident that
the learnheuristic algorithms produce superior results compared to solving the
portfolio optimisation problem without a learnheuristic approach.

###  4.2. Limitations

Despite the findings above, it should be noted that the learnheuristic
algorithms are sensitive to specific parameters, such as the _n_offsprings_ ,
_alpha_ , _beta_ , _n_max_doe_ , _n_max_infills_ and _seed_ parameters. If one
of these parameters is slightly off, the algorithms may prematurely terminate.
Furthermore, the algorithms fail to consider large portfolios as increasing
the number of stocks in the portfolio amplifies parameter sensitivity. Hence
this study only considers the results produced from a portfolio of ten assets.

###  4.3. Future Work

####  4.3.1. Implementing Search Algorithms

Since the parameters of the algorithms explored are pretty sensitive, search
algorithms, such as a grid search or a random search, can be utilised to find
the most appropriate parameters to pass into the algorithms.

####  4.3.2. Increasing The Size of The Portfolio

Once stable parameters have been obtained, future research can analyse much
larger portfolios, such as the complete S&P 500 index, instead of only the top
ten considered in the current study.

####  4.3.3. Applying Alternative Multi-Objective Metaheuristics

Future research can be extended to the current study by implementing different
multi-objective metaheuristic algorithms to solve the same problem identified
in this study. Possible algorithms to consider include Particle Swarm
Optimisation (PSO), Differential Evolution (DE) or Evolutionary Strategies
(ES).

##  5\. Conclusion

Portfolio optimisation problems must be solved while considering realistic
constraints inevitable in practical trading and investment scenarios. This
study introduces constraints, such as transaction costs, into a multi-
objective portfolio optimisation problem using taking on a hybrid approach
termed learnheuristics, which uses machine learning techniques, such as
surrogate models, to enhance the metaheuristic algorithms being used to solve
an optimisation problem, which is non-convex.

After implementing each of the four learnheuristic algorithms, it is evident
that the NSGA-II learnheuristic algorithm is the best-performing algorithm as
it reaches a higher hypervolume than the other three learnheuristic algorithms
implemented. Furthermore, the hypervolume of each learnheuristic algorithm
surpasses the hypervolume of the equivalent non-surrogate variant, which is
our baseline algorithm, despite taking much longer for the learnheuristic
algorithms to run to completion.

Based on the findings, we can conclude that the learnheuristic algorithms
outperform the baseline algorithms regarding hypervolume and convergence rate.
However, they take much longer to run to completion. Overall, using
learnheuristics to obtain the weights of each asset for backtesting purposes
results in a lower annual risk, higher annual expected return and higher
Sharpe ratio compared to backtesting without using learnheuristics. Therefore,
it can be concluded that learnheuristics produce superior and preferable
results when solving a constrained, multi-objective portfolio optimisation
problem.

This study also consists of a few limitations. Firstly, the learnheuristic
algorithms are sensitive to many input parameters, and since we aim to reach
optimal solutions, finding accurate parameters is vital. Secondly, the
implementation of the study is only run on a portfolio consisting of ten
stocks, as adding more stocks to the portfolio does not result in optimal
solutions due to the sensitivity of the learnheuristic parameters. Therefore,
in future research, implementing a search algorithm, such as a random or grid
search, to tune the parameters would be beneficial. This, in turn, may also
assist with allowing the learnheuristic algorithms to consider portfolios with
more than ten assets and still produce optimal results. Furthermore,
alternative multi-objective metaheuristic algorithms can also be considered in
future research.

## References

  * (1)
  * Abate et al. (2022) Guido Abate, Tommaso Bonafini, and Pierpaolo Ferrari. 2022\.  Portfolio Constraints: An Empirical Analysis.  _International Journal of Financial Studies_ 10, 1 (2022), 9\.
  * Blank and Deb (2020) J. Blank and K. Deb. 2020.  pymoo: Multi-Objective Optimization in Python.  _IEEE Access_ 8 (2020), 89497–89509.
  * Blank and Deb (2022a) Julian Blank and Kalyanmoy Deb. 2022a.  pysamoo: Surrogate-Assisted Multi-Objective Optimization in Python.  arXiv:2204.05855 [cs.NE]
  * Blank and Deb (2022b) Julian Blank and Kalyanmoy Deb. 2022b.  pysamoo: Surrogate-Assisted Multi-Objective Optimization in Python.  _arXiv preprint arXiv:2204.05855_ (2022).
  * Boyd et al. (2017) Stephen Boyd, Enzo Busseti, Steve Diamond, Ronald N Kahn, Kwangmoo Koh, Peter Nystrup, Jan Speth, et al. 2017\.  Multi-period trading via convex optimization.  _Foundations and Trends® in Optimization_ 3, 1 (2017), 1–76.
  * Calvet et al. (2017) Laura Calvet, Jésica de Armas, David Masip, and Angel A Juan. 2017.  Learnheuristics: hybridizing metaheuristics with machine learning for optimization with dynamic inputs.  _Open Mathematics_ 15, 1 (2017), 261–280.
  * Doering et al. (2019) Jana Doering, Renatas Kizys, Angel A Juan, Angels Fito, and Onur Polat. 2019.  Metaheuristics for rich portfolio optimisation and risk management: Current state and future trends.  _Operations Research Perspectives_ 6 (2019), 100121.
  * Filatovas et al. (2017) Ernestas Filatovas, Algirdas Lančinskas, Olga Kurasova, and Julius Žilinskas. 2017\.  A preference-based multi-objective evolutionary algorithm R-NSGA-II with stochastic local search.  _Central European Journal of Operations Research_ 25, 4 (2017), 859–878.
  * Kizys et al. (2019) Renatas Kizys, Angel A Juan, Bartosz Sawik, and Laura Calvet. 2019.  A biased-randomized iterated local search algorithm for rich portfolio optimization.  _Applied Sciences_ 9, 17 (2019), 3509.
  * Laher et al. (2021) Siddeeq Laher, Andrew Paskaramoorthy, and Terence L Van Zyl. 2021\.  Deep learning for financial time series forecast fusion and optimal portfolio rebalancing. In _2021 IEEE 24th International Conference on Information Fusion (FUSION)_. IEEE, 1–8.
  * Li and Hoi (2014) Bin Li and Steven CH Hoi. 2014.  Online portfolio selection: A survey.  _ACM Computing Surveys (CSUR)_ 46, 3 (2014), 1–36.
  * Ma et al. (2021) Yilin Ma, Ruizhu Han, and Weizhong Wang. 2021.  Portfolio optimization with return prediction using deep learning and machine learning.  _Expert Systems with Applications_ 165 (2021), 113973.
  * Markowitz ([n. d.]) HM Markowitz. [n. d.].  Portfolio Selection (1952) 7 J.
  * Martin (2021) Robert Andrew Martin. 2021\.  PyPortfolioOpt: portfolio optimization in Python.  _Journal of Open Source Software_ 6, 61 (2021), 3066\.  https://doi.org/10.21105/joss.03066
  * Paskaramoorthy et al. (2020) Andrew B Paskaramoorthy, Tim J Gebbie, and Terence L van Zyl. 2020\.  A framework for online investment decisions.  _Investment Analysts Journal_ 49, 3 (2020), 215–231.
  * Perumal and van Zyl (2020) Rylan Perumal and Terence L van Zyl. 2020.  Surrogate Assisted Methods for the Parameterisation of Agent-Based Models. In _2020 7th International Conference on Soft Computing & Machine Intelligence (ISCMI)_.
  * Qu et al. (2017) BY Qu, Q Zhou, JM Xiao, JJ Liang, and PN Suganthan. 2017.  Large-scale portfolio optimization using multiobjective evolutionary algorithms and preselection methods.  _Mathematical Problems in Engineering_ 2017 (2017).
  * Seada and Deb (2014) Haitham Seada and Kalyanmoy Deb. 2014.  U-NSGA-III: A unified evolutionary algorithm for single, multiple, and many-objective optimization.  _COIN report_ 2014022 (2014).
  * Seshadri (2006) Aravind Seshadri. 2006\.  A fast elitist multiobjective genetic algorithm: NSGA-II.  _MATLAB Central_ 182 (2006).
  * Stander et al. (2020) L Stander, M Woolway, and TL van Zyl. 2020.  Extended surrogate assisted continuous process optimisation. In _2020 7th International Conference on Soft Computing & Machine Intelligence (ISCMI)_. IEEE, 275–279.
  * Streichert et al. (2004) Felix Streichert, Holger Ulmer, and Andreas Zell. 2004\.  Evolutionary algorithms and the cardinality constrained portfolio optimization problem.  In _Operations Research Proceedings 2003_. Springer, 253–260.
  * van Zyl et al. (2021) Terence L van Zyl, Matthew Woolway, and Andrew Paskaramoorthy. 2021\.  Parden: Surrogate assisted hyper-parameter optimisation for portfolio selection. In _2021 8th international conference on soft computing & machine intelligence (ISCMI)_. IEEE, 101–107.
  * Zhang et al. (2020) Zihao Zhang, Stefan Zohren, and Stephen Roberts. 2020\.  Deep learning for portfolio optimization.  _The Journal of Financial Data Science_ 2, 4 (2020), 8–20.