11institutetext:  Institut de Physique Théorique, CEA - IPhT, CNRS, URA 2306
F-91191 Gif-sur-Yvette, France
Faculty of Engineering, Takushoku University - Hachioji, Tokyo 193-0985, Japan
Research School of Physical Sciences, Australian National University -
Canberra, ACT 0200, Australia

Conformation (statistics and dynamics) Structural transitions in nanoscale
materials Analytical theories

# Ribbon polymers in poor solvents: layering transitions in annular and
tubular condensates

Y. Y. Suzuki 1122 D. R. M. Williams 33112233

###### Abstract

We study the structures of a ribbon or ladder polymer immersed in poor
solvents. The anisotropic bending rigidity coupled with the surface tension
leads ribbon polymers to spontaneous formation of highly anisotropic
condensates in poor solvents. Unlike ordinary flexible polymers these
condensates undergo a number of distinct layering transitions as a function of
chain length or solvent quality, and the size of condensates becomes non-
monotonic function of chain length. We show that the fluctuations of the
condensates are in general small and these condensates are stable.

###### pacs:

36.20.Ey

###### pacs:

64.70.Nd

###### pacs:

87.15.ad

Recently considerable attention has been devoted to the physics of
semiflexible polymers, i.e., chains which have significant bending stiffness.
These kinds of polymers exhibit interesting physics in the form of liquid
crystalline phases and are of great importance in biology and biophysics. For
example, DNA is the most well-known semiflexible polymer. Much of the interest
has focused on semiflexible chains with isotropic bending elasticity[1, 2, 3,
4, 5, 6, 7]. Macroscopically this corresponds to modeling the chain as a
cylinder with a circular cross section. In many cases, however, the cross
section is highly anisotropic[8, 9, 10, 11, 12].

The anisotropic polymer can be modeled as a ribbon, of length $L$, width $w$,
and thickness $t$ assuming $L\gg w>t$ (fig.1) with anisotropic bending
elasticity. Indeed, if the ribbon is composed of an isotropic material, the
ratio of the elastic constants for bend in the easy and hard directions is
$\epsilon_{e}/\epsilon_{h}=(t/w)^{2}$ [13] . This strong dependence on the
ratio of $t/w$ implies that the ribbon polymer shows a very anisotropic
elastic response. On the chemical scale, the chains do not normally consist of
an isotropic material, and details of chemical bonding are important for bend,
thus, even higher ratios of $\epsilon_{e}/\epsilon_{h}$ are possible.

[width=7cm]fig1.eps

Figure 1: Model for anisotropic semiflexible polymer

Some of the solution properties, in particular the liquid-crystalline
behavior, have been studied theoretically [10]. Here we study the structure of
a single ribbon polymer immersed in a poor solvent. Assuming the system is
well below the $\Theta$ temperature, the chain collapses completely and forms
a condensate which is effectively a polymer melt.

For isotropic semiflexible chains such as DNA, the condensate often forms a
toroid [1, 2, 14, 15, 16, 17, 18, 19, 20, 21] and sometimes a globule[22, 23,
24, 25] through intermediate states[28, 29, 30]. For ribbon polymers, we show
in this paper that the condensates forms an annulus. The annulus formation
manages to avoid bend in the hard direction, while reducing contact of the
chain with the solvent. The novel character of these annular condensates is
that multilayer annuli can form and they undergo sudden changes as a function
of chain length or surface tension. This is in marked contrast to the case of
ordinary flexible or semiflexible polymers where no dramatic changes occur.
This quantization is a direct consequence of the anisotropic elasticity.

To model the system, we introduce two surface tensions $\gamma_{h}$ and
$\gamma_{e}$ which corresponds to the surface energies of contact between the
top (narrow surface) of the polymer and the solvent, and the surface energies
of contact between the side (flat and wide surface) of the polymer and the
solvent, respectively. If the chain consists of anisotropic molecules, any
ratio of $\gamma_{e}/\gamma_{h}$ is possible in principle. In this paper,
however, we assume $\gamma_{e}/\gamma_{h}\approx 1$.

A large number of parameters appear in this model. Once final results have
been obtained, it is useful to ignore numerical prefactors and to make some
crude estimates. Then, we write $w\approx t\approx a$ and
$\gamma_{e}\approx\gamma_{h}\approx gkT/a^{2}$, where $g$ is of order unity
and $a$ is a length of order of one Angstrom. It is convenient to define two
bare persistence lengths, $l_{e}\equiv\epsilon_{e}/kT$ and
$l_{h}\equiv\epsilon_{h}/kT$ which are essentially the scales on which a free
chain bends in the two directions due to thermal fluctuations.

We begin by examining the shape of condensates when we slowly lengthen the
chain. For short chains, bend costs elastic energy, then the chain more-or-
less remains as a rod. For slightly longer chains, however, the overlap
induced by bend compensates the elastic energy by reducing its surface energy,
therefore, the chain may form a ring. For a circular ring of radius $R$ it is
easy to show that the bending energy is ${1\over 2}\epsilon_{e}L/R^{2}$. The
change in surface energy in bending into a ring is $4\pi
Rw\gamma_{e}-2Lw\gamma_{e}$. Minimising it over $R$ and equating the free
energies gives the critical length,
$L^{*}=(3\sqrt{3}\pi/2)\sqrt{\epsilon_{e}/(\gamma_{e}w)}$ at which ring
formation occurs. The chain overlaps by $1/2$ turn, at $L=L^{*}$. Our interest
here is in the case $L\gg L^{*}$ so the chain winds around many times.

The simplest way the chain can pack is to form a disk with a hole in the
middle (annulus) as shown in fig.2. This annulus has height $w$ and inner and
outer radii $R_{i}$ and $R_{o}$, respectively. The total energy consists of
three terms: (1) surface energy of the top and bottom surfaces, (2) surface
energy of the inner and outer exposed sections of the annulus, (3) bending
energy of the chain. This last term is due to bending only in the easy
direction.

[width=6cm]fig2.eps

Figure 2: Annular condensate for short ribbon polymer in poor solvent

Because the volume of the chain is conserved, $R_{o}$ and $R_{i}$ are related
by

$tLw=w\pi(R_{o}^{2}-R_{i}^{2}).$ (1)

There is thus only one variable to minimize over (say $R_{i}$). Also, since
the area of the upper and lower surfaces is fixed (i.e., both are always of
area $tL$), there are really only two relevant terms in the total energy. The
lateral surface energy is easily calculated as

$2\pi\gamma_{e}w(R_{o}+R_{i}).$ (2)

The bending energy is calculated by

$\frac{\epsilon_{e}}{2}\int_{0}^{L}ds\left(\frac{d\theta}{ds}\right)^{2},$ (3)

where $s$ is the arc length and $\theta$ is the angle between the tangent of
the ribbon curve and some fixed direction in the plane of the disk.

In one complete turn of the ribbon, $\theta$ increases by $2\pi$ and $s$
increases by $2\pi r$ where $r$ is the distance of the chain from the center
of the annulus. This gives $d\theta/ds=1/r$. During the same change in $s$ the
radius increases by one ribbon thickness, i.e., $dr/ds=t/(2\pi r)$. This
allows us to write the bending energy as
$\pi\epsilon_{e}t^{-1}\ln(R_{o}/R_{i})$. The total energy is then

$U=\frac{\pi\epsilon_{e}}{t}\ln\left({R_{o}\over
R_{i}}\right)+2\pi\gamma_{e}w(R_{o}+R_{i}),$ (4)

subject to $R_{o}^{2}=R_{i}^{2}+tL/\pi$.

It is convenient at this point to introduce two characteristic lengths. One is
$\Lambda\equiv\epsilon_{e}/(2\gamma_{e}wt)$, and the other is
$\Gamma\equiv\sqrt{Lt/\pi}$. For long chains, $\Gamma$ is much larger than
$\Lambda$, and we shall see $R_{i}\approx\Lambda$ and $R_{o}\approx\Gamma$. In
this limit, $\Lambda\ll\Gamma$, the annulus is “fat”, i.e., it has a small
hole in the middle. In the opposite limit, $\Lambda\gg\Gamma$, the hole is
large and the annulus is “thin”. It looks like a ring.

In the long-chain limit, $R_{i}\ll\Gamma$, the total energy is approximated by

$U={\rm const.}+2\pi\gamma_{e}wR_{i}-\frac{\pi\epsilon_{e}}{t}\ln\left({\pi
R_{i}\over\sqrt{Lt}}\right),$ (5)

which has a minimum at $R_{i}^{*}=\Lambda$. Using our crude estimates give
$R_{i}^{*}\approx l_{e}/g$, i.e., the internal radius is somewhat less than
the bare persistence length. In this limit the internal radius of the annulus
is independent of chain length, whereas the external radius is approximately
$R_{o}\approx(tL/\pi)^{1/2}$. A more accurate expression for $R_{i}^{*}$ can
be found by expanding in ascending powers of $\alpha\equiv\Lambda/\Gamma$.
This gives $R_{i}^{*}=\Lambda(1-\alpha+{3\over 2}\alpha^{3}-...)$. From this
it is clear that, even though $R_{i}^{*}$ becomes asymptotically $\Lambda$,
the approach is rather slow as $L^{-1/2}$.

In the limit $\Lambda\gg\Gamma$, we can again expand the total energy and find
the minimum at $R_{i}^{*}=(\Gamma^{2}\Lambda/2)^{1/3}$, or more crudely,
$R_{i}^{*}\approx(Ll_{e}a/g)^{1/3}$. The crossover between the thin and fat
regimes occurs at $\Gamma=\Lambda$. A crude estimate of this gives
$L_{c}=a(L^{*}/a)^{4}$. This implies that the chain must be rather long to be
in the fat regime.

The results just derived correspond to the case of one layer. The novel point
about this system is that it can choose to form more than one layer. This
decreases the surface area exposed by the top and bottom surfaces but comes at
the expense of extra surface energy on the inside and outside surfaces and
extra bend both in the hard and easy directions. For second layer to form the
chain must bend in the hard direction. The best way of doing this is to join
the two layers together at the outer radius (fig.3).

[width=4cm]fig3.eps

Figure 3: tubular condensate for long ribbon polymer in poor solvent

We need to calculate the bending energy of the chain in the hard direction.
The chain must bend between the two layers. It must jump a distance $w$ in the
vertical direction while winding a distance $\Omega=2\pi R_{o}$ in arc length.
If we let $\phi$ be the angle made by the chain with the plane of the lower
annulus then the bending energy can be written as

$U_{h}=\frac{\epsilon_{h}}{2}\int_{0}^{\Omega}ds\left({d\phi\over
ds}\right)^{2}.$ (6)

This is to be minimized subject to a constraint that the chain jumps up a
distance $w$, i.e., $\int_{0}^{\Omega}ds\phi=w$. Here we have made the very
good approximation that the angle $\phi$ is always small. Using the method of
Lagrange multipliers and the Euler-Lagrange equations we find that the minimum
energy trajectory coresponds to the solution of the differential equation:
$d^{2}\phi/ds^{2}={\rm const}$. Using the boundary conditions
$\phi(0)=0=\phi(\Omega)$ the constraint gives

$\phi(s)=\frac{6w}{\Omega}\left\\{\frac{s}{\Omega}-\left(\frac{s}{\Omega}\right)^{2}\right\\},$
(7)

and the bending energy becomes

$U_{h}={6\epsilon_{h}w^{2}\over\Omega^{3}}.$ (8)

It is fairly clear from this formula that the chain wants to spread the bend
over as much length as possible, so that it is best to bend at the outer
radius of the annulus and $\Omega=2\pi R_{o}$.

We now examine the transition from one to two layers. The one-layer energy is

$U_{1}={\pi\epsilon_{e}\over t}\ln\left({R_{o}\over
R_{i}}\right)+2\pi\gamma_{e}w(R_{o}+R_{i})+2\gamma_{h}Lt,$ (9)

subject to $R_{o}^{2}=R_{i}^{2}+tL/\pi$. The two-layer energy is

$\displaystyle U_{2}$ $\displaystyle=$ $\displaystyle{2\pi\epsilon_{e}\over
t}\ln\left({R_{o}\over R_{i}}\right)+4\pi\gamma_{e}w(R_{o}+R_{i})$ (10)
$\displaystyle+\gamma_{h}Lt+{6\epsilon_{h}w^{2}\over(2\pi R_{o})^{3}},$

subject to $R_{o}^{2}=R_{i}^{2}+tL/2\pi$. Here $\gamma_{h}Lt$ is the surface
tension of the “top” of annulus.

The transition from one to two layers occurs when $U_{1}=U_{2}$. To solve this
equation we need to minimize $U_{1}$ over $R_{i}$ and similarly for $U_{2}$.
In principle different $R_{i}$ will be associated with the two geometries.
However, in the limit where $R_{i}\ll(Lt)^{1/2}$ the result is the same for
both, i.e., $R_{i}^{2}\approx\epsilon_{e}/(2\gamma_{e}wt)$. This gives

$\displaystyle U[n]$ $\displaystyle=$
$\displaystyle\frac{n\pi\epsilon_{e}}{t}\left\\{1+\ln\left({2\gamma_{e}tw\over\epsilon_{e}}\sqrt{Lt\over(n\pi)}\right)\right\\}$
(11) $\displaystyle+\frac{2}{n}\gamma_{h}Lt+2\gamma_{e}\pi
w\sqrt{nLt\over\pi}$
$\displaystyle+{3(n-1)\epsilon_{h}w^{2}\over\sqrt{2}\pi^{3}}\left(Lt\over\pi\right)^{-3/2},$

where $n$ is the number of layers ($1$ or $2$). Equating the two energies for
$n=1$ and $n=2$, and ignoring the logarithmic factors, gives

$\displaystyle{\pi\epsilon_{e}\over t}-\gamma_{h}Lt+2(\sqrt{2}-1)\pi
w\gamma_{e}\sqrt{Lt\over\pi}$
$\displaystyle+{3\epsilon_{h}w^{2}\over\sqrt{2}\pi^{3}}\left(Lt\over\pi\right)^{-3/2}=$
$\displaystyle 0.$ (12)

The four terms in this equation represent the four kinds of energy differences
between the two layer and one layer system. In order these are: Easy bend;
top/bottom surfaces; side surfaces; hard bend. The driving force for the
transition to two layers is the top/bottom surface energy. This is opposed by
the other three terms. We are mainly interested here in the case where the
hard bend is the dominant penalty. In this case the transition occurs when

$L^{5}\approx{\epsilon_{h}^{2}w^{4}\over(\gamma_{h}^{2}t^{5})},$ (13)

where we have ignored numerical prefactors. Even more approximately this can
be written $L\approx a(l_{h}/a)^{2/5}$ where $l_{h}$ is the persistence length
in the hard direction. Provided $l_{h}\gg l_{e}$ (it is true when $L\gg
R_{i}$) our approximations are valid.

At the transition the jump in $R_{i}$ is small, since $R_{i}$ is about the
same for both layers. However, there is a sudden jump in $R_{o}$ by a factor
of $1/\sqrt{2}$. In the case of a three layer system, the chain must jump once
on the outside radius and once on the inside radius. Since the jump on the
inside radius must take place over a shorter amount of arc length, its energy
cost is expensive. The energy for $n$ layers is then

$\displaystyle U[n]$ $\displaystyle=$ $\displaystyle{n\pi\epsilon_{e}\over
t}\ln\left(R_{o}\over R_{i}\right)$ (14)
$\displaystyle+n(2\pi\gamma_{e}w)(R_{o}+R_{i})+{2\over n}\gamma_{h}Lt$
$\displaystyle+\frac{6\epsilon_{h}w^{2}}{(2\pi)^{3}}\left(\left[{n\over
2}\right]R_{o}^{-3}+\left[{{n-1}\over 2}\right]R_{i}^{-3}\right),$

with $R_{o}=(R_{i}^{2}+n^{-1}tL/\pi)^{1/2}$, and where $[x]$ indicates the
integer part of $x$.

In order to show the layering transition of equilibrium form in a fat
condition, in fig.4 we plot the outer radius of the condesates $R_{o}$, the
radius of the inner hole of the condensates $R_{i}$, and the quantized height
of the condensate $nw$. This fat condition is always achieved in the long
chain limit. For large $L$, the size of the condensate is roughly proportional
to $\sqrt{L}$.

[width=6cm]fatdisk.eps

Figure 4: Layering transition in a fat condition $\Gamma/\Lambda\approx 10$.
The horizontal axis is the polymer length $L/w$. The upper and the lower
curves represent $R_{o}$ and $R_{i}$ respectively and the middle step line
represents the height of the tubular condensate in unit of ribbon width $w$.

Below the medium chain length, there exists a certain parameter region where
the condensate is thin. We show in fig.5 behavior of the layering transition
in a thin condition.

[width=6cm]thindisk.eps

Figure 5: Layering transition in a thin condition $\Gamma/\Lambda\approx 0.5$.
The horizontal axis is the polymer length $L/w$. The upper two curves
represent $R_{o}$ and $R_{i}$ respectively and the lower step line represents
the height of the tubular condensate in unit of ribbon width $w$.

We have assumed here that the system can form a well-defined annulus. It is
important and interesting to examine the fluctuations around the equilibrium
shape of the system. We will concentrate here on the case of a single layer.
There are two types of fluctuations. The first kind is fluctuations of the
shape in the radial direction, i.e., distortion from a circle. We discuss it
later.

The second kind is fluctuations out of the plane which cause two kinds of
energy penalty. The first is arises from hard bend and the energy penalty is
$\epsilon_{h}/2\int_{0}^{L}ds(d\phi/ds)^{2}$, where $\phi$ is the angle made
by the tangent to the ribbon out of the plane of the annulus. The second
arises from the surface area exposed due to differences in the height $z(s)$
of the layers. This gives rise to an energy penalty $\gamma_{e}\int_{2\pi
R_{1}}^{L}ds|z(s)-z\\{s-2\pi r(s)\\}|.$ Noting that
$z(s)=\int_{0}^{s}ds^{\prime}\phi(s^{\prime})$ and expanding for small
gradients gives $2\pi\ \gamma_{e}\int_{2\pi R_{1}}^{L}ds|\phi(s)|r(s).$ The
total energy penalty is then

$\frac{\epsilon_{h}}{2}\int_{0}^{L}ds\left(d\phi\over
ds\right)^{2}+2\pi\gamma_{e}\int_{2\pi R_{1}}^{L}ds\ r(s)|\phi(s)|.$ (15)

This kind of total energy is peculiar because of the form of the surface
tension term, $|\phi(s)|$. This term is first order and not differentiable
about the equilibrium position $\phi=0$. This means that the usual eigenmode
techniques we might have used to analyse the system are not applicable.

However, a good estimate of the size of the fluctuations and in particular
$\delta z\equiv\sqrt{\langle(z(L)-z(0))^{2}\rangle}\ $ can be found by mapping
the problem onto that of a semiflexible chain in a nematic field. From now on
we ignore all numerical prefactors. For our purposes it is convenient to use
an approach developed by Odijk [31, 32]. He considers a wormlike chain
directed by a parabolic potential $\mu\phi^{2}$, where $\mu$ is a constant.
The chain has a typical angular deviation
$\delta\phi\equiv\sqrt{\langle\phi^{2}\rangle}$ . Such a chain can be viewed
as a series of steps of size $\lambda_{\parallel}=l_{h}(\delta\phi)^{2}$ in
the plane of the annulus and
$\lambda_{z}=\lambda_{\parallel}\delta\phi=l_{h}(\delta\phi)^{3}$
perpendicular to the annulus. A chain of length $L$ undergoes an anisotropic
random walk of $L/\lambda_{\parallel}$ steps in the plane, with each step
travelling a distance $\lambda_{z}$ perpendicular to the plane. This implies
the deviation out of the plane is $\delta z=\sqrt{Ll_{h}}(\delta\phi)^{2}$.

Odijk shows that the strength of the potential and the angular deviations are
related by $(\delta\phi)^{4}=kT/(\mu l_{h})$, so that $\delta
z=\sqrt{kTL/\mu}$. This result cannot be applied directly to our problem,
since we have a confining potential of the form $\gamma_{e}r(s)|\phi(s)|$.
However, the details of the shape of the potential cannot matter for the gross
behaviour of the system. In particular we can replace our potential by a
harmonic one, provided we choose $\mu$ self-consistently, i.e., we choose
$\mu$ so that for a typical value of the angular deviation the two potentials
coincide. Equating the two potentials gives
$\gamma_{e}r(s)|\phi|=\mu\phi^{2}$. This together, with the relation
$(\delta\phi)^{4}=kT/(\mu l_{h})$ gives
$\mu=(\gamma_{e}r(s))^{4/3}(l_{h}/kT)^{1/3}$. Note that since $r$ varies along
the chain, the effective confining potential $\mu$ also varies along the
chain. We can now integrate over the whole annulus to get the variation in
height $\delta z$. Between $r$ and $r+dr$ there is a total length $rdr/t$ and
hence the square of the step in the $z$ direction is
$dz^{2}=(kT/(\gamma_{e}r))^{4/3}l_{h}^{-1/3}t^{-1}rdr$. Integrating over the
whole annulus gives:

$\delta
z=\frac{kT}{(\gamma_{e}l_{h})^{2/3}}\sqrt{\frac{l_{h}}{t}}\sqrt{R_{o}^{2/3}-R_{i}^{2/3}}.$
(16)

Note that this deviation grows only very slowly with the chain length. In,
particular, for very long chains we have $\delta
z=(kT/(\gamma_{e}l_{h}))^{2/3}(l_{h}/t)^{1/2}t^{1/6}L^{1/6}$. Setting $t=a$
and $\gamma_{e}=gkT/a^{2}$ gives $\delta z\approx ag^{-2/3}(L/l_{h})^{1/6}$.
Thus, in general the overall fluctuations will be less than, or of the same
order as the annular thickness, $w$, and the annulus will be well-defined and
relatively flat. This is the conclusion for one layer.

For many layers this conclusion remains true, since in that case the effective
persistence length scales with the number of layers and the radius also scales
with the number of layers. This implies that the absolute deviation in
thickness remains roughly constant and the relative variation decreases as
more layers are added.

We now discuss fluctuations of radial direction. The simplest mode is one
where the inner radius changes from its equilibrium value, $R_{i}^{*}$. In the
long-chain limit, where $R_{i}^{*}\ll\sqrt{Lt}$ the energy penalty for
deviating from equilibrium is

$\delta U={\pi\epsilon_{e}\over
2t}\left(\frac{R_{i}-R_{i}^{*}}{\Lambda}\right)^{2}.$ (17)

Equating this to $kT/2$ gives a typical size of fluctuation
$\sqrt{\langle(R_{i}-R_{i}^{*})^{2}\rangle}=\Lambda\sqrt{t/(\pi l_{e})}$. For
$l_{e}/t=10$ this gives a relative fluctuation of $13\%$. We expect that the
fluctuations within the plane will be reasonably small.

In conclusion, ribbon polymers which have large bending rigidities in lateral
direction spontaneously form anisotropic condensates in poor solvents. The
size of those peculiar condensates quantized as a variation of the length of
polymers and solvent quality.

When the material is homogeneous, the surface tensions and elastic energies
are determined by geometry of the ribbon. In general, for homogeneous
materials, we may assume $\epsilon_{e}<\epsilon_{h}$ and
$\gamma_{e}\approx\gamma_{h}$, which is the case we discussed in this paper.
On the chemical scale, however, other anisotropic conditions are possible in
principle. The ribbon polymers with $\epsilon_{e}<\epsilon_{h}$ and
$w\gamma_{e}<t\gamma_{h}$[33], form a different type of anisotropic
condensates in poor solvents, that will be discussed somewhere else.

###### Acknowledgements.

The authors thank Jean-Louis Sikorav for discussions and comments on the
manuscript. This work was initiated with support by NTT Basic Research
Laboratories. Y.Y.S. wishes to acknowledge CEA, IPhT and Takushoku University,
RISE for financial support.

## References

  * [1] Klimenko S. M., Tikkchonenko T. I. Andreev V. M. J. Mol. Biol. 23 1967 523
  * [2] Grosberg A. Yu. Biophysics 24 1979 30
  * [3] Frank-Kamenetskii M. D., Anshelevich V. V. Lukashin A. V. Sov. Phys. Usp. 30 1987 317
  * [4] Marko J. F. Siggia E. D. Macromolecules 27 1994 981
  * [5] Odijk T. J. Chem. Phys. 105 1996 1270
  * [6] Marko J. F. Phys. Rev. E 55 1997 1758
  * [7] Haijun Z., Yang Z. Zhong-can O.-Y. Phys. Rev. Lett. 82 1999 4560
  * [8] Schlüter A.-D. Adv. Mater. 3 1991 282
  * [9] Nyrkova I. A., Semenov A. N., Joanny J.-F. Khokhlov A. R. J. Phys. II (Paris) 6 1996 1411
  * [10] Nyrkova I. A., Semenov A. N. Joanny J.-F. J. Phys. II (Paris) 7 1997 825
  * [11] Nyrkova I. A., Semenov A. N. Joanny J.-F. J. Phys. II (Paris) 7 1997 847
  * [12] Golestanian R. Liverpool T. B. Phys. Rev. E 62 2000 5488
  * [13] Landau L. D. Lifshitz E. M. Theory of Elasticity Pergamon, London 1970 77
  * [14] Bloomfield V. A. Biopolymers 31 1991 1471
  * [15] Bloomfield V. A. Biopolymers 44 1997 269
  * [16] Li A.-Z., Fan T-Y. Ding M. Science in China B 35 1992 169
  * [17] Ubbink J. Odijk T. Biophys. J. 68 1995 54
  * [18] Hud T., Downing K. H. Balhorn R. Proc. Natl. Acad. Sci. USA 92 1995 3581
  * [19] Park S. Y., Harries D. Gelbart W. M. Biophys. J. 75 1998 714
  * [20] Pereira G. G. Williams D. R. M. Biophys. J. 80 2001 161
  * [21] Takenaka Y., Yoshikawa K., Yoshikawa Y., Koyama Y. Kanbe T. J. Chem. Phys. 123 2005 014902
  * [22] Grosberg A. Yu. Zhestkov A. V. J. Biomol. Struct. Dyn. 3 1986 859
  * [23] Vasilevskaya V. V., Khokhlov A. R., Kidoaki S. Yoshikawa K. Biopolymers 41 1997 51
  * [24] Jary, D. Sikorav J.-L. Biochem. 38 1999 3223
  * [25] Odijk T. Phil. Trans. R. Soc. Lond. A 362 2004 1497
  * [26] Cooke I. R. Williams D. R. M. Physica A 339 2004 45
  * [27] Yoshikawa Y.,Yoshikawa K. Kanbe T. Biophys. Chem. 61 1996 93
  * [28] Schnurr B., MacKintosh F. C. Williams D. R. M. Europhys. Lett. 51 2000 279
  * [29] Pereira G. G. Williams D. R. M. ANZIAM J. 45 2004 C163
  * [30] Yoshinaga N., Yoshikawa K. Ohta T. Euro. Phys. J. E 17 2005 485
  * [31] Odijk T. Macromolecules 16 1983 1340
  * [32] Odijk T. Macromolecules 19 1986 2313
  * [33] Nelson J. C., Saven J. G., Moore J. S. Wolynes P. G. Science 277 1997 1793


# Null Hypothesis Test for Anomaly Detection

Jernej F. Kamenik<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39,
1000 Ljubljana, Slovenia  Faculty of Mathematics and Physics, University of
Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia  Manuel Szewc
<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39, 1000 Ljubljana,
Slovenia

###### Abstract

We extend the use of Classification Without Labels for anomaly detection with
a hypothesis test designed to exclude the background-only hypothesis. By
testing for statistical independence of the two discriminating dataset
regions, we are able to exclude the background-only hypothesis without relying
on fixed anomaly score cuts or extrapolations of background estimates between
regions. The method relies on the assumption of conditional independence of
anomaly score features and dataset regions, which can be ensured using
existing decorrelation techniques. As a benchmark example, we consider the LHC
Olympics dataset where we show that mutual information represents a suitable
test for statistical independence and our method exhibits excellent and robust
performance at different signal fractions even in presence of realistic
feature correlations.

##  I Introduction

The combination of increased experimental sensitivity and no clear leading
theoretical guide for how physics beyond the standard model would manifest in
current and future particle physics experiments has resulted in increased
development of anomaly detection techniques for collider applications, see
Ref. HEP ML Community  for a living review with a continuously updated list of
references. These techniques, which make use of state of the art unsupervised
and/or weakly supervised algorithms, have the advantage of being sensitive to
a large variety of signals at the expense of losing statistical power in
comparison to dedicated searches. However, appropriately quantifying said
sensitivity is still an open problem Shanahan  _et al._ (2022), with differing
proposals, see e.g. Ref. Park  _et al._ (2022). An especially pressing
question is how to evaluate the null hypothesis exclusion sensitivity of an
anomaly detection method. The current strategy is to perform cuts using the
anomalous score and extrapolate a background model from a control region. This
can be problematic for several reasons. First, the use of the anomalous score
itself to select events is not guaranteed to yield a robust method that
disentangles the underlying processes, see e.g. Ref. Kasieczka  _et al._
(2022) for a recent discussion of how ambiguities in the data representation
can lead to different notions of anomalous events which vary in their
discriminating power. Second, even if the anomaly score is an appropriate
event selection tool, the use of cuts, which in an unsupervised search cannot
be optimized on a targeted signal model, necessarily introduces a loss in
sensitivity by discarding possible signal events. Finally, the use of a
control region potentially introduces additional biases when assuming the
absence of signal in the control region and/or employing interpolation methods
such as the fit to a monotonic mass spectrum in a Bump Hunt.

In this work we aim to address some of the shortcomings outlined above. In
particular, we propose a null hypothesis statistical test for anomaly
detection which does not rely on fixed anomaly score cuts nor requires
background model extrapolations from control regions. We apply it to a
specific anomaly detection technique, Classification Without Labels (CWoLa)
introduced as a quark/gluon tagger in Ref. Metodiev  _et al._ (2017) and as an
anomaly detection technique in Refs. Collins  _et al._ (2018, 2019), and its
extension introduced in Ref. Benkendorfer  _et al._ (2021) incorporating
simulation assisted decorrelation of features. We show that by testing for
independence between the set of features used in the anomaly score, and those
used to define signal and control regions, we can obtain a p-value which
avoids false signal-detection and is robust in presence of slight correlations
between the two sets of features.

The work is structured as follows. In Section II we review CWoLa and introduce
the proposed statistical test. In Section III we apply our method to a LHC
Olympics benchmark to demonstrate its power and limitations. We conclude in
Section IV where we also discuss possible future extensions and improvements.
All the necessary code to reproduce our results is available at GitHub cod .

##  II Method

Introduced in Ref. Metodiev  _et al._ (2017), CWoLa is a weakly-supervised
technique for anomaly detection which aims to learn a monotonic function of
the Likelihood Ratio between Signal $S$ and Background $B$ processes for a set
of features of interest $\vec{x}$,
$\mathcal{L}_{S/B}(\vec{x})=p(\vec{x}|S)/p(\vec{x}|B)$, with the help of an
additional feature $y$ uncorrelated with $\vec{x}$. The latter variable, often
but not necessarily the invariant mass of the event, can be used to define two
regions of interest: the signal region $M_{1}$ and the control (or side-band)
region $M_{2}$, where the signal-to-background ratio is assumed to be higher
in $M_{1}$ than in $M_{2}$. A weakly-supervised algorithm, CWoLa trains a
classifier to distinguish between $M_{1}$ and $M_{2}$. The obtained output
function $s(\vec{x})$ can then be mapped to
$\mathcal{L}_{M_{1}/M_{2}}(\vec{x})$ through the likelihood ratio trick. The
orthogonality of $y$ and $\vec{x}$ guarantees that
$\mathcal{L}_{M_{1}/M_{2}}(\vec{x})$ is a monotonous function of
$\mathcal{L}_{S/B}(\vec{x})$ and thus possesses in principle optimal
statistical power.

Usual applications of CWoLa use the learned optimal classifier $s(\vec{x})$ to
select events of interest and assign a certain significance to the difference
in selected events in $M_{1}$ and $M_{2}$. The difference in the resulting
selection efficiencies $\epsilon_{M_{1,2}}$ is a smoking-gun for the presence
of signal in $M_{1}$ (and also $M_{2}$). However, this is only true in the
limit of infinite statistics. In a realistic setting where the dataset is
finite, quantifying the degree to which the difference in efficiencies relates
to the presence of signal is non-trivial. One common strategy is to assume
that there is no signal in $M_{2}$ and assess the agreement between the
selected events in $M_{1}$ and a background extrapolation from $M_{2}$.

Our method constitutes an alternative to assess how the learned output
$s(\vec{x})$ encodes differences between $M_{1}$ and $M_{2}$ caused by the
presence of a signal. To introduce it, we focus on the density estimation
framing of CWoLa, which clearly defines a background-only or null hypothesis.
At its heart, CWoLa is a mixture model where $\vec{x}$ and $y$ are assumed to
be conditionally independent given the process label $z=\\{S,B\\}$. After
defining $M_{1}$ and $M_{2}$ using $y$, the trained classifier output is a
function $s(\vec{x})$ that inherits the conditional independence with respect
to $y$. The statistical model can be explicitly written as

$p(s(\vec{x}),y|\pi)=(1-\pi)\,p(s(\vec{x})|B)p(y|B)+\pi\,p(s(\vec{x})|S)p(y|S)\,,$
(1)

where $\pi$ is the signal probability. The background-only hypothesis is
explicitly written as $p(s(\vec{x}),y|\pi=0)$ and corresponds to the case
where the observed data shows independence between $s(\vec{x})$ and $y$. This
is the key observation for our strategy. For a given measured dataset of pairs
$\\{s(\vec{x_{i}}),y_{i}\\}$, one can assess whether they are statistically
independent. If statistical independence is ruled out, the background-only
hypothesis is ruled out, provided conditional independence holds. Conversely,
if statistical independence cannot be ruled out, one has a clear statement
about the incapability of CWoLa to discern whether any difference between
$M_{1}$ and $M_{2}$ originates from the presence of a signal or is due to
statistical fluctuations in the data.

Several tests of statistical independence exist for both discrete and
continuous distributions, including mutual information Bishop (2006),
Hoeffding’s D independence test Hoeffding (1948) and distance correlation
Székely  _et al._ (2007). For simplicity, in the present work we focus on the
use of the estimated mutual information (MI) $I$ of the measured probability
distribution. MI encodes the exact property we want to test as it measures the
difference between the joint distribution and the marginals:

$\displaystyle I(s,y)$ $\displaystyle=D_{\text{KL}}(p(s,y)||p(s)p(y))$ (2)
$\displaystyle=\int\,ds\,dy\ p(s,y)\log\frac{p(s,y)}{p(s)p(y)}\,,$ (3)

where $D_{\text{KL}}(p,q)$ is the Kullback-Leibler divergence between two
probability distributions, capturing how much information is lost when
approximating the distribution $p$ with the distribution $q$. The MI thus
captures how well one can approximate the joint distribution by the product of
its marginals and it is trivial to show that it vanishes for independent
variables. Conditional Independence can then be expressed as a vanishing MI
conditioned on a given process

$I(s,y|z)=\int\,ds\,dy\ p(s,y|z)\log\frac{p(s,y|z)}{p(s|z)p(y|z)}=0\,.$ (4)

On the other hand, for the full dataset the possible mixture between the two
processes encoded in $\pi\in[0,1]$ results in

$I(s,y)\geq 0\,,$ (5)

with the equality achieved when there is only one process or the two processes
have the same probability distributions.

A very nice feature of the MI is that it has well behaved asymptotic
properties in the limit of small MI and large sample size Goebel  _et al._
(2005). Thus, we can estimate it from the measured sample of $N$ events and
obtain the p-value of said estimator $\hat{I}(s,y)$ under the null hypothesis
$I(s,y)=0$. Assuming a two dimensional binning of $(s,y)$ with $d_{s}$ and
$d_{y}$ the number of chosen bins per variable, the estimator $\hat{I}$ is a
random variable that under the null hypothesis $I=0$ follows a Gamma
distribution with shape parameter $\frac{(d_{s}-1)(d_{y}-1)}{2}$ and scale
parameter $N$.

To estimate $\hat{I}$ we need to estimate $\hat{p}(s,y)$, with $\hat{p}(s)$
and $\hat{p}(y)$ obtained by marginalizing. We estimate $\hat{p}(s,y)$ through
two-dimensional histogram event counts with the aforementioned $d_{s}$ and
$d_{y}$ chosen bins. Because we are dealing with continuous variables, the use
of binning introduces additional hyperparameters. In this work we bin $s$ and
$y$ in such a way that each bin has a relative statistical uncertainty equal
or lower than $1\%$. Other criteria for statistical independence that deal
explicitly with continuous variables such as Hoeffding’s D independence test
or distance correlation could be used to avoid the introduction of binning at
the expense of increased computational cost. We choose MI as it is straight
forward to implement with a general signal-blind binning criteria and it
suffices to establish the relevance of the strategy detailed in this work.

We emphasize that the role of CWoLa is to provide a one-dimensional observable
$s(\vec{x})$ which can then be combined with $y$ to test for statistical
independence. Once $s(\vec{x})$ is obtained, the rest of the test relies only
on data without the need to introduce additional cuts or labels. If testing
for statistical dependence between $\vec{x}$ and $y$ directly was feasible,
then one would not need to introduce any learnable function. However, this is
often not the case. One in general needs a high-dimensional $\vec{x}$ to
ensure discriminative power between possible signals and the background, which
is in turn converted by our method into statistical power to exclude
statistical independence. On the other hand working directly with a high-
dimensional set of features renders any statistical test problematic either
due to the test being designed for two variables, as is the case for
Hoeffding’s D independence test and distance correlation, or due to the
necessary density estimation suffering from the course of dimensionality as is
the case for the mutual information test presented in this work.

The method relies on the assumption of conditional independence between
$\vec{x}$ and $y$. In a realistic application this is not ensured, specially
when considering highly-discriminative variables between the background and
potential signals. The presence of correlation between $\vec{x}$ and $y$ will
result in non-null MI for each process separately. Thus, the p-value obtained
from the MI estimation will be merely testing for conditional independence,
not the presence of a single process. In other words, the null hypothesis
ceases to be equal to the background-only hypothesis. This challenge is
already present in current implementations of CWoLa, with correlations
resulting in loss of classification power.

One possible strategy introduced in Ref. Benkendorfer  _et al._ (2021) is to
ensure that $s(\vec{x})$ is agnostic to the correlation between $\vec{x}$ and
$y$ through the addition of a simulated background dataset during the training
stage. In this approach, named Simulation Assisted Classification Without
Labels (SA-CWoLa), the loss function is modified with an additional term that
incorporates the simulation dataset. Following Ref. Benkendorfer  _et al._
(2021), we define the loss function as

$\displaystyle\mathcal{L}_{\text{SA-CWoLa}}[s]$ $\displaystyle=$
$\displaystyle-\left(\sum_{\vec{x}_{n}\in M^{\text{data}}_{1}}\log
s(\vec{x}_{n})+\sum_{\vec{x}_{n}\in
M^{\text{data}}_{2}}\log\left(1-s(\vec{x}_{n})\right)\right)$ (6)
$\displaystyle-$ $\displaystyle\lambda\left(\sum_{\vec{x}_{n}\in
M^{\text{sim.}}_{1}}\log\left(1-s(\vec{x}_{n})\right)+\sum_{\vec{x}_{n}\in
M^{\text{sim.}}_{2}}\log s(\vec{x}_{n})\right)\,,$

that inverts the labelling in the simulation so as to penalize learning
background differences between $M_{1}$ and $M_{2}$, with $\lambda$ the hyper-
parameter that controls the relative importance of said penalization 111We
always reweigh the events from both data and simulation during training in
such a way that the each of the four subset of events
$\\{M^{\text{data}}_{1},M^{\text{data}}_{2},M^{\text{sim.}}_{1},M^{\text{sim.}}_{2}\\}$
has the same total weight.. Note that the specific choice of the loss function
is not relevant as long as it ensures proper decorrelation. Similarly, the
simulated dataset does not need to be perfect, it only needs to encode
accurately enough the correlation between $\vec{x}$ and $y$ for the background
process. Learning to ignore the correlations in the background guarantees that
excluding the null hypothesis $I(s,y)=0$ corresponds to excluding the
background-only hypothesis $p(s,y)=p(s,y|\pi=0)$ and not merely excluding
conditional independence $p(s,y|B)=p(s|B)p(y|B)$. The main drawback of
introducing decorrelation is that the learned function $s(\vec{x})$ ceases to
be optimal and looses classification power. Thus, one should tune $\lambda$
with a given criteria that balances learning to decorrelate between $\vec{x}$
and $y$ for $B$ and learning to distinguish between $B$ and $S$ through
discriminating between $M_{1}$ and $M_{2}$.

In this proof-of-principle, we are satisfied with presenting results for fixed
$\lambda$ that is large enough to ensure decorrelation in the simulation
sample and at the same time small enough so that $s(\vec{x})$ is sensitive to
the presence of signal in the measured sample and improves over the naive
significance estimation $S/\sqrt{B}$. To verify that decorrelation is
enforced, we follow Ref. Benkendorfer  _et al._ (2021) and compute the Area-
Under-Curve (AUC) for the $s(\vec{x})$ classifier for the $M_{1}$ and $M_{2}$
samples in the simulation dataset. If the AUC is approximately $0.5$, we have
a classifier that is not better than a random classifier and thus it has
learned to ignore any possible correlations between $\vec{x}$ and $y$. We
emphasize that the sole purpose of the simulation dataset is to ensure
decorrelation, and we never compare data to simulation to obtain a
significance after training. This makes the test more robust against
background mismodelling than other unsupervised methods for anomaly detection
which avoid the use of anomaly cuts at the expense of performing data-to-
simulation hypothesis tests such as Refs. D’Agnolo and Wulzer (2019); D’Agnolo
_et al._ (2021); d’Agnolo  _et al._ (2022); Letizia  _et al._ (2022);
Krzyżańska and Nachman (2022); Chakravarti  _et al._ (2021). This is only
possible because we assume that the simulations are precise enough to capture
qualitative correlations between features in data. The simulator precision
needed for decorrelation to be effective is considerably less than what is
needed for a full multivariate comparison with measurements which often
suffers both from systematic biases in the simulations and from the
computational cost of achieving a given statistical precision.

##  III Application: LHC Olympics

In order to demonstrate its power, we apply our method to the LHC Olympics R&D
labelled dataset Kasieczka  _et al._ (2021). The dataset is comprised of dijet
events from two different sources: SM quantum chromodynamics (QCD) processes
(background $B$), and the production of a hypothetical new resonance
$W^{\prime}$ with mass $m_{W^{\prime}}=3.5$ TeV, decaying to two intermediate
particles $X$ and $Y$ with masses $m_{X}=500$ GeV and $m_{Y}=100$ GeV, which
in turn both decay promptly to pairs of quarks producing two large-radius jets
with a two-prong substructure (signal $S$). Our variable of interest $y$ is
the reconstructed dijet invariant mass $m_{jj}$ of the two hardest (in
$p_{T}$) jets in the event. Mimicking the selection criteria of Ref.
Benkendorfer  _et al._ (2021), the selected events have a reconstructed dijet
mass $m_{jj}\in[3.1,3.9]\,\mbox{\rm TeV}$. To perform CWoLa Metodiev  _et al._
(2017), we define two orthogonal regions
$M_{1}\equiv\\{m_{jj}\in[3.3,3.7]\,\mbox{\rm TeV}\\}$ and
$M_{2}\equiv\\{m_{jj}\in[3.1,3.3]\,\mbox{\rm TeV}\cup[3.7,3.9]\,\mbox{\rm
TeV}\\}$. In the following, we refer to $S$ and $B$ as the total number of
Signal and Background events in $M_{1}\cup M_{2}$.

For anomaly score input features $\vec{x}$, we choose a set of variables based
on the invariant masses and the first N-subjettiness ratios Thaler and Van
Tilburg (2011, 2012) of the two selected jets. Ordering the jets by mass, with
$j=1$ being the heavier jet, our variables are

$\vec{x}=\\{m_{1}-m_{2},m_{2},\tau_{21,1},\tau_{21,2}\\}\,.$ (7)

The correlation between $\vec{x}$ and $m_{jj}$ is mostly concentrated in the
correlations between $\\{m_{1},m_{2}\\}$ and $m_{jj}$. To illustrate how
important they are, we define

$\Delta^{ab}_{c}(m^{\text{bin}}_{jj})=\frac{\mathbb{E}[m_{c}|m_{jj}\in
m^{\text{bin}}_{jj},a]-\mathbb{E}[m_{c}|b]}{\mathbb{E}[m_{c}|b]}\,,$ (8)

where $a,b\in\\{B,S\\}$ and $c\in\\{1,2\\}$.
$\Delta^{ab}_{c}(m^{\text{bin}}_{jj})$ represents the relative difference
between the average of $m_{c}$ in a given $m_{jj}$ bin $m^{\text{bin}}_{jj}$
for process $a$ and the average of the same observable over the whole $m_{jj}$
range for process $b$. If $a=b$, this observable conveys the presence of a
correlation between $m_{c}$ and $m_{jj}$ for a given process. For $a\neq b$,
this observable conveys the difference in the $m_{c}$ probability
distributions for the two processes and its dependence on $m_{jj}$. By
comparing $\Delta^{BB}_{c}$ with $(S/B)\Delta^{SB}_{c}$ we can check whether
the correlations between $m_{c}$ and $m_{jj}$ can obscure the differences
between signal and background that CWoLa aims to learn. The prefactor $S/B$ is
introduced to account for the fact that the data contains less signal than
background and originates from comparing the full data distribution to the
background-only distribution and separating the signal and the background
contributions:

$\Delta^{S+B,B}_{c}=\frac{B}{S+B}\Delta^{BB}_{c}+\frac{S}{S+B}\Delta^{SB}_{c}\approx\Delta^{BB}_{c}+(S/B)\Delta^{SB}_{c}\,.$

We show in Fig. 1 the resulting distributions for $m_{1}$ and $m_{2}$, with
the largest $S/B$ considered in this work, $S/B=0.01$. We observe how the
correlation between features for the background process, as evidenced by the
monotonic increase of $\Delta^{BB}_{c}$ from negative to positive values
towards larger $m_{jj}$, is sizable and crucially more pronounced compared to
the $S/B$ weighted difference between $S$ and $B$ as traced by
$\Delta^{SB}_{c}$. In other words, the correlation between $m_{c}$ and
$m_{jj}$ in the background can easily mask the presence of a small signal. For
lower $S/B$, the correlations become even more dominant and can lead to a
strongly biased anomaly score. In addition, in our approach the correlations
will also induce statistical dependence between $s(\vec{x})$ and $M_{1,2}$
even in absence of a signal and thus jeopardize the validity of the null-
hypothesis test.

---
Figure 1: Distribution of $\Delta^{ab}_{c}$ defined in Eq. (8) comparing
correlations between $\\{m_{1},m_{2}\\}$ and $m_{jj}$ to the differences
between $S$ and $B$ for the largest $S/B$ considered. Vertical error bars
signify statistical uncertainties, while horizontal bars indicate $m_{jj}$ bin
width. Vertical dashed lines denote boundaries between signal ($M_{1}$) and
side-band ($M_{2}$) regions. See text for details.

In order to address this crucial issue, we follow Ref. Benkendorfer  _et al._
(2021) and incorporate a set of simulated events into the training stage. As
simulation, we consider the background provided in the labelled version of the
Black Box 1 (BB1) dataset. We use BB1 as simulation to take advantage of the
larger signal sample provided in the R&D dataset. Using the loss function
defined in Eq. (6), the classifier is then trained to distinguish $M_{1}$ and
$M_{2}$ in data but not in simulation, obtaining a $s(\vec{x})$ that is
agnostic to correlations between $\vec{x}$ and $y$ for QCD. As a classifier,
we use a very similar set-up as in Ref. Benkendorfer  _et al._ (2021): we
train a Neural Network composed of three hidden layers with 64 nodes each and
ReLU activation function with a sigmoid function applied to the output to
ensure $s(\vec{x})\in[0,1]$ for 20 epochs using the ADAM Kingma and Ba (2014)
optimizer. The Neural Network is implemented in PyTorch Paszke  _et al._
(2019). We have trained our classifier using $k$-fold cross-validation with
$k=10$ to avoid overfitting by ensuring that every $s(\vec{x}_{n})$ is
obtained by combining the data point $\vec{x}_{n}$ with a classifier which has
not seen $\vec{x}_{n}$ during training. We also perform several random weight
initializations to ensure better convergence. At the end of training, we
evaluate the AUC score between the $M_{1}$ and $M_{2}$ simulated samples to
ensure that decorrelation is achieved.

---
Figure 2: Anomaly score $s(\vec{x})$ probability density function (PDF) after
training for different event labellings. Each row corresponds to a given $S/B$
and each column to a given $\lambda$. Non-uniform binning ensures that when
considering the full dataset each bin has a relative statistical uncertainty
equal or lower than $1\%$, resulting in between 15 and 25 bins per plot. The
data is shown both labelled according to the (observable) values of $m_{jj}$
(defining $M_{1}$ and $M_{2}$) as well as according to the (unobservable)
truth labels (background and signal).

We show in Fig. 2 the learned $s(\vec{x})$ for different $S/B$ with $B=250$k
and $\lambda=\\{0.0,1.0\\}$. In each plot, the binning is chosen in such a way
that each bin has a relative statistical uncertainty lower than or equal to
$1\%$. This choice ensures good performance of the density estimation needed
for the hypothesis test. We can appreciate how as $S/B$ increases,
$s(\vec{x})$ goes from being mostly centered around $s=0.5$ to yielding higher
$s$ values, indicating improved learning of signal features. However, the
binning choice obscures somewhat how much the background and signal are
separated (within the highest $s$ bin), as the whole signal is grouped
together with the necessary background events to obtain a $1\%$ statistical
uncertainty. In absence of correlation mitigation (for $\lambda=0$), anomaly
score bias causes clearly unbalanced classification of events in $M_{1}$ and
$M_{2}$ even in absence of any signal. The introduction of $\lambda=1$ causes
the events to be even more centered around $0.5$, specially for low $S/B$.
However, more importantly, $\lambda\geq 0$ forces the training to ignore
possible correlations between $\vec{x}$ and $y$ for (simulated) background and
thus approaches a random classifier for $S/B\to 0$.

---
Figure 3: PDFs of estimated mutual information, its numerical distribution
under the null hypothesis estimated through resampling, and its asymptotic
distribution under the null hypothesis, for different considered benchmark
datasets. Each plot corresponds to a different choice of $S/B$ with
$\lambda=1$. We combine the estimated $\hat{I}$ with their asymptotic
distribution to obtain the resulting p-values. When $S/B>0$, we also compute
the discovery significance $Z=\Phi^{-1}(1-p)$, where $\Phi$ is the unit
Gaussian cumulative distribution function, and compare it to
$Z_{0}=S/\sqrt{B}$.

As expected, the impact of $\lambda$ is even more pronounced when testing for
statistical independence. For each dataset of $\\{s,m_{jj}\\}$ values, we
estimate mutual information and obtain the p-value associated with the null
hypothesis as detailed in Section II. We show in Fig. 3 the resulting
estimated $\hat{I}_{\mathrm{data}}$ and their corresponding p-values. We also
ran a series of pseudo-experiments to verify that the asymptotic limit is
appropriate. We only present results for $\lambda=1$ because for $\lambda=0$
we are able exclude $I(s,y)=0$ for all $S/B$ with very high confidence
(p-value$\,<10^{-14}$). This, as detailed in Section II, is because we are
excluding conditional independence in the background process due to
correlations between $\vec{x}$ and $y$. This is specially important for
$S/B=0$, where the effect of correlations can mislead CWoLa to falsely exclude
the background-only hypothesis.

We observe that for $\lambda=1$ the proposed test has the required behavior:
for $S/B=0$ the test yields results consistent with statistical independence,
while an increase of $S/B$ leads to an increasingly strong exclusion of the
null hypothesis. The use of SA-CWoLa thus ensures that we can identify the
null hypothesis with the background-only hypothesis. For $S/B>0$, we also
compute the discovery significance $Z=\Phi^{-1}(1-p)$, where $\Phi$ is the
unit Gaussian cumulative distribution function, and compare it to the naive
counting significance $Z_{0}=S/\sqrt{B}$. We observe how our method presents
an increased discovery significance even compared to the case of perfect (up
to statistical fluctuation) knowledge of background yields.

Overall, Fig. 3 shows how $\hat{I}$ can be used to infer whether the data
presents a deviation from the null hypothesis, defined as the case where a
single process (or an $m_{jj}$ independent mixture of processes) is present
for which $s$ and $m_{jj}$ are independent. This is the analogous to the
p-value obtained using the Bump Hunt in a “traditional” implementation of
CWoLa. However, contrary to existing approaches, here there is no selection
cut to be optimized, and no extrapolation of the background into the signal
region is required.

In principle the method could learn the likelihood-ratio between signal and
background, which is the optimal test-statistic. However, the presence of the
decorrelation term in Eq. (6) reduces the optimality of $s$ and consequently
$\hat{I}$. This can be seen by the decrease in $Z/Z_{0}$ as $S/B$ increases.
When $S/B$ is large enough, $s(\vec{x})$ will ignore the small correlations in
each individual process even for $\lambda=0.0$. A non-null $\lambda$ will thus
only worsen the performance of the algorithm. However, as we are interested in
low $S/B$ cases where anomaly detection is useful, a more conservative
approach which is robust to correlations even when there is no signal present
is warranted.

In the previous paragraphs, we have shown how the proposed method yields an
appropriate test statistic which is different from existing approaches. In
Table 1 we provide a significance comparison of the proposed method to two
traditional implementation of CWoLa which we denote as “Anomaly cuts” and
“Bump Hunt”. The former, implemented e.g. in Ref. Finke  _et al._ (2022),
assumes that no signal is present in the side-band region $M_{2}$ and
estimates the total number of background events in the signal region $M_{1}$
through the use of cuts on the anomaly score. For a fixed efficiency in the
side-band region $\epsilon_{2}$, the estimated background event yield in
signal region is $\epsilon_{2}N_{1}$. If the measured efficiency in the signal
region $\epsilon_{1}$ is larger than $\epsilon_{2}$, then there is an excess
of events with a significance of

$Z=\begin{cases}\frac{(\epsilon_{1}-\epsilon_{2})N_{1}}{\sqrt{\epsilon_{2}(N_{1}+N_{2})}},&\quad\text{if
}\epsilon_{1}\geq\epsilon_{2}\,,\\\\[8.61108pt] 0,&\quad\text{if
}\epsilon_{1}<\epsilon_{2}\,.\end{cases}$

Similarly, the Bump Hunt method also assumes no signal is present in $M_{2}$
but estimates the background event yield in $M_{1}$ differently. In this
approach, the background $m_{jj}$ distribution is explicitly modelled and a
Profile Likelihood Ratio fit of the number of signal events in $M_{1}$ is
performed. We follow Ref. Benkendorfer  _et al._ (2021) and model the
background distribution for $m_{jj}\in[3.1,3.9]$ TeV as

$\frac{d\sigma}{dm_{jj}}=\frac{p_{0}\left(1-\frac{m_{jj}}{\sqrt{s}}\right)^{p_{1}}}{\left(\frac{m_{jj}}{\sqrt{s}}\right)^{p_{2}+p_{3}\log\frac{m_{jj}}{\sqrt{s}}}}\,,$

where $\sqrt{s}$ is the center-of-mass energy and $p_{i}$ are parameters to be
fitted from the $m_{jj}$ distribution with the signal region masked. Once the
fit is performed, the expected number of background events in $M_{1}$ $b$ is
estimated from the integral of the background distribution over the signal
region.We define the Likelihood function in the signal region as

$\mathcal{L}(s,\theta)=\mathcal{P}(N_{1}|s+b+\theta)\mathcal{N}(\theta|0,\sigma)\,,$

where $\mathcal{P}$ is the Poisson probability mass function of measuring
$N_{1}$ events, $\mathcal{N}$ is the Normal probability density, $\theta$ is a
nuisance parameter for background mismodelling and $\sigma$ is the background
yield error propagated from the $p_{i}$ fit. From this Likelihood we build the
usual test statistic

$q_{0}=\begin{cases}-2\text{ Ln
}\mathcal{L}(0,\hat{\hat{\theta}})/\mathcal{L}(\hat{s},\hat{\theta}),&\quad\text{if
}\hat{s}\geq 0\,,\\\\[8.61108pt] 0,&\quad\text{if }\hat{s}<0\,,\end{cases}$

where $\hat{s},\hat{\theta}$ are the maximum likelihood estimates of $s$ and
$\theta$ and $\hat{\hat{\theta}}$ is the maximum likelihood estimate of
$\theta$ when keeping $s$ fixed to 0. The resulting significance is
$Z=\sqrt{q_{0}}$. We implement the Bump Hunt by itself and in conjunction with
the use of cuts in the anomaly score to enhance the efficiency as would be
done in a resonance search.

From Table 1 we observe how in every case the introduction of $\lambda>0$
reduces the significance. However, it does not imply resilience against
spurious signals for all strategies. Both traditional methods are highly
dependent on the arbitrary $\epsilon_{2}$ choice, showing the appearance of
spurious significance at $S/B=0$ for certain values. In general, we observe
that our method is better suited for smaller $S/B$ than existing methods. This
is mainly because it does not discard potential signal events. The Bump Hunt
with no cuts also does this, but its significance is lower for non-null $S/B$
(for the extreme $S/B=0.01$, its significance is even lower than $S/\sqrt{B}$
due to the presence of the nuisance parameter). From this comparison, we
assess that without clear criteria for an optimal anomaly cut
($\epsilon_{2}$), the mutual information-based method performs better for low
to null $S/B$ whereas the cut and count based methods outperform the mutual
information test for larger $S/B$.

Another benefit of our model compared to traditional searches is scalability
with sample size. If decorrelation is ensured, the larger the sample size the
more powerful the method. This is certainly not true for the Bump Hunt, where
the background modelling is an inherent approximation which necessarily
introduces bias for large enough sample sizes. Conversely, for smaller
datasets our method takes advantage of the full dataset in a better way than
through the use of fixed anomaly cuts. The asymptotic approximation for the
mutual information CDF, which is vital for the test, has been shown
numerically Goebel  _et al._ (2005) to be valid for $N>50$ events and true
mutual information $I\leq 0.14$, conditions that we expect to always be
satisfied in a realistic anomaly search at the LHC.

Significance |  $S/B$ = 0.0 |  $S/B$ = 0.0025 |  $S/B$ = 0.005 |  $S/B$ = 0.01
---|---|---|---|---
$S/\sqrt{B}$ | 0.0 | 1.29 | 2.55 | 6.40
Mutual Info $\lambda=0.0$ | 6.40 | 7.04 | 7.58 | 14.1
Mutual Info $\lambda=1.0$ | 0.70 | 3.03 | 5.33 | 13.0
Anomaly cuts $\epsilon_{2}=0.1$, $\lambda=0.0$ | 3.35 | 4.78 | 6.27 | 11.6
Anomaly cuts $\epsilon_{2}=0.1$, $\lambda=1.0$ | 2.48 | 2.26 | 4.49 | 10.0
Anomaly cuts $\epsilon_{2}=0.01$, $\lambda=0.0$ | 2.26 | 4.62 | 10.1 | 27.0
Anomaly cuts $\epsilon_{2}=0.01$, $\lambda=1.0$ | 0.55 | 1.66 | 10.7 | 27.1
Anomaly cuts $\epsilon_{2}=0.001$, $\lambda=0.0$ | 1.39 | 10.3 | 17.9 | 34.2
Anomaly cuts $\epsilon_{2}=0.001$, $\lambda=1.0$ | 0. | 0.57 | 13.6 | 37.0
Bump Hunt | 0.95 | 1.97 | 2.74 | 5.30
Bump Hunt $\epsilon_{2}=0.1$, $\lambda=0.0$ | 6.41 | 9.26 | 10.92 | 19.5
Bump Hunt $\epsilon_{2}=0.1$, $\lambda=1.0$ | 3.81 | 4.35 | 6.93 | 16.0
Bump Hunt $\epsilon_{2}=0.01$, $\lambda=0.0$ | 4.77 | 6.96 | 14.2 | 34.7
Bump Hunt $\epsilon_{2}=0.01$, $\lambda=1.0$ | 0.97 | 2.53 | 14.0 | 35.0
Bump Hunt $\epsilon_{2}=0.001$, $\lambda=0.0$ | 2.98 | 12.3 | 20.0 | 35.8
Bump Hunt $\epsilon_{2}=0.001$, $\lambda=1.0$ | 0.29 | 1.60 | 15.9 | 38.7
Table 1: Significances obtained with different strategies for different $S/B$
ratios, see text for details.

##  IV Discussions and outlook

In this work, we have presented a novel strategy to quantify the sensitivity
of a specific anomaly detection technique, Simulation-Assisted Classification
Without Labels, by testing for statistical independence of the learned
$\\{s(\vec{x}),y\\}$ samples. We have shown that as long as one can rely on
SA-CWoLa to enforce conditional independence of the background processes
$p(s,y|B)=p(s|B)p(y|B)$, the null hypothesis of statistical independence is
equivalent to the background-only hypothesis. Thus, testing for statistical
independence in the observed data corresponds to testing for the background-
only hypothesis.

As a proof of principle, we have considered mutual information as a test
statistic. MI has a known asymptotic distribution under the null hypothesis
for binned data and has low computational cost. We have tested our method with
LHC Olympics datasets and have shown that the test statistic yields the
expected behavior. Most importantly our proposed test statistic provides a
clear statement on the presence of signal, i.e. is capable of correctly
yielding a no-signal response. This opens the door to testing for new physics
in LHC datasets without the need for anomaly score cuts, as well as reducing
the need for accurate background modelling.

Possible extensions of the present work could consider other tests for
statistical independence such as Hoeffding’s D independence test or distance
correlation, which can be applied on unbinned $s(\vec{x})$ and $y$, at the
expense of increased computational cost. Similarly, other methods for anomaly
score training and decorrelation of features could be explored. Employing the
most suitable classification and decorrelation methods can be model and
dataset dependent, and has not been the main focus of this work.

Another possibility is to assume that no signal populates the side-bands and
identify $p(\hat{s}|M_{2})=p(\hat{s}|z=0)$. This opens the door for an optimal
analysis since one can now perform a template fit in the signal region Nachman
(2020). However, it also potentially introduces additional uncertainties
and/or biases due to background modelling that the present method avoids. We
leave a more complete study in this direction for future work.

Regarding other physics applications, we emphasize that by dispensing with the
need for explicit functional background modelling, our test is especially
useful for anomaly detection applications that do not search for predetermined
(modelled) signal shapes Shanahan  _et al._ (2022), such as invariant mass
resonances as in e.g. Refs. Alvarez  _et al._ (2020); Finke  _et al._ (2022).
However, it could also be applied in supervised searches as an additional
cross-check to control bias due to modelling at the expense of loss of
optimality.

##  ◆ ◆ ◆ Acknowledgments ◆ ◆ ◆

The authors acknowledge the financial support from the Slovenian Research
Agency (grant No. J1-3013 and research core funding No. P1-0035). MS is
grateful to the Mainz Institute for Theoretical Physics (MITP) of the Cluster
of Excellence PRISMA+ (Project ID 39083149) and to GGI for their hospitality,
support and useful discussions.

## References

  * (1) HEP ML Community, “A Living Review of Machine Learning for Particle Physics,” .
  * Shanahan  _et al._ (2022) P. Shanahan _et al._ , (2022), arXiv:2209.07559 [physics.comp-ph] .
  * Park  _et al._ (2022) S. E. Park, P. Harris, and B. Ostdiek, (2022), arXiv:2208.05484 [hep-ph] .
  * Kasieczka  _et al._ (2022) G. Kasieczka, R. Mastandrea, V. Mikuni, B. Nachman, M. Pettee, and D. Shih, (2022), arXiv:2209.06225 [hep-ph] .
  * Metodiev  _et al._ (2017) E. M. Metodiev, B. Nachman, and J. Thaler, JHEP 10, 174 (2017), arXiv:1708.02949 [hep-ph] .
  * Collins  _et al._ (2018) J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. Lett. 121, 241803 (2018), arXiv:1805.02664 [hep-ph] .
  * Collins  _et al._ (2019) J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. D 99, 014038 (2019), arXiv:1902.02634 [hep-ph] .
  * Benkendorfer  _et al._ (2021) K. Benkendorfer, L. L. Pottier, and B. Nachman, Phys. Rev. D 104, 035003 (2021), arXiv:2009.02205 [hep-ph] .
  * (9) “Null Hypothesis Test for Anomaly Detection,” https://github.com/ManuelSzewc/Null_Hypothesis_Test_for_Anomaly_Detection.
  * Bishop (2006) C. M. Bishop,  _Pattern recognition and machine learning_, Information science and statistics (Springer, New York, NY, 2006) softcover published in 2016\.
  * Hoeffding (1948) W. Hoeffding, The Annals of Mathematical Statistics 19, 546 (1948).
  * Székely  _et al._ (2007) G. J. Székely, M. L. Rizzo, and N. K. Bakirov, The Annals of Statistics 35, 2769 (2007).
  * Goebel  _et al._ (2005) B. Goebel, Z. Dawy, J. Hagenauer, and J. Mueller (2005) pp. 1102 – 1106 Vol. 2.
  * D’Agnolo and Wulzer (2019) R. T. D’Agnolo and A. Wulzer, Phys. Rev. D 99, 015014 (2019), arXiv:1806.02350 [hep-ph] .
  * D’Agnolo  _et al._ (2021) R. T. D’Agnolo, G. Grosso, M. Pierini, A. Wulzer, and M. Zanetti, Eur. Phys. J. C 81, 89 (2021), arXiv:1912.12155 [hep-ph] .
  * d’Agnolo  _et al._ (2022) R. T. d’Agnolo, G. Grosso, M. Pierini, A. Wulzer, and M. Zanetti, Eur. Phys. J. C 82, 275 (2022), arXiv:2111.13633 [hep-ph] .
  * Letizia  _et al._ (2022) M. Letizia, G. Losapio, M. Rando, G. Grosso, A. Wulzer, M. Pierini, M. Zanetti, and L. Rosasco, Eur. Phys. J. C 82, 879 (2022), arXiv:2204.02317 [hep-ph] .
  * Krzyżańska and Nachman (2022) K. Krzyżańska and B. Nachman, (2022), arXiv:2203.09601 [hep-ph] .
  * Chakravarti  _et al._ (2021) P. Chakravarti, M. Kuusela, J. Lei, and L. Wasserman, (2021), arXiv:2102.07679 [stat.AP] .
  * Kasieczka  _et al._ (2021) G. Kasieczka _et al._ , (2021), arXiv:2101.08320 [hep-ph] .
  * Thaler and Van Tilburg (2011) J. Thaler and K. Van Tilburg, JHEP 03, 015 (2011), arXiv:1011.2268 [hep-ph] .
  * Thaler and Van Tilburg (2012) J. Thaler and K. Van Tilburg, JHEP 02, 093 (2012), arXiv:1108.2701 [hep-ph] .
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,”  (2014).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Finke  _et al._ (2022) T. Finke, M. Krämer, M. Lipp, and A. Mück, JHEP 08, 015 (2022), arXiv:2204.11889 [hep-ph] .
  * Nachman (2020) B. Nachman, SciPost Phys. 8, 090 (2020), arXiv:1909.03081 [hep-ph] .
  * Alvarez  _et al._ (2020) E. Alvarez, F. Lamagna, and M. Szewc, JHEP 01, 049 (2020), arXiv:1911.09699 [hep-ph] .