FTUV-17-0508.7937

IFIC/17-24

Direct Bounds on Heavy Top-Like Quarks
With Standard and Exotic Decays

Mikael Chala

Departament de Física Tèorica, Universitat de València and IFIC, Universitat
de València-CSIC, Dr. Moliner 50, E-46100 Burjassot (València), Spain

###### Abstract

Heavy vector-like quarks with electric charge $Q=2/3$ (also called heavy tops)
appear naturally in many extensions of the Standard Model. Although these
typically predict the existence of further particles below the TeV scale,
direct searches for heavy tops have been performed assuming that they decay
only into SM particles. The aim of this paper is to overcome this situation.
We consider the most constraining experimental LHC searches for vector-like
quarks, including analyses of the 36 fb-1 of data collected in the latest run
at 13 TeV of center of mass energy, as well as searches sensitive to heavy
tops decaying into a new scalar, $S$. Combining all these, we derive bounds
for arbitrary values of the heavy top branching ratios. A simple code that
automatizes this process is also provided. At the physics level, we
demonstrate that bounds on heavy tops are not inevitably weaker in the
presence of new light scalars. We find that heavy tops with masses below $\sim
900$ GeV are excluded by direct searches, independently of whether they decay
into $Zt,Ht,Wb$ or $St$ (with $S$ giving either missing energy of bottom
quarks) or into any combination of them.

###### Contents

  1. 1 Introduction
  2. 2 Assumptions
  3. 3 Searches for vector-like quarks with standard decays
    1. 3.1 $T\rightarrow Zt$
    2. 3.2 $T\rightarrow Ht$
    3. 3.3 $T\rightarrow Wb$
  4. 4 Searches for vector-like quarks with exotic decays
    1. 4.1 $T\rightarrow St,S\rightarrow b\overline{b}$
    2. 4.2 $T\rightarrow St,S\rightarrow E_{T}^{\text{miss}}$
  5. 5 Final results
  6. 6 Conclusions
  7. A Tables
  8. B Software

##  1 Introduction

Heavy vector-like quarks with electric charge $Q=2/3$ and masses around the
TeV scale appear naturally in many extensions of the Standard Model (SM). We
will refer to them as heavy tops or top partners interchangeably. Since they
are colored, they can be produced in pairs via QCD interactions; their
signatures depending only on the way they decay. These have been widely
studied assuming that these resonances decay only into SM particles [1, 2, 3],
namely into $Zt$, $Ht$ and $Wb$ (charged-conjugated decays are also
understood). However, many models of new physics predict also the existence of
further particles, into which the heavy tops could decay. Consequently, these
models can not be easily constrained in light of current analyses.

This is the case, for example, of composite Higgs models (CHM) [4, 5, 6] with
partial compositeness [7]. Apart from the minimal setup, based on the
symmetry-breaking pattern $SO(5)/SO(4)$ [8], all realizations of the composite
Higgs paradigm contain new scalars of electroweak (EW) mass. In particular
singlets, as for example $SO(6)/SO(5)$ [9], $SO(7)/SO(6)$ [10], $SO(7)/G_{2}$
[11], $SO(6)/SO(4)$ [12] or $SU(5)/SO(5)$ [13], among others. Note also that
the minimal CHM is in no way preferred over the non-minimal realizations. More
the contrary: non-minimal CHMs can provide dark matter candidates [14, 15, 11,
16, 17, 18], explanations for the observed baryon-anti-baryon asymmetry [19,
10] and feasible UV completions [20, 21, 22, 23, 24, 25]. Therefore, the top
partner phenomenology can be totally different from what current experimental
analyses consider. As it has been anticipated in previous references [26, 27],
new decay modes must be also taken into account.

Recent works [27, 28, 29, 30, 31] have made a first attempt to address this
question by exploring the signatures of new channels. However, this approach
is valid only when new decays are dominant and hence of little help. Likewise,
standard analyses apply only when the top-partners decay mostly into SM
particles, namely into $Ht,Zt$ and $Wb$. Overall, top-partners with both
standard and exotic decays can be bounded by no means in regard of present
analyses. Our aim in what follows is to fill this gap by extending previous
efforts in two ways. First, we consider all possible sets of branching ratios
for the top partners, counting not only new channels but also elusive decays,
i.e. those that evade current searches. They have the sole (but so far rather
unexplored) effect of making the branching ratios of all observable topologies
not add to 1 (see ref. [32] for previous considerations of this possibility).
And second, we include LHC data acquired at both $8$ [33] and $13$ [34, 35,
36, 37] TeV of center of mass energy (c.m.e.). All results are given in terms
of tables and plots than can be trivially used to constrain generic top-like
resonances. To simplify further this task, we also provide a very simple code
that implements our findings. It can be found in
http://github.com/mikaelchala/vlqlimits. We intend to update this program with
the inclusion of incoming analyses and other signatures of top and also bottom
partners.

Evidently, heavy tops can decay into several new channels without conflicting
with the SM gauge symmetries. We do not explore all of them in this article,
nor are they included in the mentioned code in its current form. Instead, we
restrict ourselves to the standard channels plus the top-partners decaying
into a neutral scalar and a top quark. We assume that the former either decays
into bottom quarks or it escapes detectors. As a matter of fact, these
production modes are expected to dominate over others in concrete CHMs. This
is discussed in section 2. The rest of the article is structured as follows.
In section 3 we describe the status of experimental searches that bound the
top-partner decays into SM particles. We also explain how these can be
combined without necessarily recasting the corresponding analyses. In section
4 we concentrate on searches for heavy tops with exotic decays. Different
analyses are recast and applied to simulated events in these topologies and
their expected collider signatures are obtained. In section 5 we describe how
to combine all previous outcomes to bound heavy tops of several masses. Plots
for a broad set of branching ratios are provided. It is shown that bounds on
heavy tops are not inevitably weaker in the presence of new light scalars.
Finally, in appendix B we comment on how to obtain the code that implements
the aforementioned results.

##  2 Assumptions

Inspired by non-minimal CHMs, we assume that the relevant degrees of freedom
at the TeV scale involve the SM particles, some new fermionic resonances and
at least one real neutral scalar, $S$. The latter, as well as the Higgs, $H$,
are supposed to be composite objects resulting from the confinement of a new
strongly interacting sector at some scale $f\sim$ TeV. We will denote by
$m_{H}$, $m_{S}$ and $M$ the masses of $H$, $S$ and the fermionic resonances,
respectively. The phenomenological study presented in this article is based on
the following assumptions.

Assumption 1: The top partners, $T$, are among the lightest fermionic
resonances. This has been explained at length in the CHM literature. A sketch
of the argument reads as follows. According to the partial compositeness
paradigm [7], the elementary fermions, $q$, mix with vector-like resonances,
$Q$, through the mixing Lagrangian $\sim\Delta\overline{q}Q+\text{h.c.}$ The
latter fermions are the only which, in turn, couple to the Higgs boson,
because they are all composite. In the physical basis, however, the massless
particles (to be identified with the SM fermions) obtain Yukawa interactions
of the order $\sim(\Delta/M)^{2}$. The smallest $M$ is then that of the
fermionic resonance coupling to the SM fermion with the largest Yukawa. Namely
the top quark, as we claim. More quantitative analyses can be found, for
example, in refs. [38, 39, 40, 41, 42, 43].

Assumption 2: $m_{S}>m_{H}$, being in general of a few hundreds GeV. The
scalar potential for $H$ and $S$ can be estimated to be [44]

$V\sim\frac{3}{(4\pi)^{2}}f^{2}M^{2}\left[-\alpha\frac{H^{2}}{f^{2}}+\beta\frac{H^{4}}{f^{4}}+\gamma\frac{S^{2}}{f^{2}}\right]+...$
(2.1)

where $\alpha,\beta$ and $\gamma$ are dimensionless coefficients of order
$\mathcal{O}(1)$, and the ellipsis stand for terms with higher powers of $S$.
Experimental limits from Higgs searches [45] and EW precision data [46] impose
$f\gtrsim 900$ GeV, while natural arguments suggest that this can not be much
larger. Thus, let us take $M\sim f\sim 1$ TeV. It turns out that $m_{S}$
ranges from $\sim 150$ GeV (for $\gamma\sim 0.5$) to $\sim 200$ GeV (for
$\gamma\sim 1$). We can further strengthen this conclusion by focusing on
particular models with computable $m_{S}$. Among these, we find CHMs in which
$S$ is stable. The first such a model we consider is based on the coset
$SO(6)/SO(5)$, with the third-generation $q_{L}$ and $t_{R}$ mixing with
composite resonances transforming in the $\mathbf{20}$ and the $\mathbf{1}$
representations of $SO(6)$, respectively. The scalar potential at the leading
order adopts the form [26, 44]

$V=f^{2}\bigg{[}c_{1}-\frac{7}{4}c_{2}\bigg{]}H^{2}+(c_{2}-c_{1})H^{4}-c_{2}f^{2}S^{2}+\cdots$
(2.2)

where $c_{1}$ and $c_{2}$ are the only free parameters not constrained by the
symmetries. These can be traded by the known values of the Higgs mass term,
$\mu_{H}^{2}$, and the quartic coupling, $\lambda_{H}\sim 0.13$:

$V=\frac{1}{2}\mu_{H}^{2}H^{2}+\frac{1}{4}\lambda_{H}H^{4}+\frac{1}{3}f^{2}\left[1-2\frac{v^{2}}{f^{2}}\right]\lambda_{H}S^{2}+\cdots$
(2.3)

with $v\sim 246$ GeV the EW VEV. If we take $f\sim 1$ TeV, we obtain
$m_{S}\sim 300$ GeV. The second model we consider is based on the larger coset
$SO(7)/G_{2}$, with the third generation $q_{L}$ and $t_{R}$ mixing with the
$\mathbf{35}$ and the $\mathbf{1}$ representations of $SO(7)$, respectively.
The scalar potential can then be written as [18]

$V=\mu_{H}^{2}H^{2}+\lambda_{H}H^{4}+\frac{1}{3}f^{2}\bigg{[}1-\frac{8}{3}\frac{v^{2}}{f^{2}}\bigg{]}\lambda_{H}S^{2}+\cdots$
(2.4)

Again, for $f\sim 1$ TeV we obtain $m_{S}\sim 300$ GeV.

Assumption 3: If $S$ is not stable, it decays mostly into bottom quarks. $S$
is a composite particle, and therefore interacts stronger with the heavier
particles. Thus, with tops and bottoms. However, in light of the discussion
above, it is clear that there exists a natural regime in which $S$ has no
kinematics space to decay into $t\bar{t}$, nor into $HH$. Besides, the latter
can be further suppressed if $S$ is CP odd.

As a consequence of these three assumptions, the heavy tops can only decay
sizably into $Zt,Ht,Wb$ and $St$ (with $S$ either decaying into bottom quarks
or escaping detection). With this in mind, in next sections we obtain bounds
on $M$ for arbitrary branching ratios and for different values of $m_{S}$ in
between 100 and 400 GeV. We focus on the pair production channel, because it
is driven by model-independent QCD interactions. In this respect, however, two
comments are in order. (i) In principle, spin-1 resonances, as for example
heavy gluons, could also mediate this process [47, 48, 49]. Were this the
case, the production cross section would be larger: our bounds considering
only the QCD contribution would be conservative. (ii) Likewise, in the context
of CHMs, the cross section for single production of heavy tops might be also
large [50, 51].

##  3 Searches for vector-like quarks with standard decays

Figure 1: NNLO cross sections for pair production of top partners [33, 35] at
8 (left) and 13 (right) TeV of c.m.e.

Standard searches for QCD pair-produced top partners assume that these decay
only into SM particles, namely

$\text{BR}(T\rightarrow Ht)+\text{BR}(T\rightarrow Wb)+\text{BR}(T\rightarrow
Zt)=1~{}.$ (3.1)

In fact, bounds are typically given by plots showing excluded regions in the
plane of the first two branching ratios (see for example all searches
discussed in ref. [52]). Consequently, bounds for top partners with elusive
decays, for which this equality does not hold, cannot be read from current
plots. Some constraints can be still obtained by rescaling the ones given for
each channel separately. However, these disregard all signal events in which
the two heavy tops decay into different particles. Bounds obtained in this way
are hence highly underestimated. At first glance, any attempt of departing
from eq. 3.1 requires reinterpreting the data from scratch (e.g. recasting the
corresponding analyses).

Fortunately, this exercise can be greatly simplified in a counting experiment
if limits are provided for at least three different sets of branching ratios.
We will illustrate how to proceed in that case in the first SM channel to
which we pay attention, namely $T\rightarrow Zt$.

###  3.1 $T\rightarrow Zt$

Let us consider a counting experiment for $T\overline{T}\rightarrow Zt+X$,
with $X$ being any standard $T$ decay channel. Let us assume that bounds on
the QCD model-independent pair-production cross section are known for
$\text{BR}(T\rightarrow Zt)=1$. We will dub this quantity $\sigma_{1}$. Then,
it is clear that the maximum amount of allowed signal events after cuts, $N$,
can be written as $N=\sigma_{1}\,L\,\epsilon(ZtZt)$, where $L$ stands for the
integrated luminosity considered in the analysis and $\epsilon(ZtZt)$ is the
efficiency for selecting signal events in the $Zt+Zt$ channel.

Let us further assume that limits for the cases $\text{BR}(T\rightarrow
Ht)=\text{BR}(T\rightarrow Zt)=0.5$ (the doublet case) and
$2~{}\text{BR}(T\rightarrow Ht)=2~{}\text{BR}(T\rightarrow
Zt)=\text{BR}(T\rightarrow Wb)=0.5$ (the singlet case) are also known. We will
name their values by $\sigma_{D}$ and $\sigma_{S}$, respectively. Then, we get
that the following equations approximately hold:

$N\sim\sigma_{D}\,L\,\bigg{[}\frac{1}{4}\epsilon(ZtZt)+\frac{1}{2}\epsilon(ZtHt)\bigg{]}\sim\sigma_{S}\,L\,\bigg{[}\frac{1}{16}\epsilon(ZtZt)+\frac{1}{8}\epsilon(ZtHt)+\frac{1}{4}\epsilon(ZtWb)\bigg{]}~{}.$
(3.2)

We have (conservatively) neglected $\epsilon(HtHt),\epsilon(WbWb)$ as well as
$\epsilon(HtWb)$, given that the experimental search tags mainly $Zt+X$
events. From the equation above, all the relevant efficiencies can be
estimated, and therefore the number of expected signal events passing all the
analysis cuts, $N_{\text{sig}}$, for arbitrary values of the branching ratios:

$\displaystyle N_{\text{sig}}\sim\sigma\,L\,\bigg{[}\text{BR}(T\rightarrow
Zt)^{2}\epsilon(ZtZt)$ $\displaystyle+2~{}\text{BR}(T\rightarrow
Zt)\,\text{BR}(T\rightarrow Ht)~{}\epsilon(ZtHt)$
$\displaystyle+2~{}\text{BR}(T\rightarrow Zt)\,\text{BR}(T\rightarrow
Wb)~{}\epsilon(ZtWb)\bigg{]}~{}.$ (3.3)

In this expression, $\sigma$ represents the theoretical cross section for
pair-production of top partners. This is shown in fig. 1, assuming that is
driven by QCD interactions. The number of signal events can be compared with
the number of observed and predicted SM events. Altogether, this method can be
used to obtain bounds on top-like resonances decaying into SM particles even
when eq. 3.1 is not fulfilled. Conversely, if eq. 3.1 is satisfied, the
constraints obtained using this approach are slightly conservative. Let us
apply this procedure to a real example. This will allow us to check the
goodness of the method. In addition, the results obtained in this way will be
used in next parts of this article, when bounds combining different decay
modes (including exotic $T$ decays) will be given.

efficiency (%) $M=800$ GeV $M=1000$ GeV $M=1100$ GeV $Zt+Zt$ 0.5 1.0 1.3
$Zt+Ht$ 0.3 0.8 0.9 $Zt+Wb$ 0.2 0.4 0.5

Table 1: Relevant efficiencies of the analysis of ref. [36] for three
different masses of the heavy top as obtained using the procedure outlined in
the text.

We consider the latest search for $Zt+X$ of ref. [36]. This is a counting
experiment based on 36 fb-1 of collected data at 13 TeV of c.m.e. at the ATLAS
experiment. It focuses mainly on the invisible decay of the $Z$ boson. Among
the most important cuts, the analysis requires the presence of exactly one
light lepton (either electron or muon), missing transverse energy
$E_{T}^{\text{miss}}>350$ GeV, the presence of at least four hard small-radius
jets and the presence of two large-radius jets with high $p_{T}$. The search
is optimized for a heavy top with mass $M\gtrsim 1$ TeV, for which the
selection efficiency is reported to be $\sim 1\%$. (Note that this small
efficiency takes into account all SM decay modes.) Bounds at the 95 % C.L. for
the case $\text{BR}(T\rightarrow Zt)=1$ as well as for the singlet and the
doublet cases are provided in fig. 6 of that reference for masses between 600
and 1400 GeV. Thus, following the procedure outlined above, we obtain the
relevant efficiencies. These are written in tab. 1 for three representative
masses. Clearly, although the efficiencies in the $Zt+Zt$ channel are larger,
the efficiencies in the mixed channels are not negligible. The number of
predicted SM events passing all cuts is reported to be $N_{\text{SM}}\sim
6.5$, while the number of observed events is $N_{\text{O}}=7$. The maximum
number of allowed signal events can then be obtained using a simple CLs method
[53]. This method considers the statistic

$Q=\prod_{i}\frac{(s_{i}+b_{i})^{\tilde{n}_{i}}e^{-(s_{i}+b_{i})}}{b_{i}^{\tilde{n}_{i}}e^{-b_{i}}}=e^{-\sum_{i}s_{i}}\prod_{i}\left[1+\frac{s_{i}}{b_{i}}\right]^{\tilde{n}_{i}}=e^{-N_{\text{sig}}}\left[1+\frac{N_{\text{sig}}}{N_{\text{SM}}}\right]^{\tilde{N}_{\text{O}}}~{}.$
(3.4)

In the last equality we have used that the number of independent signal
regions, over which the index $i$ runs, is exactly 1 in a counting experiment.
Besides, $\tilde{N}_{\text{O}}$ represents a random Poisson distributed
variable with mean $N_{\text{sig}}+N_{\text{SM}}$ in the signal+background
hypothesis, and with mean $N_{\text{O}}$ in the background-only hypothesis.
Let us call $P_{\text{sig+SM}}(Q)$ and $P_{\text{SM}}(Q)$ the random
distribution followed by $Q$ in the first and second cases, respectively. Two
confidence estimators are defined in the CLs method:

$\text{CL}_{s+b}=1-\int_{Q_{\text{obs}}}^{\infty}P_{\text{sig+SM}}(Q)dQ,\qquad\text{CL}_{b}=1-\int_{Q_{\text{obs}}}^{\infty}P_{\text{SM}}(Q)dQ~{},$
(3.5)

where $Q_{\text{obs}}$ is the value of $Q$ when
$\tilde{N}_{\text{O}}=N_{\text{O}}$. A signal is said to be excluded at the 95
% C.L. if CL${}_{s}=\text{CL}_{s+b}/\text{CL}_{s}<0.05$. In this case, this
occurs for $N_{\text{sig}}>8$. Fixing eq. 3.1, we compute this ratio using the
TLimit class of ROOT [54] for the three values of $M$ in tab. 1 and for all
values of $\text{BR}(T\rightarrow Ht)$ and $\text{BR}(T\rightarrow Wb)$ in the
range $(0,1)$ with steps of $0.01$. $N_{\text{sig}}$ is computed using eq. 3.1
and the efficiencies quoted in the table. The excluded regions are shown in
orange in the three panels of fig. 2. The bounds obtained in (fig. 7 of) the
experimental reference are superimposed with green dashed lines. Clearly, the
agreement is excellent.

Figure 2: Left) Region of the plane of branching ratios excluded by the $Zt+X$
analysis of ref. [36] for $M=800$ GeV. The dashed green lines represent the
constraints provided by the experimental collaboration, while the orange area
is the one we obtain with the simplified method outlined in the text. Middle)
Same as Left) but for $M=1000$ GeV. Right) Same as Left) but for $M=1100$ GeV.

###  3.2 $T\rightarrow Ht$

We follow the same approach as before to obtain the efficiencies for selecting
$Ht+X$ events. We consider the experimental analysis of ref. [33]. It is based
on 20 fb-1 of collected luminosity at 8 TeV of c.m.e. at the ATLAS experiment.
Despite the bounds reported by the paper are based on a likelihood analysis of
several signal regions, the sensitivity of the search is driven mainly by one
of them. So, it can be approximated by a counting experiment. In this signal
region, at least six light jets have to be present, of which at least four
must be $b$-tagged. Again, exactly one light lepton is required. The analysis
reports the observation of 84 events, while 81 were expected from the SM
alone. Consequently, no more than 22 signal events are allowed. To validate
the approach in this case, we consider a heavy top mass of 750 GeV. This
choice is intended to combine the $Ht+X$ channel with the $Wb+X$ one, for
which ref. [33] provides bounds too. (For other masses, the regions excluded
combining the $Ht+X$ and $Wb+X$ channels are just the trivial overlapping of
those excluded by each search separately. We elaborate on this point in the
next section.) The bounds that we obtain are shown in the left panel of fig.
3. The agreement is worse than in the $Zt+X$ case, the reason being that this
analysis is not exactly a counting experiment. The results, however, are
conservative.

On top of this analysis, we consider the latest ATLAS search for $Ht+X$
described in ref. [35]. It is based on 13.2 fb-1 of collected luminosity at 13
TeV of c.m.e. We consider the channel that requires the presence of exactly
one light lepton. The search further imposes a stringent cut on the
multiplicity of jets and $b$-tagged jets, as well as on the scalar sum of the
transverse momenta of all final state objects. Boosted techniques are used to
reconstruct hadronically-decaying resonances giving rise to large-radius jets.
Contrary to the previous searches, the bounds obtained in this one rely on the
several bins of the effective mass, defined as the scalar sum of the
transverse momenta of the lepton, the selected jets and the
$E_{T}^{\text{miss}}$. Consequently, we can not proceed as before. Instead, we
consider just the reported excluded cross sections, $\sigma_{\text{excl}}$,
for $\text{BR}(T\rightarrow Ht)=1$. Thus, a heavy top of a given mass $M$ is
ruled out by this analysis if

$\sigma(pp\rightarrow T\overline{T})\times\text{BR}(T\rightarrow
Ht)^{2}>\sigma_{\text{excl}}~{}.$ (3.6)

Bounds obtained in this way can not be combined with the ones before, but they
will be superimposed to them. As a matter of fact, when this branching ratio
is close to 1, this search sets the most stringent constraint.

###  3.3 $T\rightarrow Wb$

A search for $Wb+X$ has been also carried out in ref. [33]. Among the most
important cuts, exactly one light lepton is required, as well as at least four
light jets, one of which must be $b$-tagged. Besides, one jet with $p_{T}>400$
GeV has to be present, too. This jet results typically from the hadronic decay
of one $W$ boson. According to the study, 27.6 SM events are expected to pass
all cuts, while 30 have been observed. This implies that no more than 15
signal events are allowed. The analysis does not provide bounds for the
doublet case, so we only obtain the efficiency for selecting $Wb+Wb$ events
using the simplified method explained above. The region excluded using this
approach for $M=750$ GeV is shown in the middle panel of fig. 3. It matches
perfectly the one obtained by the experimental collaboration, which means that
very few events result from the disregarded $Wb+Ht$ or $Wb+Zt$ channels.

Figure 3: Left) Region of the plane of branching ratios excluded by the $Ht+X$
analysis of ref. [33] for $M=750$ GeV. The dashed green lines represent the
constraints provided by the experimental collaboration, while the orange area
is the one we obtain with the simplified method outlined in the text. Middle)
same as Left) but for the $Wb+X$ analysis. Right) Same as Left) but for the
combination of both analyses.

Using the CLs method described at the beginning of this section, we obtain
bounds combining the signal region of the 8 TeV $Ht+X$ analysis with that of
the $Wb+X$ one (they are statistically independent). The result is shown in
the right panel of the figure. Such a combination has been also provided by
the experimental collaboration, and it is depicted by the green-dashed line.
Our results are again conservative. On top of this, we also consider the
latest search for $Wb+X$ described in ref. [34], based on 14.7 fb-1 of
collected luminosity at 13 TeV of c.m.e. In the final plots, we will also
superimpose constraints obtained from those reported for
$\text{BR}(T\rightarrow Wb)=1$.

##  4 Searches for vector-like quarks with exotic decays

###  4.1 $T\rightarrow St,S\rightarrow b\overline{b}$

Figure 4: Left) Normalized distribution of $m_{bb}^{\text{min}\Delta R}$ for
$M=600$ GeV and $m_{S}=200$ GeV in the $Ht+Wb$ (thick solid orange) and the
$St+Wb$ (thin solid green) channels. The cut on $m_{bb}^{\text{min}\Delta
R}>100$ GeV is depicted with a vertical dashed line. Right) Same as Left) but
for $M=800$ GeV and $m_{S}=300$ GeV.

The search for $Ht+X$ [33] relies on the Higgs decaying into bottom quarks.
Hence, it is also sensitive to $St+X$ events with $S$ decaying into the same
final state. In fact, if $m_{S}$ were of the size of the Higgs mass, the
sensitivity of this search for both channels would be similar. A thorough
inspection of this search reveals that the only variable that can distinguish
both channels is actually $m_{bb}^{\text{min}\Delta R}$, defined as the
invariant mass of the two $b$-tagged jets closest in $\Delta R$. It is
required to be larger than 100 GeV.

The boost factor of $S$ grows as $\sim M/m_{S}$. Therefore, the larger its
mass, the smaller the fraction of events for which the two bottom quarks
resulting from its decay are the closest ones. Therefore, it could be expected
that, for small $M$ and large $m_{S}$, the sensitivity of this search for
$St+X$ reduced with respect to $Ht+X$. However, this effect turns out to be
largely compensated by the smaller Higgs decay rate into bottom quarks and the
fact that the invariant mass of the closest $b$-tagged jets peaks at much
larger values in the $S$ case. This is made explicit in fig. 4 for two
different pairs of $M$ and $m_{S}$. The distributions are obtained out of QCD
pair-produced heavy tops generated with MadGraph v5 [55] whose SM final state
particles are subsequently decayed with Pythia v6 [56]. The basic cuts of ref.
[33], described in section 3.3, have been applied. The efficiency with respect
to the standard Higgs channel varies from $\sim 1.1$, for
$(M,m_{S})=(600,200)$ GeV, to $\sim 1.2$, for $(M,m_{S})=(1500,200)$ GeV.

Similarly, the efficiency of the search for $Zt+X$ [36] examined in section
3.1 for the $Zt+St$ channel matches also the one for $Zt+Ht$ when $S$ decays
to bottom quarks. As a result, we can safely add $\text{BR}(T\rightarrow
St)+\text{BR}(T\rightarrow Ht)$ into a single variable that, abusing notation,
we will still write as $\text{BR}(T\rightarrow Ht)$. Bounds considering this
new channel are then trivial to obtain. Hereafter, we shall only write $St$
explicitly for the case in which $S$ scapes detectors, to be discussed in the
next section.

###  4.2 $T\rightarrow St,S\rightarrow E_{T}^{\text{miss}}$

The search for $Zt+X$ focuses on the invisible decay of the $Z$ boson.
Therefore it is also sensitive to $St+X$ channels whenever $S$ is stable at
detector scales. This analysis puts a stringent cut on the amount of
$E_{T}^{\text{miss}}$. The efficiency depends notably on $M/m_{S}$, and hence,
contrary to the previous case, $\text{BR}(T\rightarrow St)$ can not just be
added to $\text{BR}(T\rightarrow Zt)$. Instead, we need to recast the whole
analysis. For this matter, we use homemade routines based on a combination of
MadAnalysis v5 [57, 58], ROOT v5 [54] and Fastjet v3 [59]. QCD pair-produced
heavy tops decaying into $St+Zt,St+Ht$ and $St+Wb$ are generated with MadGraph
v5 [55] and Pythia v6 [56]. Detector effects are disregarded, although an
efficiency of 0.85 (0.9) for selecting electrons (muons) has been simulated.
Among other variables, this analysis relies on
$H_{T,\text{sig}}^{\text{miss}}=(H_{T}^{\text{miss}}-100~{}\text{GeV})/\sigma_{H_{T}^{\text{miss}}}$,
where $H_{T}^{\text{miss}}$ stands for the $p_{T}$ of the vectorial sum of the
jet and lepton momenta. Likewise, $\sigma_{H_{T}^{\text{miss}}}$ represents
the resolution on $H_{T}^{\text{miss}}$. We take a conservative value of 7 %.
In order to validate the implementation, we recalculate the efficiencies for
the $Zt+Zt$ channel. We find very good agreement with those computed in
section 3.1, some of which have been shown in tab. 1.

Figure 5: Left) excluded region in the stop-neutralino mass plane by the CMS
search of ref. [37]. The area enclosed by the solid (dashed) orange (green)
line is obtained with the recast analysis (by the experimental collaboration).
Right) Bound on $\text{BR}(T\rightarrow St)$ for $m_{S}=100$ (dashed green),
$200$ (dashed orange), $300$ (solid orange) and $400$ (solid green) GeV.

Then, we obtain the efficiencies for $St+X$ for different values of $M$ and
$m_{S}$. A representative set of these efficiencies is collected in the
different tables of appendix A, together with those computed in the previous
sections.

When both top partners decay into $St$, the resulting signal is identical to
that explored by searches for stops decaying into tops and neutralinos. We
recast one of the latest analyses carried out by CMS [37] at 13 TeV of c.m.e.
in the fully-hadronic channel (we also considered the homologous search by the
ATLAS Collaboration [60], but the validation was far less convincing). The
large top decay rate into jets together with the current performance of
boosted techniques make this channel optimal for probing large values of $M$
and $m_{S}$.

The study relies on 60 statistically independent signal regions. We find that
the sensitivity is driven mainly by two of them. Roughly, the first one
requires $E_{T}^{\text{miss}}$ to be between $250$ and $350$ GeV, as well as
$n_{t}=0$ and $n_{W}=1$. $n_{t}$ (resp. $n_{W}$) is the number of fat jets
with mass within 110 and 120 (resp. 60 and 110) GeV and $p_{T}>$ 400 (resp.
200) GeV. Fat jets are clustered out of hadronic traces with the anti-$k_{t}$
algorithm [61] with $R=0.8$ and subsequently re-clustered with the Cambridge-
Aachen algorithm [62] with the same $R$. The analysis reports the observation
of 14 events, while only 9.6 are expected from the SM alone (see tab. 2).
Consequently, any model giving more than 13 events in this signal region is
excluded at the 95 % C.L. The second signal region, instead, requires
$E_{T}^{\text{miss}}>500$ GeV, $n_{t}\geq 1$ and $n_{W}\geq 1$. The number of
observed and expected events are 1 and 0.16, respectively. This sets an upper
limit on the number of allowed signal events after cuts of 5. In both cases,
at least five $R=0.4$ anti-$k_{t}$ jets and at least two $b$-tagged jets have
to be present.

$Ht+X$ $Wb+X$ $Zt+X$ CMS c.m.e. [TeV] 8 8 13 13 $HtHt$, $HtWb$, $WbWb$ $ZtZt$,
$ZtHt$, $StSt$ Channels $HtZt$ $ZtWb$, $StHt$, $StZt$, $StWb$ $N_{b}$ 81
$27.6$ $6.5$ 9.6, 0.16 $N_{d}$ 84 30 $7$ 14, 1 $N_{s}$ 22 15 $8$ 13, 5

Table 2: Experimental analyses considered in this work and channels for which
efficiencies have been estimated. The number of observed and expected SM
events in each signal region are also shown. Any signal leading to more than
$N_{s}$ events in at least one analysis is excluded, although stronger bounds
are obtained if the different analyses are combined into a shape analysis. See
the text for details.

Figure 6: Excluded regions in the space of branching ratios for $M=700$ GeV
and $m_{S}=100$ GeV. The gray area results from a combined statistical
analysis of all signal regions of tab. 2. In the green area, the $St+X$ events
are disregarded. The red region is excluded by the analysis of ref. [35].
Likewise, the blue region is excluded by the analysis of ref. [34].

For the validation we produce again pairs of stops with MadGraph v5 that are
subsequently decayed only into tops and neutralinos with Pythia v6. We then
obtain the number of survivor signal events in each of the 60 signal regions
considered in the experimental analysis. We combine all them into a single
statistics which is compared with the provided number of expected and observed
events via a CLs analysis. The excluded region in the plane of stop and
neutralino masses ($m_{\tilde{t}}-m_{\tilde{\chi}^{0}}$) is depicted by the
orange region in the left panel of fig. 5. The limits reported by the
experimental collaboration are also shown for comparison, the agreement being
apparent. We can therefore obtain the efficiency of this analysis for $St+St$
events for different values of $m_{S}$. As the rest of efficiencies, these are
reported in appendix A. In the right panel of fig 5 we plot the maximum value
of $\text{BR}(T\rightarrow St)$ that is allowed by this analysis when it is
applied to the $St+St$ channel alone taken into account only the two main
important signal regions. Note that values of $M$ as large as 1100 GeV are
excluded for $\text{BR}(T\rightarrow St)=1$. Much lighter resonances are still
allowed if the branching ratio is smaller. In that case, however, decays into
SM particles are also present. Consequently both standard and exotic $T$
decays have to be considered at once. We address this point in the next
section.

##  5 Final results

At this stage, the number of signal events passing the cuts of the analyses
reported in tab. 2 can be computed for arbitrary branching ratios of a heavy
top $T$ into $Zt,Ht,Wb$ ot $St$. We recall that $St$ stands for an invisible
$S$, while the $\text{BR}(T\rightarrow St,S\rightarrow b\overline{b})$ is
included in $\text{BR}(T\rightarrow Ht)$. The number of expected and observed
events for each analysis are also written in the table. Given this
information, bounds combining all these searches can be obtained by using the
CLs method as described in section 3.1. These are shown in figs. 6, 7, 8 and 9
for different values of $M$ and different choices of $\text{BR}(T\rightarrow
St)$. Given the small dependence of searches for $St+St$ on the value of
$m_{S}$ (see the right panel of fig. 5), we have fixed $m_{S}=100$ GeV in all
plots. We have also assumed

$\text{BR}(T\rightarrow Ht)+\text{BR}(T\rightarrow Wb)+\text{BR}(T\rightarrow
Zt)+\text{BR}(T\rightarrow St)=1~{}.$ (5.1)

The different colors have the following meanings. The gray area enclosed by
the solid black line shows the region excluded after combining all the
(statistically independent) signal regions presented in tab. 2. In the green
one, enclosed by the dashed green line, we do the same, but neglecting the
events resulting from $T\rightarrow St$. This region gives an idea of how
constraining current searches are if $T$ decays also into some elusive channel
(for example into very soft final state particles). In other words, it
reflects departures from eq. 3.1 by an amount of $\text{BR}(T\rightarrow St)$.
Clearly, in the limit in which $\text{BR}(T\rightarrow St)=0$, this region
coincides with the gray one. In the limit $\text{BR}(T\rightarrow St)=1$, this
region is empty. The red region enclosed by the dashed red line represents the
area excluded by searches for $Ht+X$ at 13 TeV of c.m.e., as explained in
section 3.2. Finally, the blue region enclosed by the dashed blue line
represents the area excluded by searches for $Wb+X$ at 13 TeV of c.m.e., as
explained in section 3.3. Let us discuss these results case by case.

Figure 7: Same as fig. 6 but for $M=800$ GeV.

$\mathbf{M=700}$ GeV: All combinations of branching ratios are excluded,
independently on the value of $\text{BR}(T\rightarrow St)$. Note that we do
not show results for values of this rate larger than 0.3. The reason is that,
according to the plot in the right panel of fig. 5, such scenarios are already
excluded by searches for stops decaying into neutralinos.

$\mathbf{M=800}$ GeV: In this case, a small region of the plane when
$\text{BR}(T\rightarrow St)<0.1$ is still allowed by the data we are
considering. We stress, however, that this region disappears when all signal
regions in the search for $Ht+X$ (see ref. [35]) are considered. It is also
worth to emphasize the power of combining, in a single statistical analysis,
the data from heavy tops with standard decays with those for top partners with
exotic decays. Indeed, let us focus on the case $\text{BR}(T\rightarrow
St)=0.2$ (left lower panel). The region excluded by searches for top partners
decaying into $Zt,Ht$ or $Wb$ is depicted by the green area. Clearly, this is
far from being the whole parameter space. On the other hand,
$\text{BR}(T\rightarrow St)=0.2$ is perfectly allowed by searches for $St+St$
(right panel of fig. 5). The combination of the various channels, however,
gives much more stringent constraints than the mere superposition of those
coming from the different searches alone.

Figure 8: Same as fig. 6 but for $M=900$ GeV.

$\mathbf{M=900}$ and 1000 GeV: In both cases, there are sizable regions of the
parameter space not excluded by current searches for $\text{BR}(T\rightarrow
St)<0.2$. For larger values of this branching ratio, however, the combination
of standard and exotic decays forbids again the existence of top partners,
independently of their branching ratio into SM channels.

Figure 9: Same as fig. 6 but for $M=1000$ GeV.

$\mathbf{M>1000}$ GeV: independently of the value of $\text{BR}(T\rightarrow
St)$, searches for $pp\rightarrow T\overline{T}\rightarrow t\overline{t}SS$
are no longer constraining (see fig. 5). Likewise, standard searches loose
sensitivity for $M>1100$ GeV. For different values of $m_{S}$ and/or different
values of the branching ratios, we refer the reader to appendix B, in which we
explain how to obtain experimental bounds straightforwardly using a simple
Python script.

##  6 Conclusions

Heavy tops, $T$, which are predicted by different extensions of the Standard
Model (SM), can decay not only into SM final states such as $Zt,Ht$ or $Wb$,
but also into new particles. Among others, new scalar singlets, $S$, are good
candidates. These appear, for example, in non-minimal realizations of the
composite Higgs model paradigm. In light of this, we have investigated the LHC
signals of top partners with sizable decay rate also into $St$. We have
focused on the pair-production mode, because it depends only on the model-
independent QCD interactions. We have considered ATLAS and CMS experimental
analyses performed at both 8 and 13 TeV of center of mass energy. These
include searches for $Ht+X$ (which are also sensitive to $St+X$ final states
when $S$ decays into bottom quarks) and $Wb+X$, as well as analyses of $Zt+X$
and supersymmetry searches for pair produced stops with neutralino decays
(both of which are also sensitive to $St+X$ events if $S$ is a stable
particle). We have obtained the expected number of signal events passing the
cuts of these analyses for arbitrary combinations of branching ratios and for
different masses of $T$ and $S$. Then, we have obtained bounds resulting from
a combined statistical study of all signal regions. Besides, these two
procedures have been implemented in a simple code that can be found in
http://github.com/mikaelchala/vlqlimits. Details about the mode of use are
provided in appendix B. On the physics side, we have found that heavy tops
with masses below $\sim 900$ GeV are excluded by direct searches for any value
of $\text{BR}(T\rightarrow St)$, while even larger masses are forbidden for
larger values of this branching ratio. All in all, we can conclude that bounds
on heavy tops are not necessarily weaker if these decay also into new
particles.

##  Acknowledgments

I am grateful to J. P. Araque and N. Castro for the useful discussions on the
experimental side of this work. I would like to thank also G. Nardini and J.
Santiago for useful comments on the manuscript. This work is partially
supported by the Spanish MINECO under grant FPA2014-54459-P and by the Severo
Ochoa Excellence Program under grant SEV-2014-0398.

##  Appendix A Tables

The efficiencies of the experimental analyses (see tab. 2) for the different
final states considered in this work for several values of $M$ are shown in
rows $2-12$ of tab. 3 ($m_{S}=100$ GeV), tab. 4 ($m_{S}=200$ GeV), tab. 5
($m_{S}=300$ GeV) and tab. 6 ($m_{S}=400$ GeV).

$M$ [GeV] 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 $Zt+Zt$ 0.004
0.004 0.006 0.009 0.010 0.013 0.015 0.016 0.017 0.017 0.017 $Zt+Ht$ 0.002
0.002 0.003 0.005 0.008 0.009 0.012 0.014 0.015 0.015 0.015 $Zt+Wb$ 0.001
0.000 0.002 0.002 0.004 0.005 0.006 0.007 0.008 0.008 0.008 $Ht+Ht$ 0.056
0.083 0.104 0.139 0.167 0.185 0.185 0.185 0.185 0.185 0.185 $Ht+Zt$ 0.039
0.068 0.105 0.151 0.179 0.193 0.193 0.193 0.193 0.193 0.193 $Ht+Wb$ 0.026
0.053 0.066 0.149 0.208 0.257 0.257 0.257 0.257 0.257 0.257 $Wb+Wb$ 0.021
0.025 0.027 0.031 0.031 0.031 0.031 0.031 0.031 0.031 0.031 $St+Ht$ 0.001
0.003 0.005 0.007 0.009 0.010 0.010 0.010 0.010 0.010 0.010 $St+Zt$ 0.002
0.004 0.008 0.011 0.014 0.017 0.019 0.020 0.020 0.021 0.021 $St+Wb$ 0.001
0.002 0.004 0.006 0.009 0.010 0.010 0.011 0.010 0.011 0.011 $St+St(1)$ 0.009
0.007 0.007 0.006 0.006 0.006 0.007 0.006 0.006 0.006 0.006 $St+St(2)$ 0.003
0.006 0.008 0.009 0.010 0.010 0.010 0.009 0.009 0.009 0.009 $Wb+Wb$ [pb] 0.170
0.130 0.070 0.048 0.025 0.019 0.017 0.014 0.013 0.013 0.013 $Ht+Ht$ [pb] 0.106
0.081 0.055 0.031 0.023 0.016 0.012 0.010 0.009 0.009 0.009

Table 3: Efficiencies of the experimental analyses (see tab. 2) for the
different final states considered in this work for $m_{S}=100$ GeV and several
values of $M$ (rows $2-12$). The last two rows show the upper limit on the
pair production of $T$ quarks with subsequent decay into the indicated final
state at 13 TeV of c.m.e. obtained from refs. [34] and [35], respectively.

$M$ [GeV] 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 $Zt+Zt$ 0.004
0.004 0.006 0.009 0.010 0.013 0.015 0.016 0.017 0.017 0.017 $Zt+Ht$ 0.002
0.002 0.003 0.005 0.008 0.009 0.012 0.014 0.015 0.015 0.015 $Zt+Wb$ 0.001
0.000 0.002 0.002 0.004 0.005 0.006 0.007 0.008 0.008 0.008 $Ht+Ht$ 0.056
0.083 0.104 0.139 0.167 0.185 0.185 0.185 0.185 0.185 0.185 $Ht+Zt$ 0.039
0.068 0.105 0.151 0.179 0.193 0.193 0.193 0.193 0.193 0.193 $Ht+Wb$ 0.026
0.053 0.066 0.149 0.208 0.257 0.257 0.257 0.257 0.257 0.257 $Wb+Wb$ 0.021
0.025 0.027 0.031 0.031 0.031 0.031 0.031 0.031 0.031 0.031 $St+Ht$ 0.001
0.003 0.006 0.008 0.011 0.012 0.012 0.014 0.013 0.014 0.014 $St+Zt$ 0.002
0.005 0.009 0.011 0.017 0.020 0.022 0.023 0.023 0.024 0.024 $St+Wb$ 0.001
0.002 0.005 0.007 0.010 0.010 0.013 0.013 0.013 0.012 0.012 $St+St(1)$ 0.009
0.008 0.007 0.006 0.005 0.006 0.006 0.005 0.006 0.005 0.005 $St+St(2)$ 0.004
0.006 0.009 0.011 0.012 0.012 0.012 0.013 0.013 0.012 0.012 $Wb+Wb$ [pb] 0.170
0.130 0.070 0.048 0.025 0.019 0.017 0.014 0.013 0.013 0.013 $Ht+Ht$ [pb] 0.106
0.081 0.055 0.031 0.023 0.016 0.012 0.010 0.009 0.009 0.009

Table 4: Same as tab. 3 but for $m_{S}=200$ GeV.

$M$ [GeV] 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 $Zt+Zt$ 0.004
0.004 0.006 0.009 0.010 0.013 0.015 0.016 0.017 0.017 0.017 $Zt+Ht$ 0.002
0.002 0.003 0.005 0.008 0.009 0.012 0.014 0.015 0.015 0.015 $Zt+Wb$ 0.001
0.000 0.002 0.002 0.004 0.005 0.006 0.007 0.008 0.008 0.008 $Ht+Ht$ 0.056
0.083 0.104 0.139 0.167 0.185 0.185 0.185 0.185 0.185 0.185 $Ht+Zt$ 0.039
0.068 0.105 0.151 0.179 0.193 0.193 0.193 0.193 0.193 0.193 $Ht+Wb$ 0.026
0.053 0.066 0.149 0.208 0.257 0.257 0.257 0.257 0.257 0.257 $Wb+Wb$ 0.021
0.025 0.027 0.031 0.031 0.031 0.031 0.031 0.031 0.031 0.031 $St+Zt$ 0.002
0.004 0.007 0.009 0.015 0.015 0.015 0.017 0.016 0.015 0.015 $St+Ht$ 0.003
0.005 0.010 0.016 0.019 0.023 0.024 0.026 0.025 0.026 0.026 $St+Wb$ 0.001
0.003 0.006 0.008 0.011 0.013 0.014 0.015 0.015 0.017 0.017 $St+St(1)$ 0.009
0.008 0.006 0.005 0.005 0.005 0.005 0.004 0.004 0.004 0.004 $St+St(2)$ 0.003
0.007 0.010 0.012 0.014 0.015 0.014 0.015 0.015 0.014 0.014 $Wb+Wb$ [pb] 0.170
0.130 0.070 0.048 0.025 0.019 0.017 0.014 0.013 0.013 0.013 $Ht+Ht$ [pb] 0.106
0.081 0.055 0.031 0.023 0.016 0.012 0.010 0.009 0.009 0.009

Table 5: Same as tab. 3 but for $m_{S}=300$ GeV.

$M$ [GeV] 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 $Zt+Zt$ 0.004
0.004 0.006 0.009 0.010 0.013 0.015 0.016 0.017 0.017 0.017 $Zt+Ht$ 0.002
0.002 0.003 0.005 0.008 0.009 0.012 0.014 0.015 0.015 0.015 $Zt+Wb$ 0.001
0.000 0.002 0.002 0.004 0.005 0.006 0.007 0.008 0.008 0.008 $Ht+Ht$ 0.056
0.083 0.104 0.139 0.167 0.185 0.185 0.185 0.185 0.185 0.185 $Ht+Zt$ 0.039
0.068 0.105 0.151 0.179 0.193 0.193 0.193 0.193 0.193 0.193 $Ht+Wb$ 0.026
0.053 0.066 0.149 0.208 0.257 0.257 0.257 0.257 0.257 0.257 $Wb+Wb$ 0.021
0.025 0.027 0.031 0.031 0.031 0.031 0.031 0.031 0.031 0.031 $St+Zt$ 0.002
0.006 0.008 0.012 0.016 0.018 0.019 0.019 0.021 0.021 0.021 $St+Ht$ 0.004
0.006 0.012 0.018 0.021 0.023 0.028 0.030 0.029 0.030 0.030 $St+Wb$ 0.002
0.003 0.007 0.010 0.012 0.015 0.016 0.018 0.019 0.018 0.018 $St+St(1)$ 0.008
0.007 0.006 0.005 0.005 0.004 0.004 0.004 0.004 0.004 0.004 $St+St(2)$ 0.003
0.007 0.011 0.014 0.015 0.016 0.017 0.018 0.017 0.015 0.015 $Ht+Ht$ [pb] 0.170
0.130 0.070 0.048 0.025 0.019 0.017 0.014 0.013 0.013 0.013 $Wb+Wb$ [pb] 0.106
0.081 0.055 0.031 0.023 0.016 0.012 0.010 0.009 0.009 0.009

Table 6: Same as tab. 3 but for $m_{S}=400$ GeV.

##  Appendix B Software

A small Python script can be found in http://github.com/mikaelchala/vlqlimits.
The usage is extremely simple. First, one needs to edit the file input.dat
(the name can be of course changed). It reads:

    % THIS IS A COMMENT
    % FORMAT: M, MS, BR(T->WB), BR(T->HT)+BR(T->ST, S->BB), BR(T->ZT), BR(T->ST, S -> INV)
      735.    210.    0.5     0.1     0.1     0.3

The non-commented line contains an example point. More points, corresponding
to other possible quarks present in the spectrum, can be added (interference
effects are disregarded). We recall that the branching ratios do not have to
add to 1. This would reflect the existence of further decays for which there
is no experimental information (for example, decays into very soft particles,
which are hard to detect). Given this, one can compute whether the set of
heavy tops indicated in the file is excluded at the 95 % C.L. or not. This can
be done in two ways:

  * •

./code_indp.py input.dat: returns 1 if the number of predicted signal events
in at least one signal region is above that quoted in tab. 2, or if the cross
sections in $Wb+Wb$ or $Ht+Ht$ exceed the ones reported by the corresponding
analyses at 13 TeV of c.m.e. [34, 35]. It returns 0 otherwise.

  * •

./code_comb.py input.dat: returns 1 if the CLs computed considering all signal
regions in tab. 2 is below $0.05$. Or if the cross sections in $Wb+Wb$ or
$Ht+Ht$ exceed the ones reported by the corresponding analyses at 13 TeV of
c.m.e. [34, 35]. It returns 0 otherwise. This code can be only run if PyROOT
is properly installed.

In both cases, the NumPy and the SciPy libraries are mandatory. The cross
sections and efficiencies for masses not considered in the Monte Carlo
simulations are linearly interpolated from those included in the hidden file
.tables. More sophisticated tools, that include also heavy fermions with
different quantum numbers, are also publicly available (see for example XQCAT
[32, 63]). However, exotic decays and signals at 13 TeV are generally not
considered. In this respect, this code complements previous works.

## References

  * [1] J. A. Aguilar-Saavedra, _Pair production of heavy Q = 2/3 singlets at LHC_ , _Phys. Lett._ B625 (2005) 234–244, [hep-ph/0506187].
  * [2] J. A. Aguilar-Saavedra, _Identifying top partners at LHC_ , _JHEP_ 11 (2009) 030, [0907.3155].
  * [3] J. A. Aguilar-Saavedra, R. Benbrik, S. Heinemeyer and M. Pérez-Victoria, _Handbook of vectorlike quarks: Mixing and single production_ , _Phys. Rev._ D88 (2013) 094010, [1306.0572].
  * [4] D. B. Kaplan and H. Georgi, _SU(2) x U(1) Breaking by Vacuum Misalignment_ , _Phys. Lett._ B136 (1984) 183–186.
  * [5] D. B. Kaplan, H. Georgi and S. Dimopoulos, _Composite Higgs Scalars_ , _Phys. Lett._ B136 (1984) 187–190.
  * [6] S. Dimopoulos and J. Preskill, _Massless Composites With Massive Constituents_ , _Nucl. Phys._ B199 (1982) 206–222.
  * [7] D. B. Kaplan, _Flavor at SSC energies: A New mechanism for dynamically generated fermion masses_ , _Nucl. Phys._ B365 (1991) 259–278.
  * [8] K. Agashe, R. Contino and A. Pomarol, _The Minimal composite Higgs model_ , _Nucl. Phys._ B719 (2005) 165–187, [hep-ph/0412089].
  * [9] B. Gripaios, A. Pomarol, F. Riva and J. Serra, _Beyond the Minimal Composite Higgs Model_ , _JHEP_ 04 (2009) 070, [0902.1483].
  * [10] M. Chala, G. Nardini and I. Sobolev, _Unified explanation for dark matter and electroweak baryogenesis with direct detection and gravitational wave signatures_ , _Phys. Rev._ D94 (2016) 055006, [1605.08663].
  * [11] M. Chala, _$h\rightarrow\gamma\gamma$ excess and Dark Matter from Composite Higgs Models_, _JHEP_ 01 (2013) 122, [1210.6208].
  * [12] V. Sanz and J. Setford, _Composite Higgses with seesaw EWSB_ , _JHEP_ 12 (2015) 154, [1508.06133].
  * [13] L. Vecchi, _The Natural Composite Higgs_ , 1304.4579.
  * [14] M. Frigerio, A. Pomarol, F. Riva and A. Urbano, _Composite Scalar Dark Matter_ , _JHEP_ 07 (2012) 015, [1204.2808].
  * [15] J. Mrazek, A. Pomarol, R. Rattazzi, M. Redi, J. Serra and A. Wulzer, _The Other Natural Two Higgs Doublet Model_ , _Nucl. Phys._ B853 (2011) 1–48, [1105.5403].
  * [16] N. Fonseca, R. Zukanovich Funchal, A. Lessa and L. Lopez-Honorez, _Dark Matter Constraints on Composite Higgs Models_ , _JHEP_ 06 (2015) 154, [1501.05957].
  * [17] Y. Wu, B. Zhang, T. Ma and G. Cacciapaglia, _Composite Dark Matter and Higgs_ , 1703.06903.
  * [18] G. Ballesteros, A. Carmona and M. Chala, _Exceptional Composite Dark Matter_ , 1704.07388.
  * [19] J. R. Espinosa, B. Gripaios, T. Konstandin and F. Riva, _Electroweak Baryogenesis in Non-minimal Composite Higgs Models_ , _JCAP_ 1201 (2012) 012, [1110.2876].
  * [20] F. Caracciolo, A. Parolini and M. Serone, _UV Completions of Composite Higgs Models with Partial Compositeness_ , _JHEP_ 02 (2013) 066, [1211.7290].
  * [21] J. Barnard, T. Gherghetta and T. S. Ray, _UV descriptions of composite Higgs models without elementary scalars_ , _JHEP_ 02 (2014) 002, [1311.6562].
  * [22] G. Ferretti, _UV Completions of Partial Compositeness: The Case for a SU(4) Gauge Group_ , _JHEP_ 06 (2014) 142, [1404.7137].
  * [23] L. Vecchi, _A dangerous irrelevant UV-completion of the composite Higgs_ , _JHEP_ 02 (2017) 094, [1506.00623].
  * [24] T. Ma and G. Cacciapaglia, _Fundamental Composite 2HDM: SU(N) with 4 flavours_ , _JHEP_ 03 (2016) 211, [1508.07014].
  * [25] G. Ferretti, _Gauge theories of Partial Compositeness: Scenarios for Run-II of the LHC_ , _JHEP_ 06 (2016) 107, [1604.06467].
  * [26] J. Serra, _Beyond the Minimal Top Partner Decay_ , _JHEP_ 09 (2015) 176, [1506.05110].
  * [27] A. Anandakrishnan, J. H. Collins, M. Farina, E. Kuflik and M. Perelstein, _Odd Top Partners at the LHC_ , _Phys. Rev._ D93 (2016) 075009, [1506.05130].
  * [28] G. Cacciapaglia, H. Cai, A. Deandrea, T. Flacke, S. J. Lee and A. Parolini, _Composite scalars at the LHC: the Higgs, the Sextet and the Octet_ , _JHEP_ 11 (2015) 201, [1507.02283].
  * [29] J. Fan, S. M. Koushiappas and G. Landsberg, _Pseudoscalar Portal Dark Matter and New Signatures of Vector-like Fermions_ , _JHEP_ 01 (2016) 111, [1507.06993].
  * [30] S. Banerjee, D. Barducci, G. Bélanger and C. Delaunay, _Implications of a High-Mass Diphoton Resonance for Heavy Quark Searches_ , _JHEP_ 11 (2016) 154, [1606.09013].
  * [31] C. Niehoff, P. Stangl and D. M. Straub, _Electroweak symmetry breaking and collider signatures in the next-to-minimal composite Higgs model_ , 1611.09356.
  * [32] D. Barducci, A. Belyaev, M. Buchkremer, G. Cacciapaglia, A. Deandrea, S. De Curtis et al., _Framework for Model Independent Analyses of Multiple Extra Quark Scenarios_ , _JHEP_ 12 (2014) 080, [1405.0737].
  * [33] ATLAS collaboration, G. Aad et al., _Search for production of vector-like quark pairs and of four top quarks in the lepton-plus-jets final state in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector_, _JHEP_ 08 (2015) 105, [1505.04306].
  * [34] ATLAS Collaboration collaboration, _Search for pair production of heavy vector-like quarks decaying to high- $p_{T}$ $W$ bosons and b quarks in the lepton-plus-jets final state in pp collisions at $\sqrt{s}$=13 TeV with the ATLAS detector_, Tech. Rep. ATLAS-CONF-2016-102, CERN, Geneva, Sep, 2016\.
  * [35] ATLAS Collaboration collaboration, _Search for new phenomena in $t\bar{t}$ final states with additional heavy-flavour jets in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, Tech. Rep. ATLAS-CONF-2016-104, CERN, Geneva, Sep, 2016.
  * [36] ATLAS Collaboration collaboration, _Search for pair production of vector-like top quarks in events with one lepton and an invisibly decaying Z boson in $\sqrt{s}=13$ TeV pp collisions at the ATLAS detector_, Tech. Rep. ATLAS-CONF-2017-015, CERN, Geneva, Mar, 2017.
  * [37] CMS collaboration, C. Collaboration, _Search for direct top squark pair production in the fully hadronic final state in proton-proton collisions at sqrt(s) = 13 TeV corresponding to an integrated luminosity of 12.9/fb_ , .
  * [38] R. Contino, L. Da Rold and A. Pomarol, _Light custodians in natural composite Higgs models_ , _Phys. Rev._ D75 (2007) 055014, [hep-ph/0612048].
  * [39] O. Matsedonskyi, G. Panico and A. Wulzer, _Light Top Partners for a Light Composite Higgs_ , _JHEP_ 01 (2013) 164, [1204.6333].
  * [40] M. Redi and A. Tesi, _Implications of a light Higgs in composite models_ , _JHEP_ 1210 (2012) 166, [1205.0232].
  * [41] D. Marzocca, M. Serone and J. Shu, _General Composite Higgs Models_ , _JHEP_ 08 (2012) 013, [1205.0770].
  * [42] A. Pomarol and F. Riva, _The composite Higgs and light resonance connection_ , _JHEP_ 1208 (2012) 135, [1205.6434].
  * [43] G. Panico, M. Redi, A. Tesi and A. Wulzer, _On the Tuning and the Mass of the Composite Higgs_ , _JHEP_ 03 (2013) 051, [1210.7114].
  * [44] M. Chala, G. Durieux, C. Grojean, L. de Lima and O. Matsedonskyi, _Minimally extended SILH_ , 1703.10624.
  * [45] ATLAS, CMS collaboration, G. Aad et al., _Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV_, _JHEP_ 08 (2016) 045, [1606.02266].
  * [46] D. Ghosh, M. Salvarezza and F. Senia, _Extending the Analysis of Electroweak Precision Constraints in Composite Higgs Models_ , _Nucl. Phys._ B914 (2017) 346–387, [1511.08235].
  * [47] M. Chala, J. Juknevich, G. Perez and J. Santiago, _The Elusive Gluon_ , _JHEP_ 01 (2015) 092, [1411.1771].
  * [48] J. P. Araque, N. F. Castro and J. Santiago, _Interpretation of Vector-like Quark Searches: Heavy Gluons in Composite Higgs Models_ , _JHEP_ 11 (2015) 120, [1507.05628].
  * [49] A. Azatov, D. Chowdhury, D. Ghosh and T. S. Ray, _Same sign di-lepton candles of the composite gluons_ , _JHEP_ 08 (2015) 140, [1505.01506].
  * [50] O. Matsedonskyi, G. Panico and A. Wulzer, _On the Interpretation of Top Partners Searches_ , _JHEP_ 12 (2014) 097, [1409.0100].
  * [51] O. Matsedonskyi, G. Panico and A. Wulzer, _Top Partners Searches and Composite Higgs Models_ , _JHEP_ 04 (2016) 003, [1512.04356].
  * [52] ATLAS collaboration, J. P. Araque, _Overview of the vector-like quark searches with the LHC data collected by the ATLAS detector_ , in _9th International Workshop on Top Quark Physics (TOP 2016) Olomouc, Czech Republic, September 19-23, 2016_ , 2016.  1611.09056.
  * [53] A. L. Read, _Presentation of search results: The CL(s) technique_ , _J. Phys._ G28 (2002) 2693–2704.
  * [54] R. Brun and F. Rademakers, _ROOT: An object oriented data analysis framework_ , _Nucl. Instrum. Meth._ A389 (1997) 81–86.
  * [55] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * [56] T. Sjostrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026, [hep-ph/0603175].
  * [57] E. Conte, B. Fuks and G. Serret, _MadAnalysis 5, A User-Friendly Framework for Collider Phenomenology_ , _Comput. Phys. Commun._ 184 (2013) 222–256, [1206.1599].
  * [58] E. Conte, B. Dumont, B. Fuks and C. Wymant, _Designing and recasting LHC analyses with MadAnalysis 5_ , _Eur. Phys. J._ C74 (2014) 3103, [1405.3982].
  * [59] M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * [60] ATLAS collaboration, T. A. collaboration, _Search for the Supersymmetric Partner of the Top Quark in the Jets+Emiss Final State at sqrt(s) = 13 TeV_ , .
  * [61] M. Cacciari, G. P. Salam and G. Soyez, _The Anti-k(t) jet clustering algorithm_ , _JHEP_ 04 (2008) 063, [0802.1189].
  * [62] Y. L. Dokshitzer, G. D. Leder, S. Moretti and B. R. Webber, _Better jet clustering algorithms_ , _JHEP_ 08 (1997) 001, [hep-ph/9707323].
  * [63] D. Barducci, A. Belyaev, M. Buchkremer, J. Marrouche, S. Moretti and L. Panizzi, _XQCAT: eXtra Quark Combined Analysis Tool_ , _Comput. Phys. Commun._ 197 (2015) 263–275, [1409.3116].


# Tackling the phase problem of diffraction for retrieval of photonic
structures formed in nanocomposite materials

Martin Fally  Faculty of Physics, University of Vienna, 1090 Wien, Austria
Yasuo Tomita  Department of Engineering Science, University of Electro-
Communications, 1-5-1 Chofugaoka, Chofu, Tokyo 182, Japan  Antonio Fimia
Department of Material Science, Optics and Electronic Technology, University
Miguel Hernández, Elche, Alicante, Spain  Roque F. Madrigal  Department of
Material Science, Optics and Electronic Technology, University Miguel
Hernández, Elche, Alicante, Spain  Jinxin Guo  Present address: Institute of
Information Photonics Technology, Faculty of Science, Beijing University of
Technology, Beijing 100124, China  Department of Engineering Science,
University of Electro-Communications, 1-5-1 Chofugaoka, Chofu, Tokyo 182,
Japan  Joachim Kohlbrecher  Laboratory for Neutron Scattering, ETH Zurich &
Paul Scherrer Institut, 5232 Villigen PSI, Switzerland  Jürgen Klepp
Corresponding author<EMAIL_ADDRESS>Faculty of Physics,
University of Vienna, 1090 Wien, Austria

###### Abstract

We experimentally demonstrate how to solve the phase problem of diffraction
using multi-wave interference with standard diffraction experimental setups
without the need for taking any auxiliary data. In particular, we show that
the phases of the Fourier components of a periodic structure can be fully
recovered by deliberately choosing a probe wavelength of the diffracting
radiation much smaller than the lattice constant. In the course of the
demonstration, we accurately determine the refractive index profile of
nanocomposite phase gratings by light and neutron diffraction measurements.

Diffraction experiments are usually a method of choice to determine the
internal structure of bulk materials. Instead of a crystal, let us consider
here the simple case of a planar one-dimensional phase grating, which is
characterized by the spatial profile of the refractive index

$\tilde{n}(x)=\sum\limits_{s=-\infty}\limits^{+\infty}\tilde{n}_{s}e^{\imath
sGx}.$ (1)

Here, $G$ is the spatial frequency, $|\tilde{n}_{s}|=n_{s}\in\mathbb{R}$ the
amplitude and $\varphi_{s}=\arg(\tilde{n}_{s})$ the relative phase of the
Fourier component at the index $s$. In the case of phase gratings
($\tilde{n}_{-s}=\tilde{n}_{s}^{*}\in\mathbb{C}$) the Fourier-series of the
real-valued refractive index reads
$n(x)=n_{0}+2\sum_{s=1}^{\infty}|\tilde{n}_{s}|\cos(sGx+\varphi_{s})$. By
determining the Fourier components, i.e. amplitudes and phases of the Fourier
components, the structure is fully retrieved. However, by measuring the
intensities of diffracted signals – as it is usually done in standard
diffraction experiments – only the magnitude of the $s$-th Fourier component
that corresponds to the $s$-th diffraction order, can be obtained, whereas the
phases $\varphi_{s}$ are lost. The latter is called the phase problem of
crystallography or diffraction (see, for instance, Hauptman (1991)). In
principle, the phase problem can be overcome by measuring the phase
differences between the wave incident to a sample of interest and each
diffracted wave interferometrically, which, however, is found to be utterly
difficult in many relevant cases.

Sophisticated techniques have been developed to solve the phase problem and
recover the phase information Hauptman (1991); Woolfson and Fan (1995) either
by so-called non-physical methods (direct methods, isomorphous replacement and
anomalous scattering) or by introducing – sometimes demanding – experimental
procedures to provide a physical solution Chapman  _et al._ (1981); Chang
(1982); Shen (1998); Wolf (2009). Discussions on available techniques and
their development are ongoing (see, for instance, Kleywegt (2000); Cowtan
(2003); Millane (2017); Fally  _et al._ (2019); Donatelli and Spence (2020);
Heacock  _et al._ (2020)). Surprisingly, sticking to traditional diffraction
experiments alone can indeed provide all necessary information in many cases
if experimental settings are chosen carefully, as we will show here. In
particular, we demonstrate by means of simple examples – one-dimensional, non-
sinusoidal phase gratings – how such a physical solution to the phase problem
of diffraction can be obtained with not more than the usual experimental
effort in light and neutron diffraction experiments. The technique is
applicable to samples generating diffraction patterns for any type of
radiation useful for structure determination (X-rays, electrons, laser light,
and neutrons, for instance). We apply our approach to the problem of
determining deviations of the periodic refractive index profiles of
holographic nanoparticle-polymer composite (NPC) gratings Tomita  _et al._
(2016) from their ideal sinusoidal form. The information obtained is useful
for our understanding of holographic formation processes Tomita  _et al._
(2005a); Chikama  _et al._ (2008) and development of versatile holographic
applications such as diffractive optical elements, holographic data storage
and display technology Ohe  _et al._ (1999); Coufal  _et al._ (2000); Crawford
(2003); An  _et al._ (2006).

Holographic gratings were prepared from a photopolymerizable NPC material. As
reported in the past Suzuki and Tomita (2004), SiO2 nanoparticles (with the
average diameter of 13 nm and the bulk refractive index $n_{n}$ of 1.46)
dispersed in a solution of methyl isobutyl ketone are mixed with methacrylate
monomers, 2-methyl-acrylic acid 2-4-[2-(2-methyl-
acryloyloxy)-ethylsulfanylmethyl]-benzylsulfanyl-ethyl ester (the formed
polymer refractive index $n_{p}$ is 1.59 at 589 nm). The doping concentration
of SiO2 nanoparticles was 34 vol.%. Photoinitiator titanocene (Irgacure 784,
Ciba) is also mixed in 1 wt.% with respect to the monomer to provide
photosensitivity in the green. The above chemical mixture is cast on a glass
plate. It is dried in an oven and finally covered with another glass plate,
separated from the first one by spacers of known thickness. We employ a two-
beam interference setup to write an unslanted transmission NPC grating by
superposition of two mutually coherent s-polarized beams of equal light
intensities from a laser diode-pumped frequency-doubled Nd:YVO4 laser
operating at 532 nm. Two NPC gratings were prepared: G1 at grating spacing
$\Lambda$ of $5\,\mu$m (with $\Lambda=2\pi/G$, see Eq. 1) and with the
thickness $d$ of about $13\,\mu$m, and G2 at spacing $\Lambda$ of $1\,\mu$m
and with $d\approx 50\,\mu$m. In the bright regions of the interference
pattern irradiating the sample, the photoinitiator triggers
photopolymerization. Monomer is consumed there by the formation of polymer.
Due to the resulting chemical potential gradient, a mutual diffusion process
of unreacted monomer and nanoparticles sets in, leaving the dark regions
enriched with nanoparticles Tomita  _et al._ (2016). It is important to note
that, while the interference pattern is sinusoidal, the phase grating will be
non-sinusoidal with nonzero higher Fourier components ($|s|>1$ in Eq. 1) due
to the interplay of nonlinear processes governed by the photopolymerization-
driven mutual diffusion process for the formation of spatial density
modulations of the formed polymer and nanoparticles. (Such characteristic
features of the nanoparticle density profile can also be assessed by various
theoretical models Karpov  _et al._ (2000); Tomita  _et al._ (2005a).) Thus,
however elaborate the described production process might be, it is subject to
continuing optimization and fine-tuning of experimental parameters (such as
chemical compositions, nanoparticle concentration, and recording intensity and
time) for the particular application one has in mind for the produced grating.
Therefore, information as accurate as possible about the process outcome is
desired, i.e. information about the redistribution of the nanoparticles
dispersed in monomer under holographic exposure, which would eventually
determine the exact form of the periodic refractive index profile given by Eq.
1. It is the latter we investigate by our proposed method here.

Figure 1:  (a) Optical micrograph of the holographic SiO2 nanoparticle
dispersed NPC grating G1 in transmission mode. (b) Spectral transmittance
distribution along the grating vector of G1 obtained by averaging over lines
of a rectangular region in Fig. 1 (a) of several microns height. (c) Two-
dimensional phase distribution of a grating similar to G1 measured by digital
holographic microscopy. (d) Phase distribution of the grating along the
grating vector obtained by averaging over 100 line scans in Fig. 1 (c). It is
clearly seen that both profiles [(b) and (d)] exhibit higher order Fourier
components.

Figure 1 (a) shows an optical micrograph of G1 taken close to the maximum
contrast position of the Talbot carpet Goodman (2017); Zhou and Burge (2010);
Wen  _et al._ (2013), which can be observed for pure phase gratings using an
optical microscope. It can be seen that SiO2 nanoparticles (narrow, dark
fringes, corresponding to lower refractive index regions) and the formed
polymer (wide, bright fringes, corresponding to higher refractive index
regions) are periodically arranged as a result of holographic assembly of
nanoparticles Tomita  _et al._ (2005b). Their distribution, shown in Fig. 1
(b), is not perfectly sinusoidal. In Fig. 1 (c), a section of a phase map of a
sibling of G1 (also at $\Lambda=5~{}\mu$m) measured by digital holographic
microscopy (DHM-R2100 by Lyncée Tec) is shown. The measurements were made with
a 20X objective at a wavelength of $\lambda=684.9$ nm. Also in this case, it
is clearly seen from the sum of about 100 lines taken along the grooves of the
phase map that the pattern is not purely sinusoidal, as shown in Fig. 1 (d).
Thus, the Fourier series describing both distribution patterns and, therefore,
the refractive-index profiles must certainly include higher order Fourier
terms, at $|s|>1$.

Figure 2:  (a) Plot of typical DE in the Bragg regime as a function of angle
of incidence $\Theta$. (b) Multi-wave interference data of the grating G1 are
obtained diffracting a He-Ne laser of a few mm beam width at a wavelength of
$\lambda=633$ nm, much shorter than the grating spacing of $5\,\mu$m. A
grating is placed on a motorized stage which is rotated about the $y$-axis.

In many diffraction experiments aiming at structure determination, diffraction
by a sample occurs in Bragg regime Moharam  _et al._ (1980): At most two waves
propagate and interfere simultaneously within the sample at all angles of
incidence. One typically measures the dependence of the diffracted intensities
upon the sample’s rotation through angles of incidence $\Theta$ on a screen or
a detector. In Bragg regime, peaks of the diffraction efficiency (DE, $\eta$)
plotted versus $\Theta$ are relatively sharp and do not overlap as can be seen
in Fig. 2 (a). Interference of the two waves propagating within the sample’s
periodic structure (the refractive index modulation given by Eq. 1) results in
energy exchange between diffraction order pairs (the zero order beam and only
one of the $\pm 1,\pm 2,\pm 3,\dots$ orders) as a function of $\Theta$. The
first Born approximation does not hold in such a case, but theories such as
dynamical diffraction theory or Kogelnik’s theory (see, for instance,
Batterman and Cole (1964); Kogelnik (1969)) can be deployed for modelling data
taken in Bragg regime. However, the relative phases $\varphi_{s}$ of Fourier
components cannot be retrieved since multi-wave interference of the
corresponding waves does not occur at any $\Theta$. Now, in order to achieve
such multi-wave interference (also referred to as multi-wave coupling) in
diffraction and, thereby, to determine the relative phases $\varphi_{s}$ of
G1’s structure, we apply visible light of wavelength $\lambda\ll\Lambda$ as a
probe in our diffraction experiment, such that $\Lambda$ is about an order of
magnitude longer than the wavelength of the probe beam. By intentionally
choosing $\lambda$ much shorter than it would seem appropriate, we can leave
the Bragg regime behind and allow for multi-wave interference to occur within
the periodic structure of the sample. From the Bragg equation
$s\lambda=2\Lambda_{s}\sin\Theta_{s}$ ( $\Theta_{s}$ is the Bragg angle for
the Fourier component at the index $s$) it is clear that the diffraction angle
decreases with a decrease in $\lambda$, so that diffraction peaks of a given
width overlap and the experiment can no longer be described in the Bragg
regime. Clearly, more than just two waves propagate and interfere inside the
sample at angles $\Theta$ for which observed peaks overlap considerably [see
Fig. 2 (b)].

The angular dependence data obtained by multi-wave interference diffraction
can only be modelled by multiwave analysis theories such as the rigorous
coupled-wave analysis (RCWA) Moharam  _et al._ (1995). Here, the strategy is
to solve Maxwell’s equations in each of three regions (input, grating, output)
such that the tangential components of the solutions of neighbouring regions
match at these interfaces. Phase information of the Fourier components of the
refractive index profile given by Eq. 1 is inherently included in the RCWA. A
set of coupled-wave equations is solved by calculating eigenvalues and
eigenvectors of a matrix directly related to those Fourier components. The
amplitudes at arbitrary diffraction orders are found by employing the boundary
conditions. In our work RCWA was used to fit the experimental DE defined as
$\eta_{s}=I_{s}/I_{\text{\scriptsize{tot}}}$, with the diffracted intensity
$I_{s}$ at the diffraction order $s$ and the sum of all diffracted intensities
behind the sample $I_{\text{\scriptsize{tot}}}$, and yield the amplitudes and
phases of the corresponding Fourier components at the index $s$ required to
fully reconstruct $n(x)$.

Figure 3:  Measured angular dependence (data points) of the DE for the $\pm
2,\pm 1,0$ orders for G1 at a wavelength of 633 nm. Error bars are much
smaller than the symbols. The $\pm 3$ orders were measured and included in the
fitting procedure, but their DE is small and not shown here. Solid curves are
RCWA fits to the data at seven diffraction orders, yielding $\varphi_{2}$ and
$\varphi_{3}$ (see Table 1).

To measure the angular dependence data, an s-polarized beam of a He-Ne laser
($\lambda=633$ nm) was used to observe diffraction signals from G1 being
rotated about the $y$-axis (perpendicular to the plane of incidence) at angles
$\Theta$ ($-55^{\circ}\ldots+55^{\circ}$), as is shown in Fig. 2 (b). Angular
dependences of the diffracted intensities $I_{s}(\Theta)$ at $s=0,\pm 1,\pm
2,\pm 3$ were recorded by Si-photodiodes placed at the locations of the
diffraction spots observed on a screen. The results are shown in Fig. 3. Due
to the small ratio $\lambda/\Lambda\approx 0.1$ in our experiment, the Bragg
angles $\Theta_{s}=\pm\,3.6^{\circ},\pm\,7.3^{\circ},\pm\,10.9^{\circ}$ at
$s=\pm 1,\pm 2,\pm 3$, respectively, are small, too. All seven observable
diffraction orders overlapped within $\pm 20^{\circ}$ and exhibited
sufficiently large DE to be readily detected. Solid curves in Fig. 3 are
least-squares fits to the RCWA. Fitting was performed including all data
points. The RCWA is in excellent agreement with the experimental data, by
which we can achieve accurate structure determination as is shown below:
Estimations for the amplitudes $|n_{1}|,|n_{2}|$ and $|n_{3}|$ of the Fourier
components as well as their phases $\varphi_{2}$ and $\varphi_{3}$ as obtained
by the RCWA fit are given in the second column of Table 1. The values of
$\varphi_{2}$ and $\varphi_{3}$ were – without loss of generality – obtained
with respect to $\varphi_{1}$ set to zero. The spatial refractive index
profile as calculated from fit parameter estimations is shown in Fig. 4
(solid, black curve). Since the refractive index modulation amplitude of a
recorded NPC grating is proportional to $n_{n}-n_{p}$ Suzuki  _et al._ (2002),
the SiO2 rich regions correspond to the low refractive index regions.
Comparison with the data of Fig. 1 (b) (filled, red circles in Fig. 4) makes
the qualitative agreement obvious, thereby demonstrating the validity of our
approach. Even the slight asymmetry in the profile (slight flattening on the
left side of each peak), which may be caused by spatially nonuniform lateral
shrinkage during holographic exposure, is captured by the RCWA analysis of the
diffraction data as can be seen from the inset in Fig. 4.

grating | G1 (light, $\lambda=633\,$nm) | G2 (neutrons, $\lambda=1.7\,$nm)
---|---|---
$|n_{1}|$ | $4.9382(37)\times 10^{-3}$ | $2.592(28)\times 10^{-6}$
$|n_{2}|$ | $1.072(15)\times 10^{-3}$ | $5.03(58)\times 10^{-7}$
$|n_{3}|$ | $1.72(48)\times 10^{-4}$ | n. a.
$\varphi_{1}$ | $:=0$ | $:=0$
$\varphi_{2}$ | $1.0581(39)\pi$ | $0.995(27)\pi$
$\varphi_{3}$ | $0.37(15)\pi$ | n. a.
Table 1:  Amplitudes and relative phases for the Fourier coefficients of
refractive index profiles for G1 and G2. The values were determined by RCWA
fits to the data with $\varphi_{1}$ set to zero (see text). Resulting profiles
obtained from light and neutron diffraction measurements are shown in Figs. 4
and 6, respectively. Figure 4:  Plotted are the data points as shown in Fig. 1
(b) (filled, red circles). The refractive index profile of G1 from the RCWA
fit-parameter estimation (as given in Table 1) was rescaled and shifted
(solid, black curve) for comparison to the data. Peaks correspond to the
polymer ridges of the structure. Inset: Third peak from the left, magnified.
Even detailed features of the micrograph data are reproduced, as shown in the
inset.

Next, the refractive index profile of G2 ($\Lambda=1\,\mu\textrm{m}$,
$d\approx 50\,\mu$m) is to be retrieved. Accurately resolving such structures
by an optical microscope is possible but not straightforward, which underlines
the potential impact of our solution to the phase problem, especially for
refractive-index profiles, which have structural dimensions of the order of
light wavelengths. In the case of G2, Kogelnik’s theory Kogelnik (1969) allows
to roughly estimate the expected peak width (the angular distance between the
minima adjacent to the Bragg angle $\Theta_{1}$ at $s=1$ in an angular
dependence plot of DE) as $2n_{0}\Lambda/d\approx 3.4^{\circ}$ while one
obtains $\Theta_{1}\approx 17^{\circ}$ at $\lambda=633$ nm and
$\Theta_{1}\approx 10^{\circ}$ at a UV wavelength of $\lambda=351\,$nm, say.
Thus, it is difficult to produce peak overlap to transfer the diffraction
process from the Bragg regime to the multi-wave interference regime for the
determination of the refractive index profile of G2, by use of table-top laser
light sources. However, a key point of our approach is that depending on the
particular material class investigated, there might be other kinds of
radiation available for diffraction experiments to obtain data similar to the
ones shown in Fig. 3. For instance, since the refractive index profile is
produced by the density modulation of one of the involved material components
(nanoparticles, in our case), it is known that small angle neutron scattering
(SANS; see, for instance, Willis and Carlile (2009)) provide well-established
tools Fally  _et al._ (2010): For a typical de Broglie wavelength of neutrons
in a SANS experiment of 1 nm, say, one may expect Bragg angles
$\Theta_{1}\approx\lambda/(2\Lambda)$ of the order of $0.03^{\circ}$, which
can be detected with state-of-the-art SANS instruments, thanks to long sample-
detector distances (up to 20 m) and sufficient spatial resolution of
detectors. Thus, considerable peak-overlap for G2 and, therefore, multi-wave
interference is achievable with neutrons.

The neutron experiment was performed by use of the instrument SANS-I of the
SINQ neutron source of Paul-Scherrer Institute in Villigen, Switzerland. Cold
neutrons at a mean wavelength of $1.7$ nm (width of wavelength distribution
$\Delta\lambda/\lambda\approx 10\%$) were diffracted from G2. The beam
divergence was limited to about 1 mrad, using collimation slits. The
diffracted intensities were measured using a two-dimensional detector of
$7.5\times 7.5\,$mm2 pixel size. To adjust the peak width (estimated by
$2\Lambda/d$, see our above discussion) and the peak height for our purpose,
G2 was tilted around its grating vector by $\zeta\approx 60^{\circ}$ which
increases the effective thickness to $d\rightarrow d/\cos\zeta$, i. e. by a
factor of two Somenkov  _et al._ (1978); Fally  _et al._ (2010). The results
are shown in Fig. 5. Similarly to the DE of G1 shown in Fig. 3, the observed
diffraction took place in the multiple-wave interference regime, that is, many
(five, in this case) diffraction orders are observable within the angular
range of $\Theta=\pm 0.25^{\circ}$. Applying the RCWA fit, we were, again,
able to extract the full information (amplitudes _and_ phases) of the
grating’s Fourier components up to the 2nd-order. They are given in Table 1
(third column). Simulations of the DE assuming $\varphi_{2}=0$ instead of the
experimentally obtained value are also plotted (red, dashed curves) in Fig. 5.
The disagreement of simulation at $\varphi_{2}=0$ and data is easily resolved
in a standard SANS experiment. The neutron-refractive index profile of G2 –
corresponding to the spatial density distribution of formed polymer and
nanoparticles – as calculated from the RCWA parameter estimation is shown in
Fig. 6 together with its counterpart, assuming $\varphi_{2}=0$, for
comparison. Note that the refractive index modulation amplitude for neutrons
is proportional to $-(b_{n}\rho_{n}-b_{p}\rho_{p})$, where $b_{n}$ ($b_{p}$)
and $\rho_{n}$ ($\rho_{p}$) are the mean bound coherent scattering length
Willis and Carlile (2009) and the atomic density of nanoparticles (the formed
polymer), respectively. Since $b_{n}\rho_{n}$ is larger than $b_{p}\rho_{p}$
for SiO2 nanoparticles Fally  _et al._ (2010), the SiO2 rich regions
correspond to lower refractive index regions in Fig. 6, as similar to Fig. 4.
No third order component was measurable for G2 and therefore the profile in
Fig. 6 is symmetric, in contrast to the profile of G1. We attribute this
difference to the mutual diffusion process, which is much more significant in
G2 due to the smaller grating spacing.

Of course, it is – in many cases – more convenient to apply optical microscopy
for the estimation of the refractive-index profile, but resolution limits of
this method are laborious to overcome in the range of typical structure
constants of several hundred nanometers. Furthermore, electron microscopy and
physico-chemical analyses often depend on sample preparation techniques –
producing thin slices or breaking samples to look at surfaces – that are
cumbersome or even unreliable in the sense that they could, in the worst case,
mechanically alter the structure to be investigated. Our technique, however,
can provide sufficient resolution in bulk so long as a suitable state-of-the-
art instrument, meeting the wavelength requirements to work in multi-wave
interference regime for a structure of interest, can be accessed at one of the
many facilities worldwide. Resolution is limited by the –usually excellent –
detector sensitivity and background suppression necessary for reliable
observation of weak higher-order diffraction signals. The RCWA analysis
(programmed in the Python programming language) carried out here on an off-
the-shelf PC can be applied to more complicated structures by using more
powerful computation infrastructure, nowadays available at many institutions.

Figure 5:  Measured angular dependence (data points) of DE for G2 at the $\pm
2,\pm 1,0$ orders and at the neutron wavelength of $1.7$ nm. Five-wave
interference was observed near the normal incidence. An RCWA fit (solid
curves) again allows to retrieve amplitudes and phases of the grating’s
Fourier components. The resulting fit parameter estimation is given in Table
1. To illustrate the importance of the correct phase value, simulations for
$\varphi_{2}=0$ (dashed, red curves) instead of the experimentally obtained
value $\varphi_{2}=0.995\,\pi$ are added. Figure 6:  The neutron-refractive
index profile of G2 calculated from the RCWA fit-parameter estimation as given
in Table 1 (solid, black curve). The neutron-refractive index profile of G2 is
also plotted for $\varphi_{2}=0$ (dashed, red curve) instead of the
experimentally obtained value $\varphi_{2}=0.995\,\pi$.

Finally, we note that the technique presented here is closely related to
previous proposals to employ multi-wave interference (see, for instance,
Chapman  _et al._ (1981); Chang (1982); Woolfson and Fan (1995); Shen (1998)).
The fundamental difference between those and our approach is, however, that
apart from a careful choice of wavelengths, only standard diffraction
procedures are applied here. No extra data needs to be taken. Our approach can
be seen as generalization of what is discussed in the recent Ref. Heacock  _et
al._ (2020).

In summary, we have demonstrated the determination of refractive-index
profiles by diffraction from one-dimensional holographic phase gratings
recorded in SiO2 nanoparticle-dispersed NPC films. We have shown that full
phase retrieval can be made without the need for extra data collection
schemes. In particular, by choosing the probe wavelength some orders of
magnitude shorter than the structural dimensions being investigated, Bragg
regime diffraction can be turned into diffraction in the regime of multi-wave
interference for a given sample. Fitting of a multi-wave coupling model such
as the well-known RCWA to the resulting angular dependence data allows for the
accurate determination not only of amplitudes but also of phases of the
Fourier components of generic profiles of refractive index modulation. When
dealing with photonic structures (like gratings) it is a common practice to
avoid multi-wave interference for ease of the analysis, in particular when it
comes to the evaluation of the underlying parameters. Here, we have shown that
getting rid of this habit and embracing the usually shunned multi-wave
interference regime can pay off at the prize of a relatively low increase in
experimental and analytical complication. Our technique is a purely physical
solution to the phase problem of diffraction.

We would like to thank G. Heuberger for preliminary experimental work. J. G.
and Y. T. would like to acknowledge the financial supports by Ministry of
Education, Culture, Sports, Science, and Technology of Japan under grant No.
25-03052 and No. 15H03576.

## References

  * Hauptman (1991) H. A. Hauptman, Rep. Prog. Phys. 54, 1427 (1991).
  * Woolfson and Fan (1995) M. Woolfson and H.-f. Fan,  _Physical and non-physical methods of solving crystal structures_ (Cambridge University Press, Cambridge, UK, 1995).
  * Chapman  _et al._ (1981) L. D. Chapman, D. R. Yoder, and R. Colella, Phys. Rev. Lett. 46, 1578 (1981).
  * Chang (1982) S.-L. Chang, Phys. Rev. Lett. 48, 163 (1982).
  * Shen (1998) Q. Shen, Phys. Rev. Lett. 80, 3268 (1998).
  * Wolf (2009) E. Wolf, Phys. Rev. Lett. 103, 075501 (2009).
  * Kleywegt (2000) G. J. Kleywegt, Acta Crystallogr. D 56, 249 (2000).
  * Cowtan (2003) K. Cowtan, “Phase problem in x-ray crystallography, and its solution,” in  _eLS_ (American Cancer Society, 2003).
  * Millane (2017) R. P. Millane, Acta Crystallogr. A 73, 140 (2017).
  * Fally  _et al._ (2019) M. Fally, J. Klepp, M. A. Ellabban, H. Eckerlebe, P. K. Pranzas, J. Guo, and Y. Tomita, Proc. SPIE 11030, 110300I (2019).
  * Donatelli and Spence (2020) J. J. Donatelli and J. C. H. Spence, Phys. Rev. Lett. 125, 065502 (2020).
  * Heacock  _et al._ (2020) B. Heacock, D. Sarenac, D. G. Cory, M. G. Huber, J. P. W. MacLean, H. Miao, H. Wen, and D. A. Pushin, IUCrJ 7, 893 (2020).
  * Tomita  _et al._ (2016) Y. Tomita, E. Hata, K. Momose, S. Takayama, X. Liu, K. Chikama, J. Klepp, C. Pruner, and M. Fally, J. Mod. Opt. 63, S1 (2016).
  * Tomita  _et al._ (2005a) Y. Tomita, N. Suzuki, K. Furushima, and Y. Endoh, Proc. SPIE 5939, 593909 (2005a).
  * Chikama  _et al._ (2008) K. Chikama, K. Mastubara, S. Oyama, and Y. Tomita, J. Appl. Phys. 103, 113108 (2008).
  * Ohe  _et al._ (1999) Y. Ohe, M. Kume, Y. Demachi, T. Taguchi, and K. Ichimura, Polym. Adv. Technol. 10, 544 (1999).
  * Coufal  _et al._ (2000) H. J. Coufal, D. Psaltis, and G. T. Sincerbox, eds., in  _Holographic Data Storage_ (Springer-Verlag Berlin Heidelberg New York, 2000).
  * Crawford (2003) G. P. Crawford, Opt. Photon. News 14, 54 (2003).
  * An  _et al._ (2006) J. W. An, D. D. Do, N. Kim, and K. Y. Lee, IEEE Photon. Technol. Lett. 18, 788 (2006).
  * Suzuki and Tomita (2004) N. Suzuki and Y. Tomita, Appl. Opt. 43, 2125 (2004).
  * Karpov  _et al._ (2000) G. Karpov, V. Obukhovsky, T. Smirnova, and V. Lemeshko, Opt. Commun. 174, 391 (2000).
  * Goodman (2017) J. W. Goodman,  _Introduction to Fourier Optics_ (W. H. Freeman and Company, New York, USA, 2017).
  * Zhou and Burge (2010) P. Zhou and J. H. Burge, Appl. Opt. 49, 5351 (2010).
  * Wen  _et al._ (2013) J. Wen, Y. Zhang, and M. Xiao, Adv. Opt. Photon. 5, 83 (2013).
  * Tomita  _et al._ (2005b) Y. Tomita, N. Suzuki, and K. Chikama, Opt. Lett. 30, 839 (2005b).
  * Moharam  _et al._ (1980) M. Moharam, T. Gaylord, and R. Magnusson, Optics Communications 32, 14 (1980).
  * Batterman and Cole (1964) B. W. Batterman and H. Cole, Rev. Mod. Phys. 36, 681 (1964).
  * Kogelnik (1969) H. Kogelnik, Bell Syst. Tech. J. 48, 2909 (1969).
  * Moharam  _et al._ (1995) M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, J. Opt. Soc. Am. A 12, 1077 (1995).
  * Suzuki  _et al._ (2002) N. Suzuki, Y. Tomita, and T. Kojima, Appl. Phys. Lett. 81, 4121 (2002).
  * Willis and Carlile (2009) B. T. M. Willis and C. J. Carlile,  _Experimental Neutron Scattering_ (Oxford University Press, Oxford, UK, 2009).
  * Fally  _et al._ (2010) M. Fally, J. Klepp, Y. Tomita, T. Nakamura, C. Pruner, M. A. Ellabban, R. A. Rupp, M. Bichler, I. Drevenšek-Olenik, J. Kohlbrecher, H. Eckerlebe, H. Lemmel, and H. Rauch, Phys. Rev. Lett. 105, 123904 (2010).
  * Somenkov  _et al._ (1978) V. A. Somenkov, S. S. Shilstein, N. E. Belova, and K. Utemisov, Solid State Commun. 25, 593 (1978).