# Spontaneous Hopf fibration in the two-Higgs-doublet model

R. A. Battye<EMAIL_ADDRESS>Jodrell Bank Centre for
Astrophysics, Department of Physics and Astronomy, University of Manchester,
Manchester, M13 9PL, U.K.  S. J. Cotterill<EMAIL_ADDRESS>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Manchester, M13 9PL, U.K.

###### Abstract

We show that energetic considerations enforce a Hopf fibration of the Standard
Model topology within the 2HDM whose potential has either an $SO(3)$ or $U(1)$
Higgs-family symmetry. This can lead to monopole and vortex solutions. We find
these solutions, characterise their basic properties and demonstrate the
nature of the fibration along with the connection to Nambu’s monopole
solution. We point out that breaking of the $U(1)_{\rm EM}$ in the core of the
defect can be a feature which leads to a non-zero photon mass there.

Two Higgs Doublet Model, Topological Defects

Introduction: The vacuum manifold of the Standard Model of Particle Physics is
a 3-sphere $({\cal M}=S^{3})$ implying that there are no stable topological
configurations in 3D since $\pi_{2}({\cal M})=I$, where $\pi_{n}({\cal M})$ is
the $n$-th homotopy group of the manifold. However, there are interesting
topological solutions with one unstable mode in 3D (sphalerons [1]) and 2D
(electroweak vortices [2, 3]) characterised by $\pi_{3}({\cal M})=Z$. Nambu
suggested a monopole-like configuration [4] which can be understood via a
local Hopf fibration $S^{3}\cong S^{2}\times S^{1}$ [5], which we discuss in
more detail in [6], such that $\pi_{2}(S^{2}\times S^{1})=Z$. However, this
configuration is known to be unstable and, if it were to be realised, it would
need to be combined with a string to form so-called “dumbbell” configurations
[7, 8]. The dynamics of such configurations has been linked with the
production of primordial magnetic fields [9, 10] and monopoles are being
actively searched for in laboratory experiments [11].

We will discuss the two-Higgs Doublet Model (2HDM) [12]; in particular one
where there are accidental symmetries. If the two doublets are $\Phi_{1}$ and
$\Phi_{2}$ with $\Phi^{T}=\begin{pmatrix}\Phi_{1}&\Phi_{2}\end{pmatrix}$ then
the potential is given by $V(\Phi_{1},\Phi_{2})=V_{2}+V_{4}$ where
$V_{2}=-\mu_{1}^{2}|\Phi_{1}|^{2}-\mu_{2}^{2}|\Phi_{2}|^{2}$ and
$V_{4}=\lambda_{1}|\Phi_{1}|^{4}+\lambda_{2}|\Phi_{2}|^{4}+\lambda_{3}|\Phi_{1}|^{2}|\Phi_{2}|^{2}+\lambda_{4}|\Phi_{2}^{{\dagger}}\Phi_{1}|^{2}$
which is $U(1)_{\rm PQ}$ symmetric. The symmetry can be enhanced to
$SO(3)_{\rm HF}$ if $\mu_{2}=\mu_{1}$, $\lambda_{2}=\lambda_{1}$ and
$\lambda_{4}=2\lambda_{1}-\lambda_{3}$[13, 14].

The particle spectrum is well understood, for example [15, 12]. There are 5
Higgs particles: two of which are CP even with masses $M_{h}$ and $M_{H}$, a
CP odd pseudo-scalar with $M_{A}$ and two charged Higgs particles with
$M_{H_{\pm}}$. In the Standard Model alignment limit the $h$ particle is that
detected by experiments at the LHC and a wide range of measurements suggest
that this limit should be close to being the case [16, 17, 18, 19, 20]. If
there is a global $U(1)_{\rm PQ}$ symmetry then $M_{A}=0$ and if this is
extended to $SO(3)_{\rm HF}$ then in addition one has $M_{H}=0$.

Topological defect solutions [21, 22] associated with these symmetries were
studied in detail in the context of the 2HDM in [14]. In particular there can
be domain wall [14, 23, 24, 25, 26], global vortex [14] and global monopole
solutions [27]. Motivating the present study is the observation, based on
field theory simulations from random initial conditions, that the vacuum is
not neutral in the core of these defects [23, 26, 27] contrary to the
assumption of [14]. We will see that this has profound implications. One
should note that there are also several papers discussing non-topological
configurations [28, 29, 30, 31] and other compound structures of topological
defects [32, 33] (without any neutral vacuum violation) in the 2HDM that can
be dynamically stable in certain regions of the parameter space. Here, we will
focus purely only on topologically stabilised objects.

Accidental symmetries can be gauged using a covariant kinetic term

$\displaystyle D_{\mu}\Phi$ $\displaystyle=$
$\displaystyle\bigg{[}(\sigma^{0}\otimes\sigma^{0})\partial_{\mu}+\textstyle\frac{1}{2}ig(\sigma^{0}\otimes\sigma^{a})W_{\mu}^{a}$
$\displaystyle+$
$\displaystyle\textstyle\frac{1}{2}ig^{\prime}(\sigma^{0}\otimes\sigma^{0})Y_{\mu}+\textstyle\frac{1}{2}ig^{\prime\prime}(\sigma^{a}\otimes\sigma^{0})V_{\mu}^{a}\bigg{]}\Phi\,,$

where $\sigma^{\mu}=(\sigma^{0},\sigma^{a})$ are the Pauli matrices including
the identity. $W_{\mu}^{a}$ and $Y_{\mu}$ are the Standard Model gauge fields,
with coupling constants $g$ and $g^{\prime}$ respectively, and $V_{\mu}^{a}$
are the new gauge fields associated with the accidental symmetries with
coupling constant $g^{\prime\prime}$ \- see [6] for more details on the
symmetry transformations. For the purposes of this work, we set $g^{\prime}=0$
in order to simplify the defect solutions, but note that we do not expect any
changes in the qualitative features for non-zero $g^{\prime}$.

Gauging the symmetries provides a natural mechanism for removing the Goldstone
modes associated with the accidental symmetries allowing for a potentially
viable model. In particular in the case of a $U(1)_{\rm PQ}$ symmetry the
Goldstone mode with $M_{A}=0$ becomes a massive gauge boson. There can be
interesting models constructed with these symmetries, for example, models that
can generate masses for neutrinos [34, 35, 36].

In this paper we will show that the Hopf fibration associated with the Nambu
monopole is realised on energetic grounds within the 2HDM when there is either
a $SO(3)_{\rm HF}$ or $U(1)_{\rm PQ}$ symmetry; something that we term
“Spontaneous Hopf Fibration” (SHF). We will find monopole and vortex solutions
for the case where the symmetries are gauged (although these can be easily
adapted to the global limit).

Parameterizations and topology: The 8 fields of the 2HDM can be
reparameterized as

$\Phi=\frac{v_{\rm
SM}}{\sqrt{2}}e^{\textstyle\frac{1}{2}i\chi}(\sigma^{0}\otimes
U_{L})\begin{pmatrix}0&f_{1}&f_{+}&f_{2}e^{i\xi}\end{pmatrix}^{T}\,,$ (2)

using 5 fields $f_{1,2,+}$, $\xi$ and $\chi$, and $U_{L}\in SU(2)_{\rm L}$
which has 3 degrees of freedom. The constant $v_{\rm SM}=246\,{\rm GeV}$ is
the Standard Model vacuum expectation value and $f_{+}=0$ corresponds to a
neutral vacuum with zero photon mass. These degrees of freedom can also be
encoded using bilinear forms,
$R^{\mu}=\Phi^{{\dagger}}(\sigma^{\mu}\otimes\sigma^{0})\Phi$,
$n^{a}=-\Phi^{\dagger}(\sigma^{0}\otimes\sigma^{a})\Phi$ and
${\tilde{R}}=2\Phi_{2}^{T}i\sigma^{2}\Phi_{1}$, which are useful for
understanding the topology of the vacuum manifold and the associated defects
[37, 14, 27]. In particular, the neutral vacuum violation discovered in the
core of defects in [23, 27] can be traced by $R_{+}=R_{\mu}R^{\mu}$, with
$R_{+}=0$ corresponding to a neutral vacuum. One finds that two of the $U_{\rm
L}$ degrees of freedom are encoded in $\hat{n}^{a}$, with the other associated
with rotations about this axis, and the hypercharge degree of freedom, $U_{\rm
Y}$, is encoded in ${\tilde{R}}\propto\exp[i\chi]$. In contrast, $R^{\mu}$ is
invariant under the Standard Model symmetries and contains the degrees of
freedom that will, in general, change the potential. The $R^{0}$ component is
$|\Phi|^{2}$ and so does not contribute to the topology of the vacuum
manifold. The topological non-triviality of the vacuum manifold can,
therefore, be most easily extracted by looking at the remaining three-vector
$R^{a}$, which will contain all of the degrees of freedom associated with any
additional symmetry transformations, $U_{H}$.

A notable difference between this topology and that of the ’t Hooft-Polyakov
monopole [38, 39] and Nielsen-Olesen vortex [40] is that the topology lives in
a space associated with these bilinear forms, rather than the fields
themselves, which means that a half twist in field space can be topologically
non-trivial and, in general, there is a factor of $2$ difference between the
topological degree of a field configuration and what one might naively expect.
A consequence of this is that simple field configurations with unit winding
often have discontinuities which must be resolved by the attachment of another
soliton. It is for this reason that the Nambu monopole [4] has a string
“emanating” from one of its poles.

In the $\text{SO(3)}_{\text{HF}}$ case $R^{a}$ only contributes to the
potential with a term $\propto R^{a}R^{a}$, so rotations between the three
components of $R^{a}$ are a symmetry of the potential, generating an $S^{2}$
component of the vacuum manifold. Note that the topology is not $S^{3}$ as one
might expect because $U_{H}$ cannot perform rotations about $R^{a}$ — this
degree of freedom is already contained within $U_{L}$ for rotations about
$n^{a}$. Similarly, in the $\text{U(1)}_{\text{PQ}}$ case, the $R^{3}$
component splits off from the other two but there remains a symmetry for
rotations between $R^{1}$ and $R^{2}$, which is responsible for an $S^{1}$
direction. In general, there is an additional $S^{1}\times S^{3}$ associated
with the hypercharge and isospin symmetries because the degeneracy between
$U_{Y}$ and one of the directions in $U_{L}$ is broken by $f_{+}$. Since this
results in a massive photon, we avoid this scenario and choose the parameters
so that $f_{+}=0$ in the vacuum, restoring the degeneracy so that the SM
symmetries only contribute $S^{3}$. Therefore, the topology of the vacuum
manifolds associated with the accidental symmetries are ${\cal M}=S^{2}\times
S^{3}$ for the case of $SO(3)_{\rm HF}$ and ${\cal M}=S^{1}\times S^{3}$ for
$U(1)_{\rm PQ}$ [14], which admit monopoles and vortices due to the non-
trivial homotopy groups $\pi_{2}(S^{2}\times S^{3})=\mathbb{Z}$ and
$\pi_{1}(S^{1}\times S^{3})=\mathbb{Z}$, respectively.

Monopole solutions: In the $\text{SO(3)}_{\text{HF}}$ symmetric model, using
3D spherical polar coordinates $(r,\theta,\phi)$, we can construct the
monopole ansatz for the scalar field [27, 6]

$\Phi(r,\theta,\phi)=\frac{v_{\text{SM}}}{2\sqrt{2}}\begin{pmatrix}-(k+k_{+})\sin\theta
e^{-i\phi}\\\ (k-k_{+})+(k+k_{+})\cos\theta\\\
-(k-k_{+})+(k+k_{+})\cos\theta\\\ (k+k_{+})\sin\theta
e^{i\phi}\par\end{pmatrix}\,,$ (3)

where $k=k(r)$ and $k_{+}=k_{+}(r)$ are functions constructed from $f_{1,2,+}$
that retain the ability to change the potential energy while the other degree
of freedom (as well as $\xi$) is used to wind around the vacuum manifold. This
ansatz has the property that $R^{a}=n^{a}=(k^{2}-k_{+}^{2})\hat{r}^{a}$. The
feature $\hat{R}^{a}=\hat{r}^{a}$ is necessary for a monopole configuration
with unit winding (or related to this by a homomorphism) and this structure,
by itself, would give rise to a 2HDM equivalent of the Nambu monopole — with a
divergence in the gradient energy that necessitates the emergence of a string
from one of the poles. However, the isospin degrees of freedom, contained
within $U_{L}$, can resolve this divergence when we also have that
$\hat{n}^{a}=\hat{r}^{a}$, which gives the appearance of an underlying
topological complexity in the structure of $n^{a}$, but it occurs for
energetic reasons that are only indirectly topological. The degree of freedom
associated with rotations about $n^{a}$ contains no structure — only the
possibility for global transformations — and it is this effect that we have
termed SHF. The $S^{3}$ of the SM becomes $S^{2}\times S^{1}$ with the $S^{2}$
part inheriting the same twists as the $S^{2}$ of the
$\text{SO(3)}_{\text{HF}}$ symmetry and the $S^{1}$ part containing nothing
except possible global rotations.

For the gauge fields, we choose to work in the temporal gauge such that the
time components of the gauge fields are zero and make the ansatz
$gW_{i}^{a}=-\frac{1}{r}h(r){\epsilon^{a}}_{ij}\hat{r}^{j}$ and
$g^{\prime\prime}V_{i}^{a}=-\frac{1}{r}H(r){\epsilon^{a}}_{ij}\hat{r}^{j}$.
After the standard rescaling to reduce the number of significant parameters in
the model by two, we find that the energy functional under this ansatz is

$\displaystyle E$ $\displaystyle=\frac{4\pi v_{\text{SM}}}{g}\int
r^{2}dr\bigg{\\{}\frac{1}{2}\bigg{(}\frac{dk}{dr}\bigg{)}^{2}+\frac{1}{2}\bigg{(}\frac{dk_{+}}{dr}\bigg{)}^{2}$
$\displaystyle+\frac{1}{8r^{2}}(k+k_{+})^{2}(h+H-2)^{2}+\frac{1}{8r^{2}}(k-k_{+})^{2}(h-H)^{2}$
$\displaystyle+\frac{1}{2r^{2}}\bigg{[}2\bigg{(}\frac{dh}{dr}\bigg{)}^{2}+\frac{1}{r^{2}}h^{2}(2-h)^{2}\bigg{]}+\frac{1}{2\tilde{g}^{2}r^{2}}\bigg{[}2\bigg{(}\frac{dH}{dr}\bigg{)}^{2}$
$\displaystyle+\frac{1}{r^{2}}H^{2}(2-H)^{2}\bigg{]}+\frac{\tilde{\lambda}}{4}\bigg{[}(k^{2}+k_{+}^{2}-1)^{2}-\zeta_{4}k^{2}k_{+}^{2}\bigg{]}\bigg{\\}}\,,$
(4)

where the remaining parameters are $\tilde{g}=g^{\prime\prime}/g$,
$\tilde{\lambda}=\lambda/g^{2}$ and $\zeta_{4}=\lambda_{4}/\lambda$. Neutral
vacuum violation occurs when both $k$ and $k_{+}$ are simultaneously non-zero
and we choose, by convention, $k_{+}$ to be zero in the vacuum. We can perform
a simple analysis (neglecting gradient energy contributions) to predict when
there will be neutral vacuum violation in the core of the monopole by looking
at the effective mass $m_{k_{+}}^{2}=\lambda[2(k^{2}-1)-\zeta_{4}k^{2}]/4$. If
a monopole were to have $k_{+}=0$ everywhere, then the effective mass at the
core of the monopole (where $k=0$) would be $-\frac{\tilde{\lambda}}{2}$,
which is always negative and independent of $\zeta_{4}$. The presence of this
negative mass term indicates that the energy would be reduced if $k_{+}\neq 0$
and therefore we expect neutral vacuum violation to be a generic effect in the
core of 2HDM monopoles.

In figure 1 we present the energy and $R_{+}(0)$ (fixing $v_{\text{SM}}=1)$ as
a function of the mass ratio
$\epsilon=M_{H\pm}/M_{h}=\frac{1}{2}\sqrt{-\zeta_{4}}$ and we also show an
example solution in the inset plot. As expected, for all values of $\epsilon$
presented here, there is neutral vacuum violation in the core of the monopole,
although it decreases as $\epsilon$ grows and appears to be approaching zero,
while conversely, the energy grows with $\epsilon$ and appears to be
approaching a maximum. Perhaps the most noticeable feature of the solution is
that $k+k_{+}=0$ at the centre — in fact this is enforced by the gradient
energy and is true for all values of $\epsilon$. The effects of
$\tilde{\lambda}$ can be broadly described as changing the length scale ratio
between the scalar fields and the vector fields and, similarly, $\tilde{g}$
changes the length scale ratio between the two gauge fields.

Figure 1: The variation of $E$ (black) and $R_{+}$ at the monopole core
(dotted orange) as a function of the mass ratio $\epsilon$, with the other
parameters fixed to $\tilde{\lambda}=1$ and $\tilde{g}^{2}=2.$ The inset plot
presents the field profiles of an example solution ($\epsilon=1$) and also
displays the energy density, $\varepsilon$, and $R_{+}$.

Vortex solutions: In the $\text{U(1)}_{\text{PQ}}$ symmetric case, using plane
polar coordinates $(r,\theta)$, we can make the vortex ansatz
$\Phi=\Phi(r,\theta)$ from equation (2) with $f_{i}=f_{i}(r)$ for $i=1,2,+$
and $\xi=\theta$ [6] which has a similar property to the monopole in that
$\hat{R}^{b}=\hat{n}^{b}=\hat{r}^{b}$, where $b\in[1,2]$, although now, unlike
the monopole case, $n^{b}=0$ in the vacuum. The remaining components are
$R^{3}=f_{1}^{2}-f_{+}^{2}-f_{2}^{2}$ and
$n^{3}=f_{1}^{2}-f_{+}^{2}+f_{2}^{2}$. We note that, in a similar way to the
Nambu monopole solution, there is a string-like configuration, characterised
by $\pi_{1}(S^{2}\times S^{1})=Z$ and with structure only in $R^{b}$, where
the divergence in the gradient energy is resolved by attaching a domain wall
to one side. This is related to an unstable (due to the tension in the wall)
configuration in the SM where $\hat{n}^{b}=\hat{r}^{b}$. Once again, for the
2HDM configuration, we can use the isospin rotations to resolve the divergence
without the domain wall. The SHF acts here, again, to split $S^{3}\to
S^{2}\times S^{1}$, but now it is the $S^{1}$ part that inherits the twists of
the $\text{U(1)}_{\text{PQ}}$.

Now we choose to make a gauge transformation that absorbs the phase winding of
the scalar field into $W_{i}^{a}$ and $V_{i}^{3}$ (the only new gauge field in
the $\text{U(1)}_{\text{PQ}}$ model) so that we can make the ansatz
$gW_{i}^{a}=\frac{1}{r}[h_{1}(r)\hat{x}^{a}+(1-h_{3}(r))\hat{z}^{a}]\hat{\theta}_{i}+h_{2}(r)\hat{y}^{a}\hat{r}_{i}$
and $g^{\prime\prime}V_{i}^{3}=\frac{1}{r}(1-H(r))\hat{\theta}_{i}$. Again, if
we make the standard rescalings then we can express the energy per unit length
of the string as

$\displaystyle E$ $\displaystyle=2\pi v_{\text{SM}}^{2}\eta_{1}^{2}\int
rdr\bigg{\\{}\frac{1}{2}\bigg{(}\frac{df_{1}}{dr}\bigg{)}^{2}+\frac{1}{2}\bigg{(}\frac{df_{+}}{dr}\bigg{)}^{2}+\frac{1}{2}\bigg{(}\frac{df_{2}}{dr}\bigg{)}^{2}$
$\displaystyle-\frac{1}{2}h_{2}\bigg{(}f_{+}\frac{df_{2}}{dr}-f_{2}\frac{df_{+}}{dr}\bigg{)}+\frac{1}{8}\bigg{(}\frac{h_{1}^{2}}{r^{2}}+h_{2}^{2}$
$\displaystyle+\frac{1}{r^{2}}(h_{3}-H)^{2}\bigg{)}(f_{1}^{2}+f_{+}^{2})-\frac{1}{2r^{2}}h_{1}(1-H)f_{+}f_{2}$
$\displaystyle+\frac{1}{8}\bigg{(}\frac{h_{1}^{2}}{r^{2}}+h_{2}^{2}+\frac{1}{r^{2}}(2-h_{3}-H)^{2}\bigg{)}f_{2}^{2}+\frac{1}{2r^{2}}\bigg{[}\bigg{(}\frac{dh_{1}}{dr}\bigg{)}^{2}$
$\displaystyle+\bigg{(}\frac{dh_{3}}{dr}\bigg{)}^{2}-2h_{2}\bigg{(}h_{1}\frac{dh_{3}}{dr}+(1-h_{3})\frac{dh_{1}}{dr}\bigg{)}+h_{1}^{2}h_{2}^{2}$
$\displaystyle+(1-h_{3})^{2}h_{2}^{2}\bigg{]}+\frac{1}{2\tilde{g}^{2}r^{2}}\bigg{(}\frac{dH}{dr}\bigg{)}^{2}+\frac{\tilde{\lambda}_{1}}{4}\bigg{[}(f_{1}^{2}-1)^{2}$
$\displaystyle+\zeta_{2}(f_{+}^{2}+f_{2}^{2}-\tilde{\eta}^{2})^{2}+\zeta_{3}f_{1}^{2}(f_{+}^{2}+f_{2}^{2})+\zeta_{4}f_{1}^{2}f_{2}^{2}\bigg{]}\bigg{\\}},$
(5)

where $\tilde{g}=g^{\prime\prime}/g$, $\tilde{\eta}=\eta_{2}/\eta_{1}$,
$\tilde{\lambda}_{1}=\lambda_{1}/g^{2}$ and
$\zeta_{i}=\lambda_{i}/\lambda_{1}$ so that the model is left with 6
parameters ($\eta_{i}$ is defined by the relationship
$\mu_{i}^{2}=\lambda_{i}\eta_{i}^{2}v_{\text{SM}}^{2}$). We can perform an
effective mass analysis here too, just as in the monopole case. The relevant
effective mass is
$m_{f_{+}}^{2}={\tilde{\lambda}}_{1}[2\zeta_{2}(f_{2}^{2}-\tilde{\eta}^{2})+\zeta_{3}f_{1}^{2}]/4$
and if the string were to have $f_{+}=0$ everywhere, as well as $f_{1}(0)=1$
and $f_{2}(0)=0$, then this would take the value
$\tilde{\lambda}_{1}[-2\zeta_{2}\tilde{\eta}^{2}+\zeta_{3}]/4$, so the neutral
vacuum violation in the core is not generic but depends upon the sign of
$\zeta_{3}-2\zeta_{2}\tilde{\eta}^{2}$. Note that, $f_{2}(0)=0$ is enforced by
the winding of the string but $f_{1}(0)=1$ is an overly simplistic assumption,
even if $f_{+}=0$ everywhere, so this approach will be less accurate for
vortices than for the monopoles.

In figure 2 we present the field profiles of a string solution with a separate
inset plot showing the energy density and $R_{+}$ for the same solution and in
figure 3 we show how $R_{+}$, at the core of the string, varies with the mass
ratios $\epsilon=M_{H\pm}/M_{h}$ and $\delta=M_{H}/M_{h}$, in the alignment
limit and with $\tan\beta=1$ (which sets the vacuum values of the fields so
that $|\Phi_{1}|^{2}=|\Phi_{2}|^{2}=v_{\Phi}^{2}$). Note that we have rescaled
$R_{+}\to R_{+}/(v_{\Phi}v_{\text{SM}})^{4}$ implicitly in these plots.

From the profiles we see that the solution has retained a feature that is
similar to one observed for the monopole solutions, namely that $f_{1}=f_{+}$
at the centre. However, this is no longer guaranteed by the gradient energy
and, in fact, is only an approximate equality that occurs in a subset of the
parameter space. The parameters $\tilde{\lambda}_{1}$ and $\tilde{g}^{2}$ play
a very similar role to the equivalent parameters from the monopole case but
the other parameters have more complicated effects on the solutions that are
difficult to broadly summarise in this way. From the contour plot we can see
that there is a clear transition across the line which is approximately
$\epsilon=\delta$, corresponding to $\zeta_{3}=2$. Due to our fixed parameter
choices we have $\zeta_{2}=\tilde{\eta}^{2}=1$ and therefore this is
consistent with the behaviour that we predicted from the effective mass
analysis.

Discussion and conclusions: The solutions that we have presented in this paper
are evidence of a new mechanism — that we have called “Spontaneous Hopf
Fibration” — at work in the 2HDM. It allows for a topologically non-trivial
subspace of the vacuum manifold to imprint itself onto another, topologically
trivial section of the manifold. This results in solitons that appear to have
topological structure in the SM degrees of freedom but, in fact, it is purely
caused by energetics. The coupling between the $S^{3}$ of the SM to the rest
of the vacuum manifold, in the gradient energy term, causes a Hopf fibration
of the space $S^{3}\to S^{2}\times S^{1}$, with the appropriate component of
this space taking on structure to match the winding around $S^{2}$ for
monopoles and $S^{1}$ for strings.

In [27] we present evidence from simulations of the global 2HDM model in which
stable monopoles form that have neutral vacuum violation in the cores and a
structure in the bilinear vectors that is the same as what one would expect
from the solutions presented here. In [6] we do the same for the case of
strings. These simulations suggest that the solutions we have found are those
most relevant to the study of topological defects in the 2HDM.

A phenomenologically relevant consequence of the SHF in the 2HDM is that it
allows for the neutral vacuum condition to be violated inside the core of the
defects, generating a non-zero mass for the photon. In the case of strings,
this is dependent upon the choice of parameters, however in the monopole case
it is predicted to always occur if the gradient energy contributions are
neglected. The interaction between photons and superconducting defects has
been analysed in [41] for a toy model but this work has opened up the
possibility for novel interactions between standard model particles and 2HDM
defects which is deserving of more investigation. In [27] it was observed that
the additional structure of global 2HDM monopoles did not affect the scaling
of their number density, but other potential cosmological consequences warrant
further studies of these defects.

Figure 2: The field profiles for a string solution with
$\tilde{\lambda}_{1}=1$, $\tilde{g}^{2}=1$, $\tilde{\eta}^{2}=1$,
$\zeta_{2}=1$, $\zeta_{3}=-0.4$ and $\zeta_{4}=-0.8$ — the last four parameter
choices corresponding to the alignment limit with $\epsilon=1$, $\delta=2$,
$\tan\beta=1$. The inset plot displays the energy density, $\varepsilon$, and
$R_{+}$ for the same solution. Figure 3: A contour plot showing how $R_{+}$ at
the string core varies with the mass ratios $\epsilon$ and $\delta$ in the
alignment limit and with $\tilde{\lambda_{1}}=\tilde{g}^{2}=\tan\beta=1$. Note
that the dark blue colour on most of the lower right region of the plot is off
the bottom of the colour scale - corresponding to a neutral vacuum.

We would like to conclude by emphasising that, although we have discussed SHF
in the context of the 2HDM, we suggest that it could be a more general effect
that can occur in other models that have vacuum manifolds constructed from
coupled subspaces, when at least one of them is topologically non-trivial.

Acknowledgements: We would like to thank Jeff Forshaw for helpful comments on
the phenomenological viability of these models. RB would like to thank
Apostolos Pilaftsis, Gary Brawn and Dominic Viatic for their collaboration on
the 2HDM which ultimately lead to this work.

## References

  * Klinkhamer and Manton [1984] F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984).
  * Vachaspati [1993] T. Vachaspati, Nucl. Phys. B 397, 648 (1993).
  * Achucarro and Vachaspati [2000] A. Achucarro and T. Vachaspati, Phys. Rept. 327, 347 (2000), arXiv:hep-ph/9904229 .
  * Nambu [1977] Y. Nambu, Nuclear Physics B 130, 505 (1977).
  * Gibbons  _et al._ [1992] G. W. Gibbons, M. E. Ortiz, F. Ruiz Ruiz, and T. M. Samols, Nucl. Phys. B 385, 127 (1992), arXiv:hep-th/9203023 .
  * [6] See Supplemental Material for more detail on the Hopf fibration, symmetry transformations, ansatz construction and simulations of strings, which includes Ref. [27].
  * Lazarides  _et al._ [2021] G. Lazarides, Q. Shafi, and T. Vachaspati, Phys. Rev. D 104, 035020 (2021), arXiv:2106.07800 [hep-ph] .
  * Patel and Vachaspati [2023] T. Patel and T. Vachaspati, Phys. Rev. D 107, 093010 (2023).
  * Vachaspati [2021] T. Vachaspati, Rept. Prog. Phys. 84, 074901 (2021), arXiv:2010.10525 [astro-ph.CO] .
  * Patel and Vachaspati [2022] T. Patel and T. Vachaspati, JHEP 2022, 59 (2022).
  * Acharya  _et al._ [2022] B. Acharya _et al._ (MoEDAL), Nature 602, 63 (2022), arXiv:2106.11933 [hep-ex] .
  * Branco  _et al._ [2012] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, Phys. Rept. 516, 1 (2012), arXiv:1106.0034 [hep-ph] .
  * Ferreira  _et al._ [2009] P. M. Ferreira, H. E. Haber, and J. a. P. Silva, Phys. Rev. D 79, 116004 (2009).
  * Battye  _et al._ [2011] R. A. Battye, G. D. Brawn, and A. Pilaftsis, JHEP 08, 020 (2011), arXiv:1106.3482 [hep-ph] .
  * Pilaftsis and Wagner [1999] A. Pilaftsis and C. E. M. Wagner, Nucl. Phys. B553, 3 (1999), arXiv:hep-ph/9902371 [hep-ph] .
  * Celis  _et al._ [2013] A. Celis, V. Ilisie, and A. Pich, JHEP 12, 095 (2013), arXiv:1310.7941 [hep-ph] .
  * Chiang and Yagyu [2013] C.-W. Chiang and K. Yagyu, JHEP 07, 160 (2013), arXiv:1303.0168 [hep-ph] .
  * Chen  _et al._ [2013] C.-Y. Chen, S. Dawson, and M. Sher, Phys. Rev. D 88, 015018 (2013), [Erratum: Phys.Rev.D 88, 039901(E) (2013)], arXiv:1305.1624 [hep-ph] .
  * Chowdhury and Eberhardt [2018] D. Chowdhury and O. Eberhardt, JHEP 05, 161 (2018), arXiv:1711.02095 [hep-ph] .
  * Haller  _et al._ [2018] J. Haller, A. Hoecker, R. Kogler, K. Mönig, T. Peiffer, and J. Stelzer, Eur. Phys. J. C78, 675 (2018), arXiv:1803.01853 [hep-ph] .
  * Vilenkin and Shellard [1994] A. Vilenkin and E. Shellard,  _Cosmic strings and other topological defects_ , Cambridge monographs on mathematical physics (Cambridge University Press, Cambridge, 1994).
  * Manton and Sutcliffe [2004] N. Manton and P. Sutcliffe,  _Topological Solitons_ (Cambridge University Press, 2004).
  * Battye  _et al._ [2021] R. A. Battye, A. Pilaftsis, and D. G. Viatic, JHEP 01, 105 (2021), arXiv:2006.13273 [hep-ph] .
  * Battye  _et al._ [2020] R. A. Battye, A. Pilaftsis, and D. G. Viatic, Phys. Rev. D 102, 123536 (2020).
  * Chen  _et al._ [2020] N. Chen, T. Li, Z. Teng, and Y. Wu, JHEP 10, 081 (2020), arXiv:2006.06913 [hep-ph] .
  * Law and Pilaftsis [2021] K. H. Law and A. Pilaftsis, (2021), arXiv:2110.12550 [hep-ph] .
  * Battye  _et al._ [2023] R. A. Battye, S. J. Cotterill, and D. G. Viatic, Physics Letters B 844, 138091 (2023).
  * Bachas and Tomaras [1996] C. Bachas and T. N. Tomaras, Phys. Rev. Lett. 76, 356 (1996).
  * Bachas  _et al._ [1999] C. Bachas, B. Rai, and T. N. Tomaras, Phys. Rev. Lett. 82, 2443 (1999).
  * Bachas  _et al._ [1996] C. Bachas, P. Tinyakov, and T. Tomaras, Physics Letters B 385, 237 (1996).
  * Earnshaw and James [1993] M. A. Earnshaw and M. James, Phys. Rev. D 48, 5818 (1993).
  * Eto  _et al._ [2020] M. Eto, Y. Hamada, M. Kurachi, and M. Nitta, JHEP 07, 004 (2020), arXiv:2003.08772 [hep-ph] .
  * Eto  _et al._ [2021] M. Eto, Y. Hamada, and M. Nitta, (2021), arXiv:2111.13345 [hep-ph] .
  * Ko  _et al._ [2012] P. Ko, Y. Omura, and C. Yu, Physics Letters B 717, 202 (2012).
  * Camargo  _et al._ [2019] D. Camargo, A. Dias, T. de Melo, and F. Queiroz, JHEP 2019 (2019), 10.1007/JHEP04(2019)129.
  * Grimus  _et al._ [2009] W. Grimus, L. Lavoura, and B. Radovčić, Physics Letters B 674, 117 (2009).
  * Ivanov [2008] I. P. Ivanov, Phys. Rev. D77, 015017 (2008), arXiv:0710.3490 [hep-ph] .
  * ’t Hooft [1974] G. ’t Hooft, Nucl. Phys. B 79, 276 (1974).
  * Polyakov [1974] A. M. Polyakov, JETP Lett. 20, 194 (1974).
  * Nielsen and Olesen [1973] H. B. Nielsen and P. Olesen, Nucl. Phys. B 61, 45 (1973).
  * Battye and Viatic [2021] R. A. Battye and D. G. Viatic, Phys. Lett. B 823, 136730 (2021), arXiv:2110.13668 [hep-ph] .


# LAST at CMCL 2021 Shared Task: Predicting Gaze Data During Reading with a
Gradient Boosting Decision Tree Approach

Yves Bestgen
Laboratoire d’analyse statistique des textes - LAST
Institut de recherche en sciences psychologiques
Université catholique de Louvain
Place Cardinal Mercier, 10 1348 Louvain-la-Neuve, Belgium
<EMAIL_ADDRESS>

###### Abstract

A LightGBM model fed with target word lexical characteristics and features
obtained from word frequency lists, psychometric data and bigram association
measures has been optimized for the 2021 CMCL Shared Task on Eye-Tracking Data
Prediction. It obtained the best performance of all teams on two of the five
eye-tracking measures to predict, allowing it to rank first on the official
challenge criterion and to outperform all deep-learning based systems
participating in the challenge.

##  1 Introduction

This paper describes the system proposed by the Laboratoire d’analyse
statistique des textes (LAST) for the Cognitive Modeling and Computational
Linguistics (CMCL) Shared Task on Eye-Tracking Data Prediction. This task is
receiving more and more attention due to its importance in modeling human
language understanding and improving NLP technology (Hollenstein et al., 2019;
Mishra and Bhattacharyya, 2018).

As one of the objectives of the organizers is to “compare the capabilities of
machine learning approaches to model and analyze human patterns of reading”
(https://cmclorg.github.io/shared_task), I have chosen to adopt a generic
point of view with the main objective of determining what level of performance
can achieve a system derived from the one I developed to predict the lexical
complexity of words and polylexical expressions (Shardlow et al., 2021). That
system was made up of a gradient boosting decision tree prediction model fed
with features obtained from word frequency lists, psychometric data, lexical
norms and bigram association measures. If there is no doubt that predicting
lexical complexity is a different problem, one can think that the features
useful for it also play a role in predicting eye movement during reading.

The next section summarizes the main characteristics of the challenge. Then
the developed system is described in detail. Finally, the results in the
challenge are reported along with an analysis performed to get a better idea
of the factors that affect the system performance.

##  2 Data and Task

The eye-tracking data for this shared task were extracted from the Zurich
Cognitive Language Processing Corpus (ZuCo 1.0 and ZuCo 2.0, Hollenstein et
al., 2018, 2020). It contains gaze data for 991 sentences read by 18
participants during a normal reading session. The learning set consisted in
800 sentences and the test set in 191 sentences.

The task was to predict five eye-tracking features, averaged across all
participants and scaled in the range between 0 and 100, for each word of a
series of sentences: (1) the total number of fixations (nFix), (2) the
duration of the first fixation (FFD), (3) the sum of all fixation durations,
including regressions (TRT), (4) the sum of the duration of all fixations
prior to progressing to the right, including regressions to previous words
(GPT), and (5) the proportion of participants that fixated the word (fixProp).
These dependent variables (DVs) are described in detail in Hollenstein et al.
(2021). The submissions were evaluated using the mean absolute error (MAE)
metric and the systems were ranked according to the average MAE across all
five DVs, the lowest being the best.

As the DVs are of different natures (number, proportion and duration), their
mean and variance are very different. The mean of fixProp is 21 times greater
than that of FFD and its variance 335 times. Furthermore, while nFix and
fixProp were scaled independently, FFD, GPT and TRT were scaled together. For
that reason, the mean and dispersion of these three measures are quite
different: FFD must necessarily be less than or equal to TRT and GPT111The
relation between TRT and GPT is not obvious to me since one can be larger or
smaller than the other in a significant number of cases.. These two factors
strongly affect the importance of the different DVs in the final ranking.

##  3 System

###  3.1 Procedure to Build the Models

The regression models were built by the 2.2.1 version of the LightGBM software
(Ke et al., 2017), a well-known implementation of the gradient boosting
decision tree approach. This type of model has the advantage of not requiring
feature preprocessing, such as a logarithmic transformation, since it is
insensitive to monotonic transformations, and of including many parameters
allowing a very efficient overfit control. It also has the advantage of being
able to directly optimize the MAE.

Sentence preprocessing and feature extraction as well as the post-processing
of the LightGBM predictions were performed using custom SAS programs running
in SAS University (still freely available for research at
https://www.sas.com/en_us/software/university-edition.html). Sentences were
first lemmatized by the TreeTagger (Schmid, 1994) to get the lemma and POS-tag
of each word. Special care was necessary to match the TreeTagger tokenization
with the Zuco original one. Punctuation marks and other similar symbols (e.g.,
"(" or "$") were simply disregarded as they were always bound to a word in the
tokens to predict. The attribution to the words of the values on the different
lists was carried out in two stages: on the basis of the spelling form when it
is found in the list or of the lemma if this is not the case.

The features used in the final models as well as the LightGBM parameters were
optimized by a 5-fold cross validation procedure, using the sentence and not
the token as the sampling unit. The number of boosting iterations was set by
using the LightGBM early stopping procedure which stops training when the MAE
on the validation fold does not improve in the last 200 rounds. The predicted
values which were outside the [0, 100] interval were brought back in this one,
which makes it possible to improve the MAE very slightly.

###  3.2 Features

To predict the five DVs, five different models were trained. The only
differences between them were in the LightGBM parameters. There were thus all
based on exactly the same features, which are described below.

#### Target Word Length.

The length in characters of the preceding word, the target word and the
following one.

#### Target Word Position.

The position of the word in the sentence encoded in two ways: the rank of the
word going from 1 to the sentence total number of words and the ratio between
the rank of the word and the total number of words.

#### Target Word POS-tag and Lemma.

The POS-tag and lemma for the target word and the preceding one.

#### Corpus Frequency Features.

Frequencies in corpora of words were either calculated from a corpus or
extracted from lists provided by other researchers. The following seven
features have been used:

  * •

The (unlemmatized) word frequencies in the British National Corpus (BNC,
http://www.natcorp.ox.ac.uk/).

  * •

The Facebook frequency norms for American English and British English in
Herdagdelen and Marelli (2017).

  * •

The Rovereto Twitter Corpus frequency norms (Herdagdelen and Marelli, 2017).

  * •

The USENET Orthographic Frequencies from Shaoul and Chris (2006).

  * •

The Hyperspace Analogue to Language (HAL) frequency norms provided by (Balota
et al., 2007) for more that 40,000 words.

  * •

The frequency word list derived from Google’s ngram corpora available at
https://github.com/hackerb9/gwordlist.

#### Features from Lexical Norms.

The lexical norms of Age of Acquisition and Familiarity were taken from the
Glasgow Norms which contain judges’ assessment of 5,553 English words (Scott
et al., 2019).

#### Lexical Characteristics and Behavioral Measures from ELP.

Twenty-three indices were extracted from the English Lexicon Project (ELP,
Balota et al., 2007; Yarkoni et al., 2008), a database that contains, for more
than 40,000 words, reaction time and accuracy during lexical decision and
naming tasks, made by many participants, as well as lexical characteristics
(https://elexicon.wustl.edu/). Eight indices come from the behavioral
measures, four for each task: average response latencies (raw and
standardized), standard deviations, and accuracies. Fourteen indices come from
the “Orthographic, Phonological, Phonographic, and Levenshtein Neighborhood
Metrics” section of the dataset. These are all the metrics provided except
Freq_Greater, Freq_G_Mean, Freq_Less, Freq_L_Mean, and Freq_Rel. These are
variables whose initial analyzes showed that they were redundant with those
selected. The last feature is the average bigram count of a word.

#### Bigram Association Measures.

These features indicate the degree of association between the target word and
the one that precedes it according to a series of indices calculated on the
basis of the frequency in a reference corpus (i.e., the BNC) of the bigram and
that of the two words that compose it, using the following association
measures (AMs): pointwise mutual information and t-score (Church and Hanks,
1990), z-score (Berry-Rogghe, 1973), log-likelihood Chi-square test (Dunning,
1993), simple-ll (Evert, 2009), Dice coefficient (Kilgarriff et al., 2014) and
the two delta-p (Kyle et al., 2018). Most of the formulas to compute these AMs
are also provided in Evert (2009) and in Pecina (2010). As these features mix
together the assets of both collocations (by using association scores) and
ngrams (by using contiguous pairs of words), Bestgen and Granger (2014) refer
to them as collgrams. They make it possible not to rely exclusively on the
frequency of the bigram in the corpus, which can be misleading because a
bigram may be observed frequently, not because of its phraseological nature,
but because it is made up of very frequent words (Bestgen, 2018). Conversely,
a relatively rare bigram, composed of rare words, may be typical of the
language. Since word frequency is already accounted for by the corpus
frequency features, it was desirable to employ indices that reduce the impact
of this factor. Originating in works in lexicography and foreign language
learning (Church and Hanks, 1990; Durrant and Schmitt, 2009; Bestgen, 2017,
2019), they have recently shown their usefulness in predicting the lexical
complexity of multi-word expressions (Bestgen, 2021). In the present case, it
is assumed that these indices can serve as a proxy of the next word
predictability (Kliegl et al., 2004).

#### Feature coverage.

Some words to predict are not present in these lists and the corresponding
score is thus missing. Based on the complete dataset provided by the
organizers, it happens in:

  * •

1% (Google ngram) to 17% (Facebook and Twitter) of the tokens for the corpus
frequency features,

  * •

9% for the ELP Lexical Characteristics, but a few features have as much as 41%
missing values,

  * •

11% for the ELP Behavioral Measures,

  * •

18% for the Bigram AMs.

In total, sixteen tokens have missing values for all these features (Corpus
Frequency, Lexical Characteristics and Behavioral Measures from ELP, and
Bigram Association Measures). These tokens have however received values for
the length and position features. All the missing values were handled by
LightGBM default procedure.

Parameters | Run 1 | Run 2
---|---|---
bagging_fraction | 0.66 | 0.70
bagging_freq | 5 | 5
feature_fraction | 0.09 | 0.85
learning_rate | 0.0095 | 0.0050
max_depth | 11 | no limit
max_bin | 64 | 64
min_data_in_bin | 2 | 5
max_leaves | 11 | 30
min_data_in_leaf | 7 | 5
n_iter | 4800 | (see text)
Table 1: LightGBM parameters for the first two runs. Team | Run | Mean | nFix | FFD | GPT | TRT | fixProp
---|---|---|---|---|---|---|---
LAST | 3 | 3.8134 | 3.879 | 0.655 | 2.197 | 1.524 | 10.812
LAST | 2 | 3.8159 | 3.886 | 0.655 | 2.199 | 1.523 | 10.817
TALEP | 1 | 3.8328 | 3.761 | 0.662 | 2.180 | 1.486 | 11.076
LAST | 1 | 3.8664 | 3.943 | 0.662 | 2.237 | 1.545 | 10.944
TorontoCL | 2 | 3.9287 | 3.944 | 0.671 | 2.227 | 1.516 | 11.286
Table 2: Performance (MAE) for the five best runs submitted to the challenge.
Best scores are bolded.

##  4 Analyses and Results

###  4.1 Models Submitted to the Challenge

During the test phase, teams were allowed to submit three runs. My three
submissions were all based on the features described above, the only
differences between them resulting from changes in the LigthGBM parameters.
They were set at their default values except those shown in Table 1. The
official performances of the top five challenge submissions are given in Table
2.

The first submission was based on the parameters selected during the
development phase. They were identical for the five DVs. For the other two
submissions, a random grid search coded in python was used to try optimizing
the parameters independently for each DV. The parameter space for this first
random search is provided in Appendix A. As the measure of the challenge is
the MAE averaged across the five DVs and as the system MAE for fixProp was up
to 15 times higher than that of the other DVs, the optimized parameters for
this variable were selected. Additional analyzes showed that they also made it
possible to improve performance on the four other DVs. Their values are given
in Table 1. Certain initial choices were only slightly modified. The value of
other parameters such as the maximum number of leaves and the feature fraction
were markedly increased, suggesting that the risk of overfit was relatively
low (see https://lightgbm.readthedocs.io/en/latest/Parameters-Tuning.html). In
this system, the number of iterations was optimized (thanks to the early
stopping procedure) for each DV and sets at the fourth highest value: 3,740
for nFix, 3,829 for TRT, 2,861 for GPT, 3,497 for FFD, and 3,305 for fixProp.

Models | MAE | %MAE | %nFix | %FFD | %GPT | %TRT | %fixProp
---|---|---|---|---|---|---|---
W/o behavioral data | 3.849 | -0.93 | -0.69 | -1.30 | -0.75 | -0.78 | -1.05
W/o ELP charact. | 3.859 | -1.19 | -0.54 | -1.36 | -0.95 | -0.59 | -1.55
W/o frequencies | 3.880 | -1.74 | -1.38 | -1.68 | -1.88 | -1.55 | -1.87
W/o bigram AM | 3.881 | -1.78 | -2.05 | -2.32 | -1.39 | -1.94 | -1.70
W/o length feat. | 3.979 | -4.35 | -5.95 | -2.92 | -3.17 | -4.43 | -4.08
W/o position feat. | 4.095 | -7.39 | -7.68 | -4.44 | -22.88 | -7.48 | -4.30
RMSE optimization | 3.847 | -0.87 | -0.43 | 0.46 | -4.73 | -0.09 | -0.43
Default Param + MAE | 3.902 | -2.32 | -2.34 | -1.54 | -3.52 | -2.12 | -2.15
Default Param + RMSE | 4.141 | -8.59 | -7.67 | -7.65 | -12.62 | -7.43 | -8.31
Linear Regression | 4.268 | -10.64 | -9.04 | -7.88 | -24.09 | -9.47 | -8.26
LGBM on Length + Position | 4.219 | -10.63 | -10.7 | -11.4 | -8.18 | -12.1 | -10.85
Table 3: Performance (MAE) of different system versions and deviation (%) from
the best run ($MAE=3.813$). Minimum and maximum values across DVs for each row
are bolded.

For the third submission, a new round of random optimization was conducted by
evaluating parameter values close to those selected for Run 2, independently
for each DV. As it only got slightly better performance than Run 2, these
parameter values are not shown to save space.

As shown in Table 2, Runs 2 and 3 ranked at the first 2 places of the
challenge. This result was largely due to their better performance for fixProp
since the TALEP system, second in the challenge, achieved significantly better
performance for three of the five DVs, but these have less impact on the
official measurement. An analysis, carried out after the end of the challenge,
showed that the system would not have been more effective (average MAE of
3.8138) if, during the first optimization step, a specific model for each DV
had been selected.

Using Pearson’s linear correlation coefficient as a measure of effectiveness,
which is unaffected by the differences in means and variability between the
five DVs, Run 3 obtains an average r of 0.812 on the test set (min = 0.792 for
GPT; max = 0.838 for fixProp). This value is relatively high, but it can only
really be interpreted by taking into account the reliability of the average
real eye-tracking feature values.

###  4.2 Feature Usefulness

The first part of Table 3 presents the main results of an ablation procedure
aimed at examining the impact of the different types of features on the system
performance. It gives the average MAE as well as the difference in percentage
between each system and the best run for the average MAE and for the five DVs.
It must be first stressed that all features based on lemmas and POS-tag, the
two Glasgow norms and the length of the token that follows the target are
useless for predicting the test set since without them the system achieves a
MAE of 3.8134. They are thus discarded in all the ablation analyses. The
target’s positions in the sentence and the length features are clearly
essential. Among the features resulting from corpora and behavioral data, it
is the bigram association measures and the frequencies in the corpora that are
the most useful.

Generally speaking, the feature sets have comparable utility for all DVs.
However, we observe that the position in the sentences is particularly
important for predicting GPT while the length of the target is more useful for
nFix.

The second part of Table 3 presents an analysis of the utility of optimizing
the LightGBM parameters, based on the best system. Optimizing RMSE instead of
MAE is especially penalizing for GPT. Using the default values of the LightGBM
parameters is particularly penalizing when RMSE is the criterion.

A final question concerns the benefits of employing LightGBM instead of
another regression algorithm when the proposed features are used. To try to
provide at least a partial answer, I trained a multiple linear regression
model on the basis of the features used, while adding for each feature, for
which the calculation was possible, a second feature containing the logarithm
of the initial value. I replaced the missing data with 0, which is probably
not optimal. A stepwise regression procedure with a threshold to enter sets at
$p=0.01$ and a threshold to exit sets at $p=0.05$ was employed to construct
for each DV a model on the learning set and apply it to the test set. The
results obtained are given in the second to last row of Table 3. The
performances are clearly less good. It is even worse than the performance
level of a LightGBM model based only on the length and position features (see
the last row of Table 3). This regression system would have been ranked 10th
in the challenge.

##  5 Conclusion

The system proposed for the 2021 CMCL Shared Task on Eye-Tracking Data
Prediction was particularly effective, obtaining the first place in the
challenge, but it must be kept in mind that the system that came second is
superior to it for three of the five DVs. The analyzes carried out to
understand its pros and cons indicate that optimizing the LightGBM parameters
is quite beneficial to it as well as the different sets of features derived
from corpora and behavioral data, including bigram AMs which, to my knowledge,
have never been employed for this type of task.

It would have been interesting to relate these observations to the
psycholinguistic literature on the factors that influence eye fixations, but
this is unfortunately not possible here, for lack of space. In addition, this
would first require deepening the ablation analyzes by simultaneously
considering several feature sets. For instance, the lack of usefulness of the
POS-tags could simply result from the links (at least partial) between them
and the frequency and length of the tokens. Likewise, some of the bigram AMs
are relatively sensitive to the frequency of the words that compose them
(e.g., the t-score favors frequent bigrams which are usually composed of
frequent words). It is thus highly probable that some of the features in the
different sets (frequencies, behavioral data…) are redundant and can be
removed without impairing the performance of the system. This is a potential
development path.

## Acknowledgements

The author wishes to thank the organizers of this shared task for putting
together this valuable event and the reviewers for their very constructive
comments. He is a Research Associate of the Fonds de la Recherche Scientifique
- FNRS (Fédération Wallonie Bruxelles de Belgique). Computational resources
were provided by the supercomputing facilities of the UCLouvain (CISM/UCL) and
the Consortium des Equipements de Calcul Intensif en Fédération Wallonie
Bruxelles (CECI).

## References

  * Balota et al. (2007) David A. Balota, Melvin J. Yap, Keith A. Hutchison, Michael J. Cortese, Brett Kessler, Bjorn Loftis, James H. Neely, Douglas L. Nelson, Greg B. Simpson, and Rebecca Treiman. 2007.  The English lexicon project.  _Behavior Research Methods_ , 39:445–459.
  * Berry-Rogghe (1973) Godelieve L. M. Berry-Rogghe. 1973.  The computation of collocations and their relevance in lexical studies.  In Adam J Aitken, Richard W. Bailey, and Neil Hamilton-Smith, editors, _The Computer and Literary Studies_. Edinburgh University Press.
  * Bestgen (2017) Yves Bestgen. 2017.  Beyond single-word measures: L2 writing assessment, lexical richness and formulaic competence.  _System_ , 69:65–78.
  * Bestgen (2018) Yves Bestgen. 2018.  Evaluating the frequency threshold for selecting lexical bundles by means of an extension of the Fisher’s exact test.  _Corpora_ , 13:205–228.
  * Bestgen (2019) Yves Bestgen. 2019.  Evaluation de textes en anglais langue ètrangère et séries phraséologiques : comparaison de deux procédures automatiques librement accessibles.  _Revue française de linguistique appliquée_ , 24:81–94.
  * Bestgen (2021) Yves Bestgen. 2021.  LAST at SemEval-2021 Task 1: improving multi-word complexity prediction using bigram association measures.  In _Proceedings of SemEval-2021_.
  * Bestgen and Granger (2014) Yves Bestgen and Sylviane Granger. 2014.  Quantifying the development of phraseological competence in L2 English writing: An automated approach.  _Journal of Second Language Writing_ , 26:28–41.
  * Church and Hanks (1990) Kenneth Ward Church and Patrick Hanks. 1990.  Word association norms, mutual information, and lexicography.  _Computational Linguistics_ , 16(1):22–29.
  * Dunning (1993) Ted Dunning. 1993.  Accurate methods for the statistics of surprise and coincidence.  _Computational Linguistics_ , 19(1):61–74.
  * Durrant and Schmitt (2009) Philip Durrant and Norbert Schmitt. 2009.  To what extent do native and non-native writers make use of collocations?  _International Review of Applied Linguistics in Language Teaching_ , 47:157–177.
  * Evert (2009) Stefan Evert. 2009.  Corpora and collocations.  In Anke Lüdeling and Merja Kytö, editors, _Corpus Linguistics. An International Handbook_ , pages 1211–1248. Mouton de Gruyter.
  * Herdagdelen and Marelli (2017) Amac Herdagdelen and Marco Marelli. 2017.  Social media and language processing: How Facebook and Twitter provide the best frequency estimates for studying word recognition.  _Cognitive Science_ , 41:976–995.
  * Hollenstein et al. (2021) Nora Hollenstein, Emmanuele Chersoni, Cassandra Jacobs, Yohei Oseki, Laurent Prévot, and Enrico Santus. 2021.  CMCL 2021 shared task on eye-tracking prediction.  In _Proceedings of the Workshop on Cognitive Modeling and Computational Linguistics_.
  * Hollenstein et al. (2019) Nora Hollenstein, Antonio de la Torre, Nicolas Langer, and Ce Zhang. 2019.  CogniVal: A framework for cognitive word embedding evaluation.  In _Proceedings of the 23rd Conference on Computational Natural Language Learning (CoNLL)_ , pages 538–549, Hong Kong, China. Association for Computational Linguistics.
  * Hollenstein et al. (2018) Nora Hollenstein, Jonathan Rotsztejn, Marius Troendle, Andreas Pedroni, Ce Zhang, and Nicolas Langer. 2018.  ZuCo, a simultaneous EEG and eye-tracking resource for natural sentence reading.  _Scientific Data. 5:180291_ , 5(180291):1–13.
  * Hollenstein et al. (2020) Nora Hollenstein, Marius Troendle, Ce Zhang, and Nicolas Langer. 2020.  ZuCo 2.0: A dataset of physiological recordings during natural reading and annotation.  In _Proceedings of The 12th Language Resources and Evaluation Conference, LREC 2020, Marseille, France, May 11-16, 2020_ , pages 138–146. European Language Resources Association.
  * Ke et al. (2017) Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. 2017.  LightGBM: A highly efficient gradient boosting decision tree.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 3146–3154. Curran Associates, Inc.
  * Kilgarriff et al. (2014) Adam Kilgarriff, Vít Baisa, Jan Bušta, Miloš Jakubíček, Vojtěch Kovář, Jan Michelfeit, Pavel Rychlý, and Vít Suchomel. 2014.  The Sketch Engine: ten years on.  _Lexicography_ , 1(1):7–36.
  * Kliegl et al. (2004) Reinhold Kliegl, Ellen Grabner, Martin Rolfs, and Ralf Engbert. 2004.  Length, frequency, and predictability effects of words on eye movements in reading.  _European Journal of Cognitive Psychology_ , 16:262–284.
  * Kyle et al. (2018) Kristopher Kyle, Scott Crossley, and Cynthia Berger. 2018.  The tool for the automatic analysis of lexical sophistication (TAALES): version 2.0.  _Behavior Research Methods_ , 50:1030–1046.
  * Mishra and Bhattacharyya (2018) Abhijit Mishra and Pushpak Bhattacharyya. 2018.  _Applications of Eye Tracking in Language Processing and Other Areas_, pages 23–46. Springer Singapore, Singapore.
  * Pecina (2010) Pavel Pecina. 2010.  Lexical association measures and collocation extraction.  _Language Resources & Evaluation_, 44:137–158.
  * Schmid (1994) Helmutt Schmid. 1994.  Probabilistic part-of-speech tagging using decision trees.  In _International Conference on New Methods in Language Processing_ , pages 44–49.
  * Scott et al. (2019) Graham G. Scott, Anne Keitel, Marc Becirspahic, Bo Yao, and Sara C. Sereno. 2019\.  The Glasgow norms: Ratings of 5,500 words on nine scales.  _Behavior Research Methods_ , 51:1258–1270.
  * Shaoul and Chris (2006) Cyrus Shaoul and Westbury Chris. 2006.  USENET orthographic frequencies for 111,627 English words.
  * Shardlow et al. (2021) Matthew Shardlow, Richard Evans, Gustavo Paetzold, and Marcos Zampieri. 2021.  SemEval-2021 task 1: Lexical complexity prediction.  In _Proceedings of the 14th International Workshop on Semantic Evaluation (SemEval-2021)_.
  * Yarkoni et al. (2008) Tal Yarkoni, David A. Balota, and Melvin Yap. 2008.  Moving beyond Coltheart’s N: A new measure of orthographic similarity.  _Psychonomic Bulletin & Review_, 16:971–979.

##  Appendix A Appendix

At the request of a reviewer, the parameter space for the first random search
is provided below. Those for the second random search are not provided as they
did not allow to really improve the performances.

param_grid = { ’max_bin’: [16,32,48,64,80,96,112,128, 160,192,224,256],
’min_data_in_bin’:[2,3,4,5,6,8,10,12, 15,20], ’num_leaves’:
[4,5,6,7,8,9,10,11,12,13, 15,18,21,25,30], ’learning_rate’:
[0.005,0.007,0.009, 0.011,0.014,0.018,0.022,0.026,0.03, 0.035,0.05],
’min_data_in_leaf’: [2,3,4,5,6,7,8,9,10, 11,12,13,15,18,21,25,30],
’max_depth’: [3,4,5,6,7,8,9,10,11,12,13, -1], ’feature_fraction’:
list(np.linspace( 0.01, 0.90, 91)), ’bagging_freq’: list(range(3, 7, 1)),
’bagging_fraction’: list(np.linspace( 0.50, 0.90, 9))}


Next, $h_{g}^{C}$, $h_{\mu_{0}}^{C}$, and $h_{\mu_{1}}^{C}$ denote the
estimating functions used to obtain $\hat{\beta}$, $\hat{\gamma}_{0}$, and
$\hat{\gamma}_{1}$ assuming complete data. To account for incomplete sampling
of $S$, we use inverse probability sampling (IPS) weighting, which gives us
our final estimating equations:

$0=\sum_{i=1}^{n}h_{g}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\beta,\alpha)$
$0=\sum_{i=1}^{n}h_{\mu_{0}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{0})$
$0=\sum_{i=1}^{n}h_{\mu_{1}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{1})$

where

$h_{g}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\beta,\alpha):=\frac{Z_{i}\epsilon_{Si}h_{g}^{C}(X_{i},S_{i},\tilde{T}_{i},\Delta_{i};\beta)}{\pi(X_{i},Z_{i},A_{i},\tilde{T}_{i},\Delta_{i};\alpha)}$
$h_{\mu_{0}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{0}):=\frac{(1-Z_{i})(1-A_{i})\epsilon_{Si}h_{\mu_{0}}^{C}(g(X_{i},S_{i},\tilde{T}_{i},\Delta_{i};\beta),X_{i};\gamma_{0})}{\pi(X_{i},Z_{i},A_{i},\tilde{T}_{i},\Delta_{i};\alpha)}$
$h_{\mu_{1}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{1}):=\frac{(1-Z_{i})A_{i}\epsilon_{Si}h_{\mu_{1}}^{C}(g(X_{i},S_{i},\tilde{T}_{i},\Delta_{i};\beta),X_{i};\gamma_{1})}{\pi(X_{i},Z_{i},A_{i},\tilde{T}_{i},\Delta_{i};\alpha)}.$

The plug in estimators $\widehat{\theta}_{0,\mbox{\text{plug-in}}}$ and
$\widehat{\theta}_{1,\mbox{\text{plug-in}}}$ can be written as the solutions
$\phi_{0}$, $\phi_{1}$, respectively, of the following estimating equations:

$0=\sum_{i=1}^{n}(1-Z_{i})\\{\mu_{0}(X_{i};\gamma_{0})-\phi_{0}\\}$
$0=\sum_{i=1}^{n}(1-Z_{i})\\{\mu_{1}(X_{i};\gamma_{1})-\phi_{1}\\}.$

We define the parameter vector
$\nu:=(\beta,\gamma_{0},\gamma_{1},\alpha,\phi_{0},\phi_{1})$. Our final set
of stacked estimating functions are

$h_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\nu)=\left(\begin{array}[]{c}h_{\pi}(\epsilon_{Si},X_{i},Z_{i},A_{i},\tilde{T}_{i},\Delta_{i};\alpha)\\\
h_{g}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\beta,\alpha)\\\
h_{\mu_{0}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{0})\\\
h_{\mu_{1}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\alpha,\beta,\gamma_{1})\\\
(1-Z_{i})\\{\mu_{0}(X_{i};\gamma_{0})-\phi_{0}\\}\\\
(1-Z_{i})\\{\mu_{1}(X_{i};\gamma_{1})-\phi_{1}\\}\end{array}\right)$

along with the stacked estimating equation

$0=\sum_{i=1}^{n}h_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\nu).$

From this estimating equation, we obtain an estimator $\hat{\nu}$. Now, we
define the following:

$\displaystyle h^{\prime}_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\nu)$
$\displaystyle=\partial h_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\nu)/\partial\nu,$
$\displaystyle
W_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\\{-h^{\prime}_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu})\right\\}$
$\displaystyle
Q_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}h_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu})h_{\text{stack
}}(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu})^{T}.$

The sandwich variance estimator for $\hat{\nu}$ is

$\displaystyle
V_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)=$
$\displaystyle
W_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)^{-1}Q_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)$
$\displaystyle\left\\{W_{n}\left(X_{i},S_{i},Z_{i},A_{i},\epsilon_{Si},\tilde{T}_{i},\Delta_{i};\hat{\nu}\right)^{-1}\right\\}^{T}.$

##  10 Appendix D: Simulation study details

###  10.1 Design of simulation study conditions

The simulation study was designed to roughly match published characteristics
of GBS and its risk factors.

  1. 1.

Probability of IGbsD by 90 days of age for colonized mothers $\approx 0.005$.
From Vekemans et al. [2019b], approximately 20% of mothers are colonized and
0.1% of all infants have IGbsD. We obtain an approximate incidence of $0.005$
by assuming that only infants of colonized mothers can develop invasive GBS
disease.

  2. 2.

Geometric mean (95% CI) of cord-blood IgG concentration 0.01 (0.01-0.02) in
IGbsD cases and 0.04 (0.03 - 0.06) in controls, which was observed for infant
RibN IgG in Dangor et al. [2023].

  3. 3.

Covariates $X_{1},X_{2},X_{3}$, where $X_{1}$ represents pre-term birth (less
than 37 weeks gestational age), $X_{2}$ represents maternal age (younger age a
risk factor), and $X_{3}$ represents a continuous covariate unrelated to the
outcome.

$X_{1}|Z,A\sim\text{Bernoulli}(0.05),X_{2}|Z,A\sim\text{Uniform}(18,40),X_{3}|Z,A\sim\text{Normal}(0,1)$

     * •

Pre-term birth is associated with roughly $2.56\times$ odds of early onset
IGbsD (Puopolo et al., 2011)

     * •

Maternal age $<25$ years is associated with roughly $1.94\times$ odds of IGbsD
(Parente et al., 2017)

Based on these constraints, we set the following to be our true $A=0$ data
generating conditional regression function:

$P(Y=1|S,X_{1},X_{2},X_{3},A=0,Z)=\beta_{0}+\beta_{1}S+\beta_{2}X_{1}+\beta_{3}X_{2}+\beta_{4}X_{3}$

where $\beta_{0}=-17.1$, $\beta_{1}=-8.2$, $\beta_{2}=0.69$,
$\beta_{3}=-0.03$, $\beta_{4}=0$, and we set our distribution of
$S\sim\text{Normal}(-1.45,0.0225)$, where $S$ represents the log IgG
biomarker, $X_{1}\sim\text{Bernoulli}(0.05)$,
$X_{2}\sim\text{Uniform}(18,40)$, and $X_{3}\sim\text{Normal}(0,1)$.

With this data-generating function, we have that in our observational study
(using empirical measures from a simulated $n=$10,000,000 size data set):

  * •

Incidence of IGbsD by 90 days $\approx 0.005$

  * •

Geometric mean cord-blood IgG biomarker in IGbsD cases is 0.01 and in controls
is 0.04

  * •

Pre-term birth associated with $2\times$ odds of disease

  * •

Each one-year increase in maternal age associated with 3% lower odds of
disease

###  10.2 Results from the larger sample simulation study

In this section, we report results for a simulation study with a higher case
rate. We set the following to be our true $A=0$ data generating conditional
regression function:

$P(Y=1|X_{1},X_{2},X_{3},A=0,Z,S)=\beta_{0}+\beta_{1}S+\beta_{2}X_{1}+\beta_{3}X_{2}+\beta_{4}X_{3}$

where $\beta_{0}=-14$, $\beta_{1}=-7$, $\beta_{2}=0.69$, $\beta_{3}=-0.03$,
$\beta_{4}=0$.

We preserve the same distribution of $S$ as in Simulation Study 1, with
$S\sim\text{Normal}(-1.45,0.225)$, resulting in a baseline case rate of 0.016
(about 3 times higher than 0.005).

To generate VE $=\\{0,0.5,0.9\\}$, we manipulate the distribution of the
biomarker in the vaccine arm $S|A=1,Z=0$:

  * •

To create VE $=0$, set $S|A=1,Z=0\sim\text{Normal}(-1.45,0.0225)$

  * •

To create VE $=0.5$, set $S|A=1,Z=0\sim\text{Normal}(-1.29,0.04)$

  * •

To create VE $=0.9$, set $S|A=1,Z=0\sim\text{Normal}(-1.04,0.0441)$.

Table 3: Results for the Supplemental Simulation Study with approximately 3 times higher IGbsD outcome rate on empirical bias, average standard error (sandwich and bootstrap), standard deviation of estimates, and 95% confidence interval coverage of $E[Y(0)|Z=0],E[Y(1)|Z=0]$, and VE. SE (bs) = bootstrap standard error, SE (sw) = sandwich standard error, SD = standard deviation of estimates, Cov = 95% confidence interval coverage using sandwich standard error, Sampled = number with $S$ data per arm in phase 3 study. The standard errors for $\log(1-\widehat{{\rm VE}})$ are displayed for VE. | Bias | SE (bs) | SE (sw) | SD | Cov
---|---|---|---|---|---
True VE = 0 |  |  |  |  |
$E[Y(0)|Z=0]$ | 0 | 0.000592 | 0.000846 | 0.000781 | 0.97
$E[Y(1)|Z=0]$ | 0 | 0.000592 | 0.000844 | 0.000793 | 0.97
VE | -0.000535 | 0.0338 | 0.0335 | 0.0338 | 0.96
True VE = 0.5 |  |  |  |  |
$E[Y(0)|Z=0]$ | 0 | 0.000596 | 0.000845 | 0.000777 | 0.98
$E[Y(1)|Z=0]$ | -6.00E-05 | 0.000367 | 0.000473 | 0.000459 | 0.96
VE | 0.004 | 0.0487 | 0.0482 | 0.0487 | 0.95
True VE = 0.9 |  |  |  |  |
$E[Y(0)|Z=0]$ | 0 | 6.00E-04 | 0.00085 | 0.000832 | 0.96
$E[Y(1)|Z=0]$ | -1.00E-05 | 0.000157 | 0.000167 | 0.000163 | 0.96
VE | 0.001 | 0.11 | 0.109 | 0.108 | 0.96

Finally, we sample all 6,200 participants in the phase 3 study ($Z=1$) in this
supplemental simulation study.

Results are shown in Table 3. In this higher case rate scenario, we observe
minimal bias in estimating VE, approximately nominal confidence interval
coverage, and close agreement between bootstrap, sandwich, and empirical
standard errors.