# Spontaneous collapse by entanglement suppression

Eyal Buks<EMAIL_ADDRESS>Andrew and Erna Viterbi Department of
Electrical Engineering, Technion, Haifa 32000, Israel

###### Abstract

We study a recently proposed modified Schrödinger equation having an added
nonlinear term, which gives rise to disentanglement. The process of quantum
measurement is explored for the case of a pair of coupled spins. We find that
the deterministic time evolution generated by the modified Schrödinger
equation mimics the process of wavefunction collapse. Added noise gives rise
to stochasticity in the measurement process. Conflict with both principles of
causality and separability can be avoided by postulating that the nonlinear
term is active only during the time when subsystems interact. Moreover, in the
absence of entanglement, all predictions of standard quantum mechanics are
unaffected by the added nonlinear term.

Introduction \- In standard quantum mechanics a measurement is described by a
two-step process. The first step is governed by the standard Schrödinger
equation. To avoid a possible paradoxical outcome of a description based only
on the first step (undefined cat state Schrodinger_807 ), a second step is
postulated, in which the state vector collapses. However, it has remained
unknown how such a second step can be self-consistently added Penrose_4864 ;
Leggett_939 ; Leggett_022001 . This difficulty has became known as the problem
of quantum measurement.

In this work we explore an alternative to the collapse postulate, which is
based on a modified Schrödinger equation that has an added nonlinear term
giving rise to disentanglement Buks_355303 ; Buks_025302 . The proposed
equation can be constructed for any physical system whose Hilbert space has
finite dimensionality, and it does not violate norm conservation of the time
evolution. We explore the dynamics of a system made of two coupled spins, and
find that disentanglement gives rise to a process similar to state vector
collapse.

Other types of nonlinear extensions of quantum mechanics Geller_2200156  have
been previously proposed and studied Weinberg_336 ; Weinberg_61 ; Doebner_397
; Doebner_3764 ; Gisin_5677 ; Kaplan_055002 ; Munoz_110503 . Most previously
proposed extensions give rise to a spontaneous collapse Bassi_471 ; Pearle_857
; Ghirardi_470 ; Bassi_257 ; Arnquist_080401 . In some cases, however, the
proposed nonlinear models are inconsistent with well-established physical
principles. Moreover, many predictions of standard quantum mechanics, that
have been experimentally verified to very high precision, are significantly
altered by some of the proposed nonlinear extensions. Such difficulties are
discussed below in the final part of this paper for the case of our proposed
modified Schrödinger equation. We find that possible conflicts with the
principles of causality and separability, and with many experimentally
confirmed predictions of standard quantum mechanics, can be avoided by
postulating that disentanglement is active only when subsystems interact.

Figure 1: Dipolar measurement. The spin numbers are $S_{1}=1/2$ and
$S_{2}=21/2$, the rates are $\gamma=\omega_{\mathrm{d}}=1$, and
$\mathbf{\hat{n}}_{2}=-\mathbf{\hat{z}}$ (initial direction of the $S_{2}$
spin, which is labeled by a cyan star symbols). For the plots labeled by the
numbers 1, 2 and 3 the dipolar unit vector $\mathbf{\hat{u}}_{\mathrm{d}}$ is
given by $\mathbf{\hat{u}}_{\mathrm{d}}=\mathbf{\hat{x}}$ (i.e.
$\mathbf{\hat{u}}_{\mathrm{d}}$ is perpendicular to $\mathbf{\hat{n}}_{2}$),
whereas
$\mathbf{\hat{u}}_{\mathrm{d}}=\left(\sin\left(3\pi/8\right)\cos\left(3\pi/4\right),\sin\left(3\pi/8\right)\sin\left(3\pi/4\right),\cos\left(3\pi/8\right)\right)$
for the plots labeled by the number 4 (i.e.
$\mathbf{\hat{n}}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\neq 0$). At time $t=0$
the spin 1/2 is pointing in the direction
$\mathbf{\hat{n}}_{1}=\left(\sin\theta_{1}\cos\varphi_{1},\sin\theta_{1}\sin\varphi_{1},\cos\theta_{1}\right)$,
where for (1) $\theta_{1}=0.55\pi$ and $\varphi_{1}=0.45\pi$, for (2)
$\theta_{1}=0.55\pi$ and $\varphi_{1}=0.55\pi$, for (3) $\theta_{1}=0.55\pi$
and $\varphi_{1}=0.75\pi$, and for (4) $\theta_{1}=0.5\pi$, and
$\varphi_{1}=0.5\pi$. Red star symbols label the initial points
$\mathbf{\hat{n}}_{1}$, and the blue solid (dashed) lines connect the origin
and the unit vectors $\mathbf{\hat{u}}_{\mathrm{d}}$
($-\mathbf{\hat{u}}_{\mathrm{d}}$). Figure 2: Basins of attraction. All
parameters are the same as those used to generate the plots of Fig. 1 labeled
by the number 4. Initial direction of the $S_{2}$ spin
$\mathbf{\hat{n}}_{2}=-\mathbf{\hat{z}}$ is labeled by a cyan tick, and the
dipolar unit vector
$\mathbf{\hat{u}}_{\mathrm{d}}=\left(\sin\left(3\pi/8\right)\cos\left(3\pi/4\right),\sin\left(3\pi/8\right)\sin\left(3\pi/4\right),\cos\left(3\pi/8\right)\right)$
is labeled by a blue tick. At time $t=0$ the spin 1/2 is pointing in the
direction $\mathbf{\hat{n}}_{1}$. The yellow (purple) colored region is the
basin of attraction lying in the hemisphere
$\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}>0$
($\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}<0$), and the
corresponding attractor is $\mathbf{\hat{u}}_{\mathrm{d}}$
($-\mathbf{\hat{u}}_{\mathrm{d}}$).

Disentanglement \- Consider a system composed of two subsystems labeled as ’1’
and ’2’, respectively. The dimensionality of the Hilbert spaces of both
subsystems, which is denoted by $N_{1}$ and $N_{2}$, respectively, is assumed
to be finite. The system is in a normalized pure state vector
$\left|\psi\right\rangle$ given by

$\left|\psi\right\rangle=\mathcal{K}_{1}C\otimes\mathcal{K}_{2}^{\mathrm{T}}\;,$
(1)

where $C$ is a $N_{1}\times N_{2}$ matrix having entries $C_{k_{1},k_{2}}$,
matrix transposition is denoted by $\mathrm{T}$,
$\mathcal{K}_{1}=\left(\left|k_{1}\right\rangle_{1},\left|k_{2}\right\rangle_{1},\cdots,\left|k_{N_{1}}\right\rangle_{1}\right)$,
$\mathcal{K}_{2}=\left(\left|k_{1}\right\rangle_{2},\left|k_{2}\right\rangle_{2},\cdots,\left|k_{N_{2}}\right\rangle_{2}\right)$,
and $\left\\{\left|k_{1}\right\rangle_{1}\right\\}$
($\left\\{\left|k_{2}\right\rangle_{2}\right\\}$) is an orthonormal basis
spanning the Hilbert space of subsystem ’1’ (’2’).

The purity $P_{1}$ ($P_{2}$) is defined by
$P_{1}=\operatorname{Tr}\rho_{1}^{2}$ ($P_{2}=\operatorname{Tr}\rho_{2}^{2}$),
where $\rho_{1}=\operatorname{Tr}_{2}\rho$
($\rho_{2}=\operatorname{Tr}_{1}\rho$) is the reduced density operator of the
first (second) subsystem. By employing the Schmidt decomposition one finds
that $P_{1}=P_{2}\equiv P$, where
$P=1-\left\langle\mathcal{Q}\right\rangle=1-\left\langle\psi\right|\mathcal{Q}\left|\psi\right\rangle$,
the operator $\mathcal{Q}$ is given by [see Eq. (27) of appendix A, and Ref.
Buks_QMLN ]

$\mathcal{Q}=\frac{1}{2}{\displaystyle\sum\limits_{k_{1}^{\prime}<k_{1}^{\prime\prime}}}{\displaystyle\sum\limits_{k_{2}^{\prime}<k_{2}^{\prime\prime}}}\left|\Psi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right\rangle\left\langle\Psi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|\;,$
(2)

and the state
$\left\langle\Psi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|$,
which depends on the matrix $C$ corresponding to a given state
$\left|\psi\right\rangle$, is given by (note that
$\left\langle\Psi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|$
is not normalized)

$\left\langle\Psi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|=C_{{}_{\mathrm{d}}}\left\langle\mathrm{a}\right|+C_{{}_{\mathrm{a}}}\left\langle\mathrm{d}\right|-C_{\mathrm{c}}\left\langle\mathrm{b}\right|-C_{\mathrm{b}}\left\langle\mathrm{c}\right|\;,$
(3)

where $\mathrm{a}=k_{1}^{\prime},k_{2}^{\prime}$,
$\mathrm{b}=k_{1}^{\prime},k_{2}^{\prime\prime}$,
$\mathrm{c}=k_{1}^{\prime\prime},k_{2}^{\prime}$ and
$\mathrm{d}=k_{1}^{\prime\prime},k_{2}^{\prime\prime}$. Note that
$\left\langle\mathcal{Q}\right\rangle=0$ for a product state. In standard
quantum mechanics $\left\langle\mathcal{Q}\right\rangle$ is time independent
when the subsystems are decoupled (i.e. their mutual interaction vanishes).

As an example, consider a two spin 1/2 system (i.e. $N_{1}=N_{2}=2$) in a pure
state $\left|\psi\right\rangle$ given by
$\left|\psi\right\rangle=a\left|--\right\rangle+b\left|-+\right\rangle+c\left|+-\right\rangle+d\left|++\right\rangle$.
For this case the sum in Eq. (2) contains a single term with
$\left\langle\Psi\right|=d\left\langle-,-\right|-c\left\langle-,+\right|-b\left\langle+,-\right|+a\left\langle+,+\right|$,
and thus $P=1-2\left|ad-bc\right|^{2}$. Note that for this case
$\left\langle\mathcal{Q}\right\rangle\leq 1/2$ (provided that
$\left|\psi\right\rangle$ is normalized) Wootters_2245 .

Consider a modified Schrödinger equation for the ket vector
$\left|\psi\right\rangle$ having the form

$\frac{\mathrm{d}}{\mathrm{d}t}\left|\psi\right\rangle=\left[-i\hbar^{-1}\mathcal{H}-\gamma\left(\mathcal{Q}-\left\langle\mathcal{Q}\right\rangle\right)\right]\left|\psi\right\rangle\;,$
(4)

where $\hbar$ is the Planck’s constant, $\mathcal{H}=\mathcal{H}^{{\dagger}}$
is the Hamiltonian, the rate $\gamma$ is positive, and the operator
$\mathcal{Q}$ is given by Eq. (2). The added nonlinear term proportional to
$\gamma$ gives rise to disentanglement, however, it has no effect when
$\left|\psi\right\rangle$ represents a product state. Note that the norm
conservation condition
$0=\left(\mathrm{d}/\mathrm{d}t\right)\left\langle\psi\right.\left|\psi\right\rangle$
is satisfied by the modified Schrödinger equation (4).

Dipolar interaction \- As an example, the dynamics generated by the modified
Schrödinger equation (4) is explored for the case of dipolar interaction
between two spins having spin quantum numbers $S_{1}$ and $S_{2}$,
respectively. The dipolar interaction is represented by the operator
$V_{\mathrm{d}}=\hbar^{-1}\omega_{\mathrm{d}}\left(\mathbf{S}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right)\left(\mathbf{S}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right)$,
where the rate $\omega_{\mathrm{d}}$ is positive, $\mathbf{S}_{n}$ is the spin
angular momentum vector operator of the $n$’th spin
($n\in\left\\{1,2\right\\}$), and
$\mathbf{\hat{u}}_{\mathrm{d}}=\left(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta\right)$
is a unit vector.

Time evolution examples for the case $S_{1}=1/2$ and $S_{2}=21/2$ are shown by
the plots in Fig. 1. The initial state at time $t=0$ is a product state, for
which the spin 1/2 is pointing in the direction of the unit vector
$\mathbf{\hat{n}}_{1}$ (labeled by a red star symbol), and the spin 21/2 is
pointing in the direction of the unit vector
$\mathbf{\hat{n}}_{2}=-\mathbf{\hat{z}}$ (labeled by a cyan star symbol). The
overlaid blue solid (dashed) lines connect the origin and the dipolar coupling
unit vectors $\mathbf{\hat{u}}_{\mathrm{d}}$
($-\mathbf{\hat{u}}_{\mathrm{d}}$). The spin 1/2 Bloch vector
$\mathbf{k}=\left(\hbar/2\right)^{-1}\left\langle\mathbf{S}_{1}\right\rangle$
is numerically calculated by integrating the modified Schrödinger equation (4)
for the case $\mathcal{H}=V_{\mathrm{d}}$. The black solid lines in Fig.
1(a1), (a2), (a3) and (a4) represent the spin 1/2 Bloch vector $\mathbf{k}$
evolving from its initial value $\mathbf{\hat{n}}_{1}$ at time $t=0$. The
single-spin purity $P=1-\left\langle\mathcal{Q}\right\rangle$ as a function of
time $t$ is shown in Fig. 1(b1), (b2), (b3) and (b4).

For the plots in Fig. 1 labeled by the numbers 1, 2 and 3, the dipolar unit
vector $\mathbf{\hat{u}}_{\mathrm{d}}$ is given by
$\mathbf{\hat{u}}_{\mathrm{d}}=\mathbf{\hat{x}}$ (i.e.
$\mathbf{\hat{u}}_{\mathrm{d}}$ is perpendicular to
$\mathbf{\hat{n}}_{2}=-\mathbf{\hat{z}}$). These plots, which differ by the
initial direction $\mathbf{\hat{n}}_{1}$ of the spin 1/2 (labeled by red star
symbols), demonstrate that the Bloch sphere is divided into two basins of
attraction. The first (second) basin is the hemisphere
$\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}>0$
($\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}<0$), and the
corresponding attractor is $\mathbf{\hat{u}}_{\mathrm{d}}$
($-\mathbf{\hat{u}}_{\mathrm{d}}$).

While $\mathbf{\hat{n}}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}=0$ for the plots
in Fig. 1 labeled by the numbers 1, 2 and 3, the behavior when the initial
spin $S_{2}$ direction $\mathbf{\hat{n}}_{2}$ is not perpendicular to the
dipolar coupling unit vector $\mathbf{\hat{u}}_{\mathrm{d}}$ is demonstrated
by the plots labeled by the number 4. The plot in Fig. 1(a4) shows that the
Bloch vector trajectory, from the initial value $\mathbf{\hat{n}}_{1}$
(labeled by the red star symbol) towards the attractor at
$\mathbf{\hat{u}}_{\mathrm{d}}$ becomes spiral-like when
$\mathbf{\hat{n}}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\neq 0$. The basins of
attraction for this case (i.e. plots in Fig. 1 labeled by the number 4) are
shown in Fig. 2. This example demonstrates that the dipolar unit vector
$\mathbf{\hat{u}}_{\mathrm{d}}$ determines the spin 1/2 component that is
being measured. The measurement process is deterministic however the outcome,
which is either $+1$ (when
$\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}>0$) or $-1$ (when
$\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}<0$) is quantized. This
behavior is demonstrated by the green dash-dotted line in Fig. 3, in which the
probability $p_{+}$ that the measurement outcome is $+1$ is plotted as a
function of the angle
$\theta_{1}=\cos^{-1}\left(\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right)$.
For comparison, the red solid line represents the Born rule of standard
quantum mechanics, for which
$p_{+}\left(\theta_{1}\right)=\cos^{2}\left(\theta_{1}/2\right)$. A simplified
model is employed below to explore noise-induced stochasticity.

Noise \- The effect of external noise is taken into account by applying a
random rotation to the initial spin 1/2 Block vector $\mathbf{\hat{n}}_{1}$.
The random rotation is characterized by an axis normal to
$\mathbf{\hat{n}}_{1}$, and by a rotation angle $\phi_{\mathrm{r}}$. As an
example, consider the case where the rotation angle $\phi_{\mathrm{r}}$ has a
wrapped Cauchy probability distribution $f\left(\phi_{\mathrm{r}}\right)$
given by

$f\left(\phi_{\mathrm{r}}\right)=\frac{1}{2\pi}\frac{\sinh\phi_{0}}{\cosh\phi_{0}-\cos\phi_{\mathrm{r}}}\;,$
(5)

where $\phi_{0}>0$ is a scale factor. Consider a rotated frame, in which the
dipolar unit vector $\mathbf{\hat{u}}_{\mathrm{d}}$ is parallel to the unit
vector $\mathbf{\hat{z}}$. The unit vector $\mathbf{\hat{n}}_{1}$ in this
frame is denoted by $\mathbf{\hat{n}}_{1\mathrm{R}}$. The probability $p_{+}$
that the measurement outcome is $+1$ is calculated by spherical integration
over the hemisphere $z^{\prime}\geq 0$

$p_{+}=\frac{1}{4\pi}\int_{0}^{\pi/2}\mathrm{d}\theta^{\prime}\;\sin\theta^{\prime}\int_{0}^{2\pi}\mathrm{d}\varphi^{\prime}\;\frac{4f\left(\theta_{1\mathrm{R}}\right)}{\sin\theta_{1\mathrm{R}}}\;,$
(6)

where
$\theta_{1\mathrm{R}}=\cos^{-1}\left(\mathbf{\hat{n}}_{1\mathrm{R}}\cdot\mathbf{\hat{n}}^{\prime}\right)$,
and where
$\mathbf{\hat{n}}^{\prime}=\left(\sin\theta^{\prime}\cos\varphi^{\prime},\sin\theta^{\prime}\sin\varphi^{\prime},\cos\theta^{\prime}\right)$.
As can be seen from the blue dashed line in Fig. 3, which is calculated using
Eq. (6) with a scale factor of $\phi_{0}=0.5$, noise-induced stochasticity
mimics the behavior predicted by the Born rule (red solid line).

Figure 3: Noise. The probability $p_{+}$ is plotted as a function of the angle
$\theta_{1}=\cos^{-1}\left(\mathbf{\hat{n}}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right)$
for the noiseless case (green dash-dotted line), the case $\phi_{0}=0.5$ (blue
dashed line), and the Born rule (red solid line).

The measurement time \- For the examples shown in Fig. 1, initially at time
$t=0$, the ket vector $\left|\psi\right\rangle$ represents a product state
having single-spin purity $P=1$. The time dependency of $P$ is shown in Fig.
1(b1), (b2), (b3) and (b4). In the short time limit of
$\omega_{\mathrm{d}}t\ll 1$ the effect of the disentanglement term in the
modified Schrödinger equation (4) is relatively weak (since
$\left\langle\mathcal{Q}\right\rangle$ is initially small), and consequently
$P$ rapidly drops due to entanglement generated by the dipolar interaction
$V_{\mathrm{d}}$. At latter times, when disentanglement becomes sufficiently
efficient, the single-spin purity $P$ starts increasing. Interaction-induced
generation of entanglement becomes inefficient when the spin 1/2 becomes
nearly parallel or nearly anti-parallel to the dipolar unit vector
$\mathbf{\hat{u}}_{\mathrm{d}}$, and consequently the single-spin purity $P$
approaches unity in the long time limit.

For sufficiently short times after turning on the interaction (i.e. after
$t=0$), time evolution is dominated by the effect of the dipolar interaction.
When the effect of the disentanglement term is disregarded, one finds that in
the short time limit the following holds
$\mathrm{d}\left\langle\mathbf{S}_{n}\right\rangle/\mathrm{d}t\simeq\omega_{n}\mathbf{\hat{u}}_{\mathrm{d}}\times\left\langle\mathbf{S}_{n}\right\rangle$,
where $n\in\left\\{1,2\right\\}$,
$\omega_{1}=\omega_{\mathrm{d}}\hbar^{-1}\left\langle\mathbf{S}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right\rangle$
and
$\omega_{2}=\omega_{\mathrm{d}}\hbar^{-1}\left\langle\mathbf{S}_{1}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right\rangle$.
Thus, in the short time limit, the purity $P$ is roughly given by $P\simeq
1-\left(2^{-3/2}S_{2}\left|\mathbf{\hat{n}}_{1}\times\mathbf{\hat{u}}_{\mathrm{d}}\right|\left(\mathbf{\hat{n}}_{2}\cdot\mathbf{\hat{u}}_{\mathrm{d}}\right)\omega_{\mathrm{d}}t\right)^{2}$
[see Eqs. (6.192) and (8.701) of Ref. Buks_QMLN , and note that it is assumed
that in the short time limit the spin states are nearly spin coherent states
Radcliffe_313 ]. The above-derived expression for the purity time evolution
$P\left(t\right)$ reveals the dependence of short-time dynamics on the
macroscopicity of the measuring apparatus (i.e. the second spin), which is
represented by the spin number $S_{2}$.

Vanishing Hamiltonian \- To gain further insight into the disentanglement
process generated by the nonlinear term
$-\gamma\left(\mathcal{Q}-\left\langle\mathcal{Q}\right\rangle\right)$ added
to the Schrödinger equation (4), consider for simplicity the case where the
Hamiltonian vanishes, i.e. $\mathcal{H}=0$. The Schmidt decomposition of a
general state vector $\left|\psi\right\rangle$ is expressed as

$\left|\psi\right\rangle={\displaystyle\sum\limits_{l=1}^{\min\left(N_{1},N_{2}\right)}}q_{l}\left|l,l\right\rangle\;,$
(7)

where $q_{l}$ are non-negative real numbers, the tensor product$\
\left|l\right\rangle_{1}\otimes\left|l\right\rangle_{2}$ is denoted by
$\left|l,l\right\rangle\ $, and $\left\\{\left|l\right\rangle_{1}\right\\}$
($\left\\{\left|l\right\rangle_{2}\right\\}$) is an orthonormal basis spanning
the Hilbert space of subsystem ’1’ (’2’). Note that for a product state
$q_{l}=\delta_{l,l_{0}}$, where
$l_{0}\in\left\\{1,2,\cdots,\min\left(N_{1},N_{2}\right)\right\\}$. The
normalization condition reads
$\left\langle\psi\right.\left|\psi\right\rangle=L_{2}=1$, where the $n$’th
moment $L_{n}$ is defined by

$L_{n}={\displaystyle\sum\limits_{l=1}^{\min\left(N_{1},N_{2}\right)}}q_{l}^{n}\;.$
(8)

Note that for a product state $L_{n}=1$ for any positive integer $n$ (provided
that $\left|\psi\right\rangle$ is normalized).

Figure 4: Vanishing Hamiltonian. The plot shows an example solution of the set
of equations (10) for the case $\min\left(N_{1},N_{2}\right)=10$ and
$\gamma=1$. The solution for $q_{l_{0}}\left(t\right)$ is represented by the
red line, whereas the blue lines represent the solutions for
$q_{l}\left(t\right)$ with $l\neq l_{0}$. For this example,
$q_{l}\left(t=0\right)\simeq\left(\min\left(N_{1},N_{2}\right)\right)^{-1/2}$,
i.e. the initial value of the purity $P$ is close to its smallest possible
value of $1/\min\left(N_{1},N_{2}\right)$. The corresponding initial entropy
$\sigma$ is close to its largest possible value of
$\log\left(\min\left(N_{1},N_{2}\right)\right)$. In the limit
$t\rightarrow\infty$ the purity $P\rightarrow 1$ (largest possible value) and
the entropy $\sigma\rightarrow 0$ (smallest possible value).

In the Schmidt basis, the following holds [see Eqs. (2) and (3)]

$\mathcal{Q}\left|\psi\right\rangle={\displaystyle\sum\limits_{l=1}^{\min\left(N_{1},N_{2}\right)}}q_{l}\left(1-q_{l}^{2}\right)\left|l,l\right\rangle\;,$
(9)

and $\left\langle\mathcal{Q}\right\rangle=1-L_{4}$, and thus [see Eq. (4)]

$\frac{\mathrm{d\log}q_{l}}{\mathrm{d}t}=\gamma\left(q_{l}^{2}-L_{4}\right)\;.$
(10)

An example solution of the set of equations (10) for the case
$\min\left(N_{1},N_{2}\right)=10$ and $\gamma=1$ is shown in Fig. 4.

The time evolution of the $n$’th moment $L_{n}$ is governed by [see Eqs. (8)
and (10)]

$\frac{\mathrm{d}L_{n}}{\mathrm{d}t}=n\gamma\left(L_{n+2}-L_{n}L_{4}\right)\;.$
(11)

For the case of $n=2$, Eq. (11) yields the norm conservation condition
$0=\left(\mathrm{d}/\mathrm{d}t\right)\left\langle\psi\right.\left|\psi\right\rangle$,
which is satisfied provided that $\left|\psi\right\rangle$ is normalized, i.e.
$L_{2}=1$ [see Eq. (7)]. For the case $n=4$ Eq. (11) yields an evolution
equation for the purity $P=L_{4}$, which is given by
$\mathrm{d}L_{4}/\mathrm{d}t=4\gamma\left(L_{6}-L_{4}^{2}\right)$. Using the
Cauchy–Schwarz inequality one finds that $L_{4}^{2}\leq L_{2}L_{6}$ [see Eq.
(8)], hence $\mathrm{d}P/\mathrm{d}t\geq 0$ (recall the normalization
condition $L_{2}=1$), i.e. the purity $P$ monotonically increases with time.
The same conclusion can alternatively be drawn from Eq. (10), which can be
expressed as $\mathrm{d}q_{l}/\mathrm{d}t=\partial H/\partial q_{l}$, where
$H=\left(\gamma/4\right)\left(3-2L_{2}\right)L_{4}$ [see Eq. (8), and note
that $H=\left(\gamma/4\right)L_{4}=\left(\gamma/4\right)P$ when $L_{2}=1$].

For any two integers
$l^{\prime},l^{\prime\prime}\in\left\\{1,2,\cdots,\min\left(N_{1},N_{2}\right)\right\\}$
the following holds [see Eq. (10)]

$\frac{\mathrm{d\log}\frac{q_{l^{\prime}}}{q_{l^{\prime\prime}}}}{\mathrm{d}t}=\gamma\left(q_{l^{\prime}}^{2}-q_{l^{\prime\prime}}^{2}\right)\;.$
(12)

The above relation (12) implies that the ratio
$q_{l^{\prime}}/q_{l^{\prime\prime}}$ monotonically increases with time,
provided that $q_{l^{\prime}}>q_{l^{\prime\prime}}$ (recall that $\gamma>0$).
This behavior gives rise to disentanglement. Consider the case where
initially, at time $t=0$, $q_{l_{0}}=\max\left\\{q_{l}\right\\}$ for a unique
positive integer
$l_{0}\in\left\\{1,2,\cdots,\min\left(N_{1},N_{2}\right)\right\\}$. For this
case, $\left|\psi\right\rangle$ evolves into the product state
$\left|l_{0},l_{0}\right\rangle$ in the long time limit, i.e.
$q_{l}\rightarrow\delta_{l,l_{0}}$ in the limit $t\rightarrow\infty$ (see Fig.
4). Note, however, that in the long time limit the state can be strongly
affected by noise when initially the set $\left\\{q_{l}\right\\}$ doesn’t have
a unique member significantly larger than all others.

Discussion \- As was already mentioned above, several types of nonlinear
extensions of quantum mechanics have been proposed and explored Bassi_471 ;
Bennett_170502 ; Kowalski_1 ; Fernengel_385701 ; Kowalski_167955 . However, it
was found that for some cases, the proposed nonlinear extension gives rise to
the violation of the causality principle by enabling superluminal signaling
Bassi_055027 ; Jordan_022101 ; Polchinski_397 ; Helou_012021 . More recently,
it was shown that when a condition called ’convex quasilinearity’ is satisfied
by a given nonlinear master equation, the violation of the causality principle
becomes impossible Rembielinski_012027 ; Rembielinski_420 . Some of the
proposed nonlinear extensions are inconsistent with the principle of
separability Hejlesen_thesis ; Jordan_022101 ; Jordan_012010 . Moreover, any
proposed extension must be ruled out if it alters predictions of standard
quantum mechanics that have been experimentally confirmed.

The modified Schrödinger equation given by Eq. (4) has an important advantage
compared to other proposals: the added nonlinear term
$-\gamma\left(\mathcal{Q}-\left\langle\mathcal{Q}\right\rangle\right)$ has no
effect on product states. This implies that in the absence of entanglement,
the added term does not vary any prediction of standard quantum mechanics.
Moreover, possible conflicts with both principles of causality and
separability can be avoided by postulating that
$\gamma\simeq\hbar^{-1}\left\langle\psi\right|V^{{\dagger}}V\left|\psi\right\rangle^{1/2}$,
where $V$ is the coupling term in the Hamiltonian giving rise to the
interaction between subsystems [$\gamma$ is the disentanglement rate in Eq.
(4)]. This postulate implies that the added nonlinear term is active only when
subsystems interact, and that time evolution is governed by the standard
Schrödinger equation when subsystems are remote (i.e. decoupled). Note that
for the examples shown in Fig. 1, the calculations are performed for the case
$\gamma=\omega_{\mathrm{d}}$. This demonstrates that a disentanglement rate
$\gamma$ having the order of
$\hbar^{-1}\left\langle\psi\right|V^{{\dagger}}V\left|\psi\right\rangle^{1/2}$
is sufficiently large to allow full suppression of entanglement.

Summary \- Further theoretical study is needed to check whether quantum
mechanics can be self-consistently reformulated based on the proposed modified
Schrödinger equation (4). We find that conflict with some well-established
physical principles, as well as many experimental observations, can be avoided
by postulating that
$\gamma\simeq\hbar^{-1}\left\langle\psi\right|V^{{\dagger}}V\left|\psi\right\rangle^{1/2}$.

The expression given by Eq. (2) for the operator $\mathcal{Q}$ is applicable
for the bipartite case, for which the entire system is divided into two
subsystems. The multipartite case, however, for which the entire system is
divided into more than two subsystems, requires a generalization of Eq. (2).
Such generalization is discussed in Ref. Buks_2306_05853 . The generalization
of the above discussed postulate (regarding the disentanglement rate $\gamma$)
for the multipartite case states that disentanglement between two given
subsystems is active only during the time when they interact.

Further insight can be gained from experimental study of entanglement in the
region where environmental decoherence is negligible Buks_014421 . Upper
bounds imposed upon the disentanglement rate $\gamma$ in Eq. (4) can be
derived from lifetime measurements of entangled states. Experimental
observations of deviation from the Born rule may provide supporting evidence
for nonlinearity (see Fig. 3).

Acknowledgments \- We thank Jakub Rembielinski, Pawel Caban, Joakim Bergli and
Klaus Molmer for useful discussions. This work was supported by the Israeli
science foundation, the Israeli ministry of science, and by the Technion
security research foundation.

##  Appendix A The Schmidt decomposition

The system’s normalized pure state vector $\left|\psi\right\rangle$ is given
by
$\left|\psi\right\rangle=\mathcal{K}_{1}C\otimes\mathcal{K}_{2}^{\mathrm{T}}$
[see Eq. (1) in the main text]. Consider the unitary transformations (the
letter $k$ is used to label the states of the original basis, whereas the
transformed states are labeled by the letter $l$)

$\displaystyle\mathcal{K}_{1}^{\mathrm{T}}$
$\displaystyle=u_{1}\mathcal{L}_{1}^{\mathrm{T}}=u_{1}\left(\left|l_{1}\right\rangle_{1},\left|l_{2}\right\rangle_{1},\cdots,\left|l_{N_{1}}\right\rangle_{1}\right)^{\mathrm{T}}\;,$
(13) $\displaystyle\mathcal{K}_{2}^{\mathrm{T}}$
$\displaystyle=u_{2}\mathcal{L}_{2}^{\mathrm{T}}=u_{2}\left(\left|l_{1}\right\rangle_{2},\left|l_{2}\right\rangle_{2},\cdots,\left|l_{N_{2}}\right\rangle_{2}\right)^{\mathrm{T}}\;,$
(14)

where $u_{1}$ ($u_{2}$) is a $N_{1}\times N_{1}$ ($N_{2}\times N_{2}$) unitary
matrix (i.e. $u_{1}^{{\dagger}}u_{1}=1$ and $u_{2}^{{\dagger}}u_{2}=1$). The
state vector $\left|\psi\right\rangle$ in the transformed basis is expressed
as

$\displaystyle\left|\psi\right\rangle$
$\displaystyle=\mathcal{L}_{1}\hat{C}\otimes\mathcal{L}_{2}^{\mathrm{T}}$
$\displaystyle={\displaystyle\sum\limits_{l_{1},l_{2}}}\hat{C}_{l_{1},l_{2}}\left|l_{1}\right\rangle_{1}\otimes\left|l_{2}\right\rangle_{2}\;,$
(15)

where the transformed matrix $\hat{C}$ is given by

$\hat{C}=u_{1}^{\mathrm{T}}Cu_{2}\;,$ (16)

and the corresponding density operator
$\rho=\left|\psi\right\rangle\left\langle\psi\right|$ is expressed as

$\rho={\displaystyle\sum\limits_{l_{1}^{\prime},l_{2}^{\prime},l_{1}^{\prime\prime},l_{2}^{\prime\prime}}}\hat{C}_{l_{1}^{\prime},l_{2}^{\prime}}\hat{C}_{l_{1}^{\prime\prime},l_{2}^{\prime\prime}}^{\ast}\left|l_{1}^{\prime},l_{2}^{\prime}\right\rangle\left\langle
l_{1}^{\prime\prime},l_{2}^{\prime\prime}\right|\;.$ (17)

The following holds

$\displaystyle\operatorname{Tr}\rho$
$\displaystyle={\displaystyle\sum\limits_{l_{1},l_{2}}}\left|\hat{C}_{l_{1},l_{2}}\right|^{2}$
$\displaystyle=\operatorname{Tr}S_{1}=\operatorname{Tr}S_{2}=\operatorname{Tr}\left(CC^{{\dagger}}\right)=\operatorname{Tr}\left(C^{{\dagger}}C\right)\;,$
(18)

where the $N_{1}\times N_{1}$ ($N_{2}\times N_{2}$) matrix $S_{1}$ ($S_{2}$)
is given by (recall that $u_{1}^{{\dagger}}u_{1}=1$ and
$u_{2}^{{\dagger}}u_{2}=1$)

$\displaystyle S_{1}$
$\displaystyle=\hat{C}\hat{C}^{{\dagger}}=u_{1}^{\mathrm{T}}Cu_{2}u_{2}^{{\dagger}}C^{{\dagger}}u_{1}^{\mathrm{T{\dagger}}}=u_{1}^{\mathrm{T}}CC^{{\dagger}}u_{1}^{\mathrm{T{\dagger}}}\;,$
(19) $\displaystyle S_{2}$
$\displaystyle=\hat{C}^{{\dagger}}\hat{C}=u_{2}^{{\dagger}}C^{{\dagger}}u_{1}^{\mathrm{T{\dagger}}}u_{1}^{\mathrm{T}}Cu_{2}=u_{2}^{{\dagger}}C^{{\dagger}}Cu_{2}\;,$
(20)

hence $\operatorname{Tr}\rho=1$ provided that $\left|\psi\right\rangle$ is
normalized. The matrix $S_{1}$ ($S_{2}$) is Hermitian and positive definite,
hence the unitary matrix $u_{1}$ ($u_{2}$) can be chosen to diagonalize
$S_{1}$ ($S_{2}$), and the eigenvalues, which are denoted by $q_{l}$, are non-
negative. For this transformation, which is called the Schmidt decomposition,
the transformed matrix $\hat{C}$ has a diagonal form

$\hat{C}_{l_{1},l_{2}}=q_{l_{1}}\delta_{l_{1},l_{2}}\;.$ (21)

The purity $P_{1}$ ($P_{2}$) is defined by
$P_{1}=\operatorname{Tr}\rho_{1}^{2}$ ($P_{2}=\operatorname{Tr}\rho_{2}^{2}$),
where $\rho_{1}=\operatorname{Tr}_{2}\rho$
($\rho_{2}=\operatorname{Tr}_{1}\rho$) is the reduced density operator of the
first (second) subsystem. With the help of the Schmidt decomposition (21), one
finds that $P_{1}=P_{2}\equiv P$, where

$\displaystyle P$ $\displaystyle={\displaystyle\sum\limits_{l}}q_{l}^{4}$
$\displaystyle=\operatorname{Tr}S_{1}^{2}=\operatorname{Tr}\left(CC^{{\dagger}}\right)^{2}=\operatorname{Tr}S_{2}^{2}=\operatorname{Tr}\left(C^{{\dagger}}C\right)^{2}\;.$

Note that $P=1\ $for a product state, and $P$ obtains its minimum value of
$1/\min\left(N_{1},N_{2}\right)$ for a maximally entangled state. The purity
$P$ is independent on the local transformations $u_{1}$ and $u_{2}$, hence it
is a constant when the subsystems are decoupled (i.e. when the interaction
between the subsystems vanishes). Using the relations

$\operatorname{Tr}\left(C^{{\dagger}}C\right)={\displaystyle\sum\limits_{k_{1}^{\prime}=1}^{N_{1}}}{\displaystyle\sum\limits_{k_{2}^{\prime}=1}^{N_{2}}}C_{{}_{k_{1}^{\prime},k_{2}^{\prime}}}^{\ast}C_{{}_{k_{1}^{\prime},k_{2}^{\prime}}}\;,$
(23)

and

$\displaystyle\operatorname{Tr}\left(C^{{\dagger}}C\right)^{2}$
$\displaystyle={\displaystyle\sum\limits_{k_{1}^{\prime},k_{1}^{\prime\prime}=1}^{N_{1}}}{\displaystyle\sum\limits_{k_{2}^{\prime},k_{2}^{\prime\prime}=1}^{N_{2}}}C_{k_{1}^{\prime},k_{2}^{\prime}}^{\ast}C_{k_{1}^{\prime},k_{2}^{\prime\prime}}C_{k_{1}^{\prime\prime},k_{2}^{\prime\prime}}^{\ast}C_{k_{1}^{\prime\prime},k_{2}^{\prime}}\;,$

one finds that the level of entanglement $1-P$ is given by

$\displaystyle 1-P$
$\displaystyle=\left(\operatorname{Tr}\left(C^{{\dagger}}C\right)\right)^{2}-\operatorname{Tr}\left(C^{{\dagger}}C\right)^{2}$
$\displaystyle=\frac{1}{2}{\displaystyle\sum\limits_{k_{1}^{\prime},k_{1}^{\prime\prime}=1}^{N_{1}}}{\displaystyle\sum\limits_{k_{2}^{\prime},k_{2}^{\prime\prime}=1}^{N_{2}}}\left|\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|^{2}\;,$
(25)

where

$\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}=C_{{}_{k_{1}^{\prime},k_{2}^{\prime}}}C_{{}_{k_{1}^{\prime\prime},k_{2}^{\prime\prime}}}-C_{k_{1}^{\prime},k_{2}^{\prime\prime}}C_{k_{1}^{\prime\prime},k_{2}^{\prime}}\;.$
(26)

Note that the term
$\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}$
vanishes unless $k_{1}^{\prime}\neq k_{1}^{\prime\prime}$ and
$k_{2}^{\prime}\neq k_{2}^{\prime\prime}$, and the following holds
$\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}=\phi_{k_{1}^{\prime\prime},k_{1}^{\prime},k_{2}^{\prime\prime},k_{2}^{\prime}}$,
thus Eq. (25) can be rewritten as

$1-P=2{\displaystyle\sum\limits_{k_{1}^{\prime}<k_{1}^{\prime\prime}}}{\displaystyle\sum\limits_{k_{2}^{\prime}<k_{2}^{\prime\prime}}}\left|\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}\right|^{2}\;.$
(27)

Note that for any product state
$\phi_{k_{1}^{\prime},k_{1}^{\prime\prime},k_{2}^{\prime},k_{2}^{\prime\prime}}=0$
[see Eq. (26)]. The above result (27) implies that
$P=1-\left\langle\mathcal{Q}\right\rangle$, where the operator $\mathcal{Q}$
is given by Eq. (2) in the main text.

## References

  * (1) E. Schrodinger,  “Die gegenwartige situation in der quantenmechanik”,  Naturwissenschaften, vol. 23, pp. 807, 1935.
  * (2) Roger Penrose,  “Uncertainty in quantum mechanics: faith or fantasy?”,  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 369, no. 1956, pp. 4864–4890, 2011.
  * (3) A. J. Leggett,  “Experimental approaches to the quantum measurement paradox”,  Found. Phys., vol. 18, pp. 939–952, 1988.
  * (4) A. J. Leggett,  “Realism and the physical world”,  Rep. Prog. Phys., vol. 71, pp. 022001, 2008.
  * (5) Eyal Buks,  “Disentanglement and a nonlinear schrödinger equation”,  Journal of Physics A: Mathematical and Theoretical, vol. 55, no. 35, pp. 355303, 2022.
  * (6) Eyal Buks,  “Thermalization and disentanglement with a nonlinear schrödinger equation”,  in J. Phys. A: Math. Theor., 2023, vol. 56, p. 025302.
  * (7) Michael R Geller,  “Fast quantum state discrimination with nonlinear positive trace-preserving channels”,  Advanced Quantum Technologies, p. 2200156, 2023.
  * (8) Steven Weinberg,  “Testing quantum mechanics”,  Annals of Physics, vol. 194, no. 2, pp. 336–386, 1989.
  * (9) Steven Weinberg,  “Precision tests of quantum mechanics”,  in THE OSKAR KLEIN MEMORIAL LECTURES 1988–1999, pp. 61–68. World Scientific, 2014.
  * (10) H-D Doebner and Gerald A Goldin,  “On a general nonlinear schrödinger equation admitting diffusion currents”,  Physics Letters A, vol. 162, no. 5, pp. 397–401, 1992.
  * (11) H-D Doebner and Gerald A Goldin,  “Introducing nonlinear gauge transformations in a family of nonlinear schrödinger equations”,  Physical Review A, vol. 54, no. 5, pp. 3764, 1996.
  * (12) Nicolas Gisin and Ian C Percival,  “The quantum-state diffusion model applied to open systems”,  Journal of Physics A: Mathematical and General, vol. 25, no. 21, pp. 5677, 1992.
  * (13) David E Kaplan and Surjeet Rajendran,  “Causal framework for nonlinear quantum mechanics”,  Physical Review D, vol. 105, no. 5, pp. 055002, 2022.
  * (14) Manuel H Muñoz-Arias, Pablo M Poggi, Poul S Jessen, and Ivan H Deutsch,  “Simulating nonlinear dynamics of collective spins via quantum measurement and feedback”,  Physical review letters, vol. 124, no. 11, pp. 110503, 2020.
  * (15) Angelo Bassi, Kinjalk Lochan, Seema Satin, Tejinder P Singh, and Hendrik Ulbricht,  “Models of wave-function collapse, underlying theories, and experimental tests”,  Reviews of Modern Physics, vol. 85, no. 2, pp. 471, 2013.
  * (16) Philip Pearle,  “Reduction of the state vector by a nonlinear schrödinger equation”,  Physical Review D, vol. 13, no. 4, pp. 857, 1976.
  * (17) Gian Carlo Ghirardi, Alberto Rimini, and Tullio Weber,  “Unified dynamics for microscopic and macroscopic systems”,  Physical review D, vol. 34, no. 2, pp. 470, 1986.
  * (18) Angelo Bassi and GianCarlo Ghirardi,  “Dynamical reduction models”,  Physics Reports, vol. 379, no. 5-6, pp. 257–426, 2003.
  * (19) IJ Arnquist, FT Avignone III, AS Barabash, CJ Barton, KH Bhimani, E Blalock, B Bos, M Busch, M Buuck, TS Caldwell, et al.,  “Search for spontaneous radiation from wave function collapse in the majorana demonstrator”,  Physical Review Letters, vol. 129, no. 8, pp. 080401, 2022.
  * (20) Eyal Buks,  Quantum mechanics - Lecture Notes,  http://buks.net.technion.ac.il/teaching/, 2023.
  * (21) William K Wootters,  “Entanglement of formation of an arbitrary state of two qubits”,  Physical Review Letters, vol. 80, no. 10, pp. 2245, 1998.
  * (22) JM Radcliffe,  “Some properties of coherent spin states”,  Journal of Physics A: General Physics, vol. 4, no. 3, pp. 313, 1971\.
  * (23) Charles H Bennett, Debbie Leung, Graeme Smith, and John A Smolin,  “Can closed timelike curves or nonlinear quantum mechanics improve quantum state discrimination or help solve hard problems?”,  Physical review letters, vol. 103, no. 17, pp. 170502, 2009.
  * (24) Krzysztof Kowalski,  “Linear and integrable nonlinear evolution of the qutrit”,  Quantum Information Processing, vol. 19, no. 5, pp. 1–31, 2020\.
  * (25) Bernd Fernengel and Barbara Drossel,  “Bifurcations and chaos in nonlinear lindblad equations”,  Journal of Physics A: Mathematical and Theoretical, vol. 53, no. 38, pp. 385701, 2020.
  * (26) K Kowalski and J Rembieliński,  “Integrable nonlinear evolution of the qubit”,  Annals of Physics, vol. 411, pp. 167955, 2019.
  * (27) Angelo Bassi and Kasra Hejazi,  “No-faster-than-light-signaling implies linear evolution. a re-derivation”,  European Journal of Physics, vol. 36, no. 5, pp. 055027, 2015.
  * (28) Thomas F. Jordan,  “Assumptions that imply quantum dynamics is linear”,  Phys. Rev. A, vol. 73, pp. 022101, Feb 2006.
  * (29) Joseph Polchinski,  “Weinberg?s nonlinear quantum mechanics and the einstein-podolsky-rosen paradox”,  Physical Review Letters, vol. 66, no. 4, pp. 397, 1991.
  * (30) Bassam Helou and Yanbei Chen,  “Extensions of born?s rule to non-linear quantum mechanics, some of which do not imply superluminal communication”,  in Journal of Physics: Conference Series. IOP Publishing, 2017, vol. 880, p. 012021.
  * (31) Jakub Rembieliński and Paweł Caban,  “Nonlinear evolution and signaling”,  Physical Review Research, vol. 2, no. 1, pp. 012027, 2020.
  * (32) Jakub Rembieliński and Paweł Caban,  “Nonlinear extension of the quantum dynamical semigroup”,  Quantum, vol. 5, pp. 420, 2021.
  * (33) Zakarias Laberg Hejlesen,  “Nonlinear quantum mechanics”,  Master’s thesis, 2019.
  * (34) Thomas F Jordan,  “Why quantum dynamics is linear”,  in Journal of Physics: Conference Series. IOP Publishing, 2009, vol. 196, p. 012010.
  * (35) Eyal Buks,  “Partial disentanglement in a multipartite system”,  arXiv:2306.05853, 2023.
  * (36) Eyal Buks and Banoj Kumar Nayak,  “Quantum measurement with recycled photons”,  Physical Review B, vol. 105, no. 1, pp. 014421, 2022.


# Explainable expected goal models for performance analysis in football
analytics

Mustafa Cavus
Warsaw University of Technology
Faculty of Mathematics and Information Science
Warsaw, Poland
Eskisehir Technical University
Department of Statistics
Eskisehir, Turkey
<EMAIL_ADDRESS>

Przemysław Biecek
Warsaw University of Technology
Faculty of Mathematics and Information Science
Warsaw, Poland
<EMAIL_ADDRESS>
The work on this paper is financially supported by the NCN Sonata Bis-9 grant
2019/34/E/ST6/00052

###### Abstract

The expected goal provides a more representative measure of the team and
player performance which also suit the low-scoring nature of football instead
of score in modern football. The score of a match involves randomness and
often may not represent the performance of the teams and players, therefore it
has been popular to use the alternative statistics in recent years such as
shots on target, ball possessions, and drills. To measure the probability of a
shot being a goal by the expected goal, several features are used to train an
expected goal model which is based on the event and tracking football data.
The selection of these features, the size and date of the data, and the model
which are used as the parameters that may affect the performance of the model.
Using black-box machine learning models for increasing the predictive
performance of the model decreases its interpretability that causes the loss
of information that can be gathered from the model. This paper proposes an
accurate expected goal model trained consisting of 315,430 shots from seven
seasons between 2014-15 and 2020-21 of the top-five European football leagues.
Moreover, this model is explained by using explainable artificial intelligence
tool to obtain an explainable expected goal model for evaluating a team or
player performance. To best of our knowledge, this is the first paper that
demonstrates a practical application of an explainable artificial intelligence
tool aggregated profiles to explain a group of observations on an accurate
expected goal model for monitoring the team and player performance. Moreover,
these methods can be generalized to other sports branches.

_Keywords_ Football $\cdot$ Expected goal $\cdot$ Machine learning $\cdot$
Explainable artificial intelligence $\cdot$ Aggregated profiles

##  1 Introduction

In recent years, the exponential speed of improvement in the technologies
supporting the collection, storage, and analysis of data has a revolutionary
effect on football analytics as well as many fields. The easy accessibility of
data provides a great potential to propose several key performance metrics
measuring several aspects of the play such as pass evaluation, quantifying
controlled space, evaluating shots, and goal-scoring opportunities through
possession values. One of these prominent metrics is _expected goal_ (xG)
which is the most notable one in football talkshows in TV and end-of-match
statistics nowadays. It is proposed by Green [1] to quantify the probability
of a shot being the goal. The reason behind the development of such a metric
is to propose a metric to represent the low-scoring nature of football rather
than the other sports. It is an ordinary story in football that a team
dominated the game -they had many scoring opportunities- but could not goal,
and the opponent won the match by converting one of the few goal opportunities
to goal they created. In this case, the xG is used as a useful indicator of
the score. It can be defined as the mean of a large number of independent
observations of a random variable which is the shots from the statistical
point of view. Besides being a good representative of the score, it is also a
good indicator which is used to predict the future team performance [2]. There
are many studies is conducted to train a machine learning model learned from
the predictors such as shot type, distance to goal, angle to goal for
predicting the value of xG [3, 4, 5, 6, 7, 8, 9, 10]. The data used to train
the xG model is usually highly unbalanced. It causes an important problem that
seen in some of these papers is the poor prediction performance of the models
on minority class [11]. In this paper, we aimed to propose an accurate xG
model in terms of both majority and minority classes.

One of the practical applications of the xG models, which is the main focus of
this paper, is performance evaluation [12, 13]. Brechot and Flepp [14]
proposed to use the xG models for performance evaluation instead of match
outcomes which may easily be influenced by randomness in short-term results.
They introduced a chart built upon the concept of the xG by plotting the
teams’ ranking in the league table against their rankings based on xG.
Moreover, they proposed some useful metrics calculated based on xG such as
offensive and defensive ratios. Fairchild et al. [15] focused on ways for
evaluating the xG model goes beyond the accuracy, which is second from a
player and team evaluation perspective on offensive and defensive efficiency
by comparing the xG metric with the actual goals. They created the xG model
for Major League Soccer in the USA and Canada. However, these papers consider
only the output of the xG model. Thanks to the XAI tools, it is possible to
explain a black-box machine learning model’s behavior at the local and global
levels. In this way, we can gather more information from the model not only
its prediction and also its behavior. The one of most commonly used tool at
local-level is the ceteris-paribus (CP) profiles that show the change of model
prediction would change for the value of a feature on a single observation
[16]. Its use to explain only one observation is a limitation, but by
aggregating these profiles, it is possible to explain more than one
observation at the same time. Actually, partial dependence profile (PDP) is
used to explain the relationship between a feature and the response [17].
However, the PDP is the estimation of the mean of the CP profiles for all
observations in a dataset, not just some of the observations. In football,
offensive performance can be measured through the shots taken by a team or
player. In this paper, we introduced a practical application of an XAI tool
based on the aggregation of the CP profiles which are used for the local-level
explanation of the model behavior. By this approach, we can evaluate a player
or team’s performance and answer the what-if type questions about the
performance.

The main contributions of this paper are: (1) proposing the most accurate xG
model, in terms of both majority and minority classes, trained on the data
consists 315,430 shots from seven seasons between 2014-15 and 2020-21 of the
top-five European football leagues, and (2) introducing a novel team/player
performance evaluation approach which is a practical application of the XAI
tools in football based on the aggregation of the CP profiles. We believe that
the approaches given in this paper can be generalized for the other branches
of sport.

The remainder of this paper is structured as follows: Sec. II introduces the
mathematical background of the xG models and xG model training. The
performance of the trained xG models are investigated in Sec. III. Lastly, in
Sec. IV, we introduce how the aggregated profiles are created and used to
evaluate the performance of a player or team. Moreover, we demonstrate a
practical application of the aggregated profiles based on the xG model for
player and team levels.

##  2 Expected goal models

Consider $\mathbf{X}\subseteq\mathbb{R}^{d}$ is a $d$-dimensional feature
vector, $Y\in\\{0,1\\}$ is the label vector of response variable. The dataset
is denoted by $D=\\{(\mathbf{x}_{i},y_{i})\\}^{n}_{i=1}$ where each sample
$(x_{i},y_{i})$ is independently sampled from the joint distribution with
density $p(\mathbf{x},y)$ which includes an instance
$\mathbf{x}_{i}\in\mathbf{X}$ and a label $y_{i}\in Y$. The goal of a binary
classifier is to train an optimal mapping function as follows:

$f:\mathbf{X}\rightarrow Y$ (1)

by minimizing a loss function is $L(f)=P[Y\neq f(\mathbf{X})]$. The xG model
is a special case of supervised classification task which has a binary outcome
that takes the values are goal or not. Here $Y$ is the target variable which
shows the goal or not of a shot, and the X are the features that are used to
predict the value of $Y$. The most commonly used features are distance to
goal, angle to goal, shot type, last action, etc. The calculation of the xG
value from the xG model can be easily algorithmized:

  1. 1.

the individual scoring probabilities of the shots are calculated,

  2. 2.

these probabilities are summed over for a player, or a team to derive the
cumulative xG value [14].

For example, assuming that a team had three shots in a match with
probabilities of 0.50, 0.20, and 0.05 means that the team has generated
chances worth 0.75 xG. It can be calculated not only for a match, but also for
different time periods such as season(s). This creates many different
practical applications’ opportunities.

In the following subsections, we describe the data we used to train xG models,
and give the steps about the pre-processing of data. Then, we introduce the
tools we used in model training and explanation. Lastly, the problem in the xG
model is imbalanced data is discussed and the solution way is mentioned.

###  2.1 Description of the Data

The issue that needs to be discussed, before the training of a xG model, is
the characteristics of the data used to train the xG models. It is expected
that the style of play changes over time and varies between leagues from the
football enthusiasts’ point of view. The answer of "How it can be determined
whether this situation has occurred or not?" can be found in Robberechts and
Davis [18]. They conducted an extensive experimental study to investigate the
frequently asked data-related questions such as "How much data is needed to
train an accurate xG model?", "Are xG models league-specific?", and "Does data
go out of date?" that may affect the performance of an xG model. Their results
show that five seasons of data are needed to train a complex xG model, the
data does not go out of date, and using league-based xG models does not
increase the accuracy significantly. We determined our model development
strategy considering these findings in this paper.

We focus in our paper on 315,430 shots-related event data (containing 33,656
goals $\sim 10.66\%$ of total shots) from the 12,655 matches in 7 seasons
between 2014-15 and 2020-21 from the top-five European football leagues which
are Serie A, Bundesliga, La Liga, English Premier League, Ligue 1. The dataset
is collected from Understat111https://understat.com by using the R-package
worldfootballR [19] and excluded the 1,012 shots resulting in own goals due to
their unrelated pattern from the concept of the model. The package provides
useful functions to gather and handle the shots data by matches, seasons, and
leagues from the various data sources such as FBref222https://fbref.com/en/,
Transfermarkt333https://www.transfermarkt.com, and
Fotmob444https://www.fotmob.com. The detailed information, such as type and
description, about the features used in the model, are given in Table 1.

Table 1: Details of the variables used to train our xG model Features |  Type |  Description
---|---|---
status |  categorical |  situation that the shot is being a goal (0: no goal, 1: goal)
minute |  continuous |  minute of shot between 1 and 90 + possible extra time
home and away |  categorical |  status of the shooting team (home or away)
situation |  categorical |  situation at the time of the event (Direct freekick, From corner, Open play, Penalty, Set play)
shot type |  categorical |  type based on the limb used by the player to shot (Head, Left foot, Right foot, Other part of the body)
last action |  categorical |  last action before the shot (Pass, Cross, Rebound, Head Pass, and 35 more levels)
distance to goal |  continuous |  distance from where the shot was taken to the goal line $([0.295,84.892]$ in meters)
angle to goal |  continuous |  angle of the throw to the goal line $([0.10^{\circ},90^{\circ}])$

The summary statistics of the shots and goals, such as the number ($\\#$) of
matches, shots, goals, the mean ($\mu$) of shots and goals per match, and the
conversion percent ($\%$) of a shot to goal, per league over seven seasons are
given in Table 2.

Table 2: The summary statistics of the shots and goals, such as the number ($\\#$) of matches, shots, goals, the mean ($\mu$) of shots and goals per match and the conversion percent ($\%$) of a shot to goal for per league over seven seasons League |  $\\#$Match |  $\\#$Shot | $\mu_{Shot}$ |  $\\#$Goal | $\mu_{Goal}$ | $\%$
---|---|---|---|---|---|---
Bundesliga | 2,141 | 55,129 | 25.7 | 6,161 | 2.88 | 11.2
EPL | 2,650 | 66,605 | 25.1 | 6,951 | 2.62 | 10.4
La Liga | 2,648 | 62,028 | 23.4 | 6,854 | 2.59 | 11.0
Ligue 1 | 2,557 | 61,053 | 23.9 | 6,438 | 2.52 | 10.5
Serie A | 2,659 | 70,615 | 26.6 | 7,252 | 2.73 | 10.3
Mean | 2,531 | 63,086 | 24.9 | 6,371 | 2.67 | 10.7
Total | 12,655 | 315,430 | - | 33,656 | - | -

According to the summary statistics, the conversion percent of all leagues is
0.107 and the Bundesliga has the highest rate is 0.112 while the Serie A has
the lowest rate of 0.103. The interesting statistics related to Serie A is the
percent of conversion to goals is the lowest, however it is the league with
the highest number of shots per match.

Figure 1: The distribution of angle to goal and distance to goal of shots
regarding goal status in the last seven seasons of top-five European football
leagues

As expected, when the angle of the shot decreases, the probability of goal
decreases, and the distance to goal increases, and the probability of goal
decreases. Fig 1 shows that: (1) the distribution of the distance to goal of
shots seems similar for each league, (2) the range of distance to goal is
between 0-25 meters, (3) the optimal angle to goal is about $30^{\circ}$. The
distribution of angle to goal and distance to goal of shots regarding goal
status in the last seven seasons of the leagues seems similar. It is observed
that the distributions of the two most important features affecting the
probability of goal between leagues are similar.

###  2.2 Pre-processing of the Data

The pre-processing steps are necessary before modeling such as the
transformation of some features. The location of shots is given in the
coordinates system in the dataset as $L_{i}\in[0,1]$ and $W_{i}\in[0,1]$ as in
Fig 2. We must calculate the distance and angle to goal of the shots which are
the two most important features in xG models based on the coordinate values
because the coordinates are not meaningful in the interpretation of the model.
Before calculating these features, we standardized a football pitch is
$L=105$m $\times$ $W=68$m in size, however some pitches may have different
dimensions in reality. The size of the pitch is that the length should be
between 90 and 120 meters and the width should be between 45 and 90 meters are
limited by the rules of The International Football Association
Board555https://digitalhub.fifa.com/m/5371a6dcc42fbb44/original/d6g1medsi8jrrd3e4imp-
pdf.pdf.

Figure 2: The standard dimension of a football pitch

The following transformation are used to calculate the distance ($X^{DTG}$)
and angle to goal ($X^{ATG}$) features:

$X^{DTG}_{i}=\sqrt{[105-(L_{i}\times 105)]^{2}+[32.5-(W_{i}\times 68)^{2}]}$
(2)

where $L$ and $W$ are the coordinates of a shot.

$X^{ATG}_{i}=\left|\frac{a_{i}}{b_{i}}\times\frac{180}{\pi}\right|$ (3)

where $a_{i}=\arctan[7.32\times[105-(L_{i}\times 105)]]$ and
$b_{i}=[105-(L_{i}\times 105)]^{2}+[32.5-(W_{i}\times 68)]^{2}-(7.32/2)^{2}$.
We used $X^{DTG}$ and $X^{ATG}$ features in model training instead of the
original coordinates $L_{i}$ and $W_{i}$ given in the raw data.

###  2.3 Model Training

We use the forester [20] AutoML tool to train various tree-based
classification models from XGBoost [21], randomForest [22], LightGBM [23], and
CatBoost [24] libraries. These models do not provide any pre-processing steps
like missing data imputation, encoding, or transformation and show quite good
performance in the presence of outlier(s) in the dataset which is used to
train models. We use the train-test split (80-20) to train and validate the
models. Moreover, another advantage of the forester is that provides an easy
connection to DALEX [25] model explanation and exploration ecosystem.

###  2.4 Data balancing

Imbalancedness is a kind of problem when one of the classes of the target
variable is rare over the sample in the classification task. In this case, the
model learns more from the majority class that which lead to poor
classification performance of the minority class. The imbalance problem
between the class of target feature is considered a separate field called
imbalanced learning in machine learning. In imbalanced learning, there are
three main strategies to train the models [26]: (1) balancing the dataset, by
using over or under-sampling methods (2) using cost-sensitive learners, and
(3) using ensemble learning models.

The target feature we used in the model is imbalanced (90$\%$-10$\%$), and to
handle this problem we prefer to use a balancing strategy which is the random
over-sampling method provided by the ROSE package in R [27]. It consists of a
smoothed bootstrap-based technique which is proposed by Menardi and Torelli
[28].

##  3 Exploration and explanation of the proposed xG model

In this section, the performance of the trained xG models are investigated in
terms of several metrics under different sampling strategies. Then, the best
performing xG model is compared with the alternatives in the literature. The
aggregated profiles that used to analysis model behavior of the proposed xG
model is introduced for evaluating the performance at player and team levels
in the last part. Moreover, the practical applications of aggregated profiles
are given on the xG model.

###  3.1 Model Performance

For the model level exploration, the first step is usually related to model
performance. Different measures may be used such as precision, recall, F1,
accuracy, and AUC. However, these measures do not measure the performance of a
classification model in case of imbalanced data. The Mathews correlation
coefficient, brier score, log-loss and balanced accuracy are used to measure
the classification ability of both of the classes in this task. Unfortunately,
there are only a few papers have used appropriate measures regarding this
problem in the literature related to the xG model. Considering this situation,
we reported the results in terms of both of these two groups of measures.
During the study of this paper, we restricted the analysis to a comparison of
model performance between train and test data. The performance measures of the
random forest, catboost, lightgbm, and xgboost models trained on the train,
over-sampled train, and under-sampled train set are calculated and given in
Table 3.

Table 3: Performance of trained xG models Model |  Sampling |  Recall |  Precision |  F1 |  Accuracy |  AUC |  MCC |  Brier Score |  Log-loss |  Balanced Accuracy
---|---|---|---|---|---|---|---|---|---|---
random forest |  over |  0.955 |  0.921 |  0.938 |  0.937 |  0.985 |  0.874 |  0.073 |  0.275 |  0.937
under |  0.860 |  0.882 |  0.871 |  0.872 |  0.954 |  0.745 |  0.106 |  0.357 |  0.872
original |  0.297 |  0.891 |  0.446 |  0.921 |  0.975 |  0.488 |  0.051 |  0.175 |  0.646
catboost |  over |  0.739 |  0.756 |  0.747 |  0.751 |  0.835 |  0.501 |  0.167 |  0.501 |  0.750
under |  0.722 |  0.753 |  0.737 |  0.741 |  0.823 |  0.483 |  0.171 |  0.512 |  0.513
original |  0.187 |  0.718 |  0.296 |  0.905 |  0.820 |  0.334 |  0.075 |  0.264 |  0.589
xgboost |  over |  0.710 |  0.751 |  0.730 |  0.737 |  0.816 |  0.476 |  0.175 |  0.523 |  0.737
under |  0.717 |  0.751 |  0.734 |  0.738 |  0.818 |  0.477 |  0.174 |  0.521 |  0.739
original |  0.175 |  0.709 |  0.281 |  0.904 |  0.816 |  0.321 |  0.076 |  0.266 |  0.583
lightgbm |  over |  0.713 |  0.748 |  0.730 |  0.736 |  0.814 |  0.473 |  0.176 |  0.526 |  0.736
under |  0.714 |  0.748 |  0.730 |  0.736 |  0.814 |  0.471 |  0.176 |  0.525 |  0.736
original |  0.179 |  0.689 |  0.284 |  0.904 |  0.813 |  0.318 |  0.077 |  0.267 |  0.585
*The best value is given in bold for each metric.

###  3.2 Comparison of Model Performance

We compared the performance of our proposed models with the models in the
literature [10], [5], [4], [11], [29], [30], [31] in terms of precision,
recall, accuracy, F1, AUC, log-loss, Brier score, and mean absolute error
(MAE) in Table 4. We decided to use these measures because the authors of
these papers reported the performance of the models. The reason for the empty
cells in the table is that these values for the relevant models have not been
reported in these papers.

Table 4: Performance Comparison of the proposed xG model with the models in the literature Paper |  Model |  Precision |  Recall |  F1 |  AUC |  Log-loss |  Brier score |  MAE
---|---|---|---|---|---|---|---|---
Eggels et al. (2016) |  random forest |  0.785 |  0.822 |  0.800 |  0.814 |  - |  - |  -
decision tree |  0.698 |  0.678 |  0.676 |  0.677 |  - |  - |  -
logistic regression |  0.715 |  0.650 |  0.673 |  0.697 |  - |  - |  -
ada-boost |  0.624 |  0.773 |  0.688 |  0.670 |  - |  - |  -
Pardo (2020) |  logistic regression |  - |  - |  - |  - |  0.261 |  - |  -
xgboost |  - |  - |  - |  - |  0.257 |  - |  -
neural network |  - |  - |  - |  - |  0.260 |  - |  -
Tippana (2020) |  Poisson regression |  - |  - |  - |  - |  - |  - |  6.5
Anzer and Bauer (2021) |  gbm |  0.646 |  0.181 |  - |  0.822 |  - |  - |  -
logistic regression |  0.611 |  0.108 |  - |  0.807 |  - |  - |  -
ada-boost |  0.548 |  0.201 |  - |  0.816 |  - |  - |  -
random forest |  0.611 |  0.163 |  - |  0.794 |  - |  - |  -
Haaren (2021) |  boosting machine |  - |  - |  - |  0.793 |  - |  0.082 |  -
Umami et al. (2021) |  logistic regression |  - |  0.967 |  - |  - |  - |  - |  -
Fernandez et al. (2021) |  xgboost |  - |  - |  - |  - |  0.254 |  - |  -
Our models |  random forest |  0.921 |  0.955 |  0.938 |  0.985 |  0.275 |  0.073 |  2.0
catboost |  0.756 |  0.739 |  0.747 |  0.835 |  0.501 |  0.167 |  3.0
xgboost |  0.751 |  0.710 |  0.730 |  0.816 |  0.523 |  0.175 |  3.1
lightgbm |  0.748 |  0.713 |  0.730 |  0.814 |  0.526 |  0.176 |  3.1
*The best value is given in bold for each metric.

Eggels et al. [10], Pardo [5], and Anzer and Bauer [11] trained several xG
models and reported their performances. The random forest, xgboost, and gbm
models outperform others, respectively in these papers. It is seen that our
proposed random forests model outperforms the others in terms of precision,
F1, AUC, Brier score, and MAE. The model proposed by Fernandez et al. [31] has
a lower log-loss value than the model that we proposed. However, no
information is reported about the performance of their model in terms of
recall and precision, so it is difficult to be certain of its performance. It
is same for the model proposed by Umami et al. [30].

###  3.3 Model Behavior: The Aggregated Profiles

The XAI tools are classified under two main sections are local and global
levels. The local-level explanations are used to explain the behavior of a
black-box model for a single observation while the global-level explanations
are used for an entire dataset. However, our need is to explain the model for
a group of observations, i.e. for a player or team. That’s why we introduce
the aggregated profiles (AP) which can be used for a group of observations.

The idea behind AP is the aggregation of the CP profiles that show how the
change of a model’s prediction would change regarding the value of a feature.
In other words, the CP profile is a function that describes the dependence of
the conditional expected value of the response on the value z of the $j^{th}$
feature [16]. The AP can be defined simply as the averaging of the CP profiles
which are considered. The value of an AP for model $f(.)$ and feature $X_{j}$
at $z$ is defined as follows:

$g_{AP}^{j}(z)=E_{\textbf{X}}^{-j}[f(\textbf{X}^{j|z})]$ (4)

where $g_{AP}$ is the expected value of the model predictions when $X_{j}$ is
fixed at $z$ over the marginal distribution of $\textbf{X}_{j|z}$. The
distribution of $\textbf{X}_{j|z}$ can be estimated by using the mean of CP
profiles for $X_{j}$ as an estimator of the AP:

$\hat{g}^{j}_{AP}(z)=\frac{1}{k}\sum_{i=1}^{k}f(\textbf{x}^{ij|z})$ (5)

where $k$ is the number of profiles that are aggregated. The difference
between the AP and PDP is the number of aggregated profiles. The PDP is the
aggregation of all profiles which are calculated on the entire dataset while
the AP is the aggregation of several profiles. The innovative part of this
paper is to propose a performance evaluation by aggregating CP profiles of
shots (i.e. observations) taken by a team or player during the period of
interest. This evaluation can be defined as post-game analysis if the period
is a game or games. Also, it can be used after the first half of a match to
decide the second-half strategy. In this way, answers to what-if questions can
be created based on teams or players, not based on shots that are not meant to
evaluate.

###  3.4 Practical Applications of the AP in Football

In this part, we demonstrate the practical applications of the aggregated
profiles, which is constructed on the proposed xG model, for team and player
levels in football. Firstly, consider the match of Schalke 04 vs. Bayern
Munich which is played in Bundesliga on Jan 24, 2021. The end-of-match
statistics such as number of goals ($\\#$Goal), expected goal (xG), number of
shots ($\\#$Shot), mean angle to goal ($\mu_{ATG}$), and mean distance to goal
($\mu_{DTG}$) for each team over the match are given in Table 5.

Table 5: The End-of-match statistics of Schalke 04 vs. Bayern Munich in the match is played on Jan 24, 2021 Team |  $\\#$Goal | xG |  $\\#$Shot | $\mu_{ATG}$ | $\mu_{DTG}$
---|---|---|---|---|---
Schalke 04 | 0 | 2.58 | 13 | 25.18 | 18.03
Bayern Munich | 4 | 9.92 | 31 | 27.56 | 16.99

According to the match statistics, Bayern Munich took 31 shots while Schalke
04 took 13 shots in this match, 4 of them being goals and Bayern won the game
0-4. The xG values of the teams show the created goal chances over the shots.
It is expected that the final score of the match is 3-10 in terms of expected
goal. However, Schalke 04 could not find any goal while Bayern Munich found
four goals. The offensive efficiency (actual goals - expected goal) of both
teams are not good, because they only found four goals considering 12.5 xG.
The observed mean angle to goal is about $25^{\circ}$ for Schalke 04 and
$27.56^{\circ}$ for Bayern Munich, and the observed mean distance to goal is
18 meters for Schalke 04 and 17 meters for Bayern Munich.

The AP of the Schalke 04 and Bayern Munich for angle and distance to goal are
given in Fig 3. In the figure of AP, X-axis represents the value of the
interested feature, and y-axis represents the average prediction of the xG for
per-shot.

Figure 3: The aggregated xG profiles of Schalke 04 and Bayern Munich for angle
and distance to goal in the match is played on Jan 24, 2021

As known, if the shots are taken at a steeper angle and closer to the goal,
the value of xG increases. It is seen that the average of the xG is about
constant after 30 meters according to the AP for distance to goal. Evaluating
the average distances to goal, Schalke 04 is farther from the goal than
Bayern. This may be one of the reasons why they did not find any goal. They
could have increased the average xG for pre-shot by around 40 percent if they
reduced the average distance to goal from 18 meters to 15 meters. Moreover, if
Schalke had taken the shots with an average angle of $35^{\circ}$ instead of
an average angle of $25^{\circ}$, the average value of xG for per-shot would
have increased by $20\%$. This potential change can be considered to improve
the team’s performance for the next match(s). As seen, AP provides to evaluate
team performance after a match and determine how the team can increase the
value of xG with possible improvements during the game.

Secondly, consider the player performance of Burak Yilmaz is the striker of
Lille OSC from Ligue 1, Lionel Messi is the midfielder of FC Barcelona from La
Liga, and Robert Lewandowski is the striker of Bayern Munich from Bundesliga
in the season of 2020-21. According to the end-of-season statistics given in
Table 6, they took 66, 195, 132 shots, and 16, 30, 40 of them being goals
during the season, respectively. The player with the highest ability to
convert shots into goals is Robert Lewandowski with $30\%$, Burak Yilmaz is
second with $25\%$ and Lionel Messi is third with $15\%$. The ability isn’t
just about skill, it may also related to how players defend against and the
average angle and distance from which they use the shots. Robert Lewandowski
shot from the shorter distance and more right angle on average, while the
other players shot from the relatively long distance and narrow angle. It can
be roughly said that the percent of the players converting shots into goals,
in other words the created xG values, are correlated with distance and angle
to goal.

Table 6: The End-of-season statistics of Burak Yilmaz (BY), Lionel Messi (LM), and Robert Lewandowski (RL) in the season of 2020-21 Player |  $\\#$Game |  $\\#$Goal | xG |  $\\#$Shots | $\mu_{ATG}$ | $\mu_{DTG}$
---|---|---|---|---|---|---
BY | 24 | 16 | 24.48 | 66 | 21.72 | 19.67
LM | 35 | 30 | 69.83 | 195 | 21.33 | 19.39
RL | 28 | 40 | 62.52 | 132 | 33.97 | 13.07

The AP of Burak Yilmaz, Lionel Messi, and Robert Lewandowski for angle and
distance to goal in the season of 2020-21 are given in Fig 4.

Figure 4: The aggregated profiles of Burak Yilmaz, Lionel Messi, and Robert
Lewandowski for angle and distance to goal in the season of 2020-21

The average value of the xG is about similar for each player for distance to
goal. Only the average value of Robert Lewandowski is slightly higher than the
others after 10 meters. It means that if the players try to take shots at 15
meters from the goal instead of the observed mean distances seen in Table 6,
the average expected goal of the player per shot may increase about $20\%$.
When the AP is examined in terms of the feature is angle to goal, the average
xG of the player is about same from the start point to $25^{\circ}$, the AP of
Robert Lewandowski is getting increase after this angle. It can be said that
if the players try to take shots at $50^{\circ}$, the highest average xG of
Robert Lewandowski will be about 0.5. At the player level as well as at the
team level, the AP provides pretty practical information about performance
evaluation and provisioning, and it is also very useful for comparing the
players that play at similar position.

##  4 Conclusion

The papers to date aim that to evaluate a team or player’s performance that
only consider the output of the xG model. Comparing the actual goals and
expected goals, they provide some statistics based on these difference for
evaluating the defensive and offensive performance of a team or a player.
However, we focused to use the xG model’s behavior to observe the relationship
between features and response which is the xG value. In this way, we can
suggest to how a team’s performance can be improved by changing the strategies
based on the features that effect the xG value such as distance to goal, angle
to goal, and others. To do this, we first proposed an accurate xG model which
is trained a random forests model on the data consist seven seasons of top-
five European leagues. This model predicts both possible outcomes of the
output of the xG model, goal and no goal, better than the alternatives in the
literature. To obtain this model, we balanced the data by using the random
over-sampling method to solve the imbalance problem that is often ignored in
similar papers. Thus, the model we proposed learned quite well from both of
the classes. The interesting thing in this process is that we detected some
changes in the behavior of the model trained on the data obtained after over-
sampling through PDP curves. We have included a detailed discussion on this
problem in Discussion section.

We evaluated the performance at the level of the team and the player by using
the accurate xG model we proposed, in the practical application part. The AP
of the model on the team level show to relationship between the interested
features and the target variable is the average of xG values predicted over
the observations. The usage of AP is practical to see the effect of the
changes on the values of interested variable on the xG value, and also
comparison of the players who are play in similar positions. For example, it
can be extracted from Fig 4 that the average xG of Robert Lewandowski’s shots
may be higher than other players if he shoots from a steeper angle. As seen,
detailed information about performance evaluation can be obtained with the
practical use of AP, which is one of the XAI tools, on the xG models.

##  5 Discussion

Domain-specific applications of XAI tools enable key insights to be extracted
from a black-box model. In this context, this paper focuses the practical
application of the AP for explaining more than one observation, not an
observation or entire dataset. These tools may be referred to as semi-global
explainers for easier understanding in the XAI domain. It is seen from the
examples discussed: AP can be used to extract provisions for performance
analysis of a team or a player from xG models, which have been frequently used
in football in recent years. In addition, comparisons of similar players and
teams can be made in terms of interested features. Since these introduced
methods can be generalized to other sports branches such as ice hockey which
is used the xG models, its widespread effect will not be limited to football.

Another discussion we want to mention in this paper is the effect of balancing
methods on the model’s behavior. It is known that these methods for imbalanced
datasets in binary classification tasks is a very commonly used solution to
improve the prediction performance of ML models in both the classes. However,
how balancing the observation of classes causes effects on the model behavior
is a subject that has not been discussed yet. The only paper in the literature
about this is Patil et al. [32] discussed that whether the model is reliable
or not after balancing. They decided to verify the reliability of the models
and oversampling method with the help of the feature importance which is one
of the XAI techniques. Their results demonstrate that the higher accuracy
obtained by the over-sampled dataset while ensuring that the oversampling does
not alter the feature correlation of the original dataset. This paper does not
satisfy to decide the change of model behavior, because it only examined the
change of order of features’ importance in the model. The point we want to
raise is the model behavior in terms of PDP’s values, because it provides more
detailed information than the feature importance to detecting the change in
the model behavior. Thus, we compared the behavior of the original, over-
sampled, and under-sampled versions of our proposed model using PDP curves for
some features in Fig 5.

Figure 5: The behavior comparison of the random forest models trained on
original, over-sampled, and under-sampled data in terms of PDP curves

It is seen that behavior of the PDP curves for the feature is distance to goal
are not same for the models. For the model trained on original data, the PDP
curve is increasing from start point to 20 meters, then slowly increasing
after some fluctations between 20 and 30 meters. However, the fluctations are
seen on different range of the feature for the models trained on over and
under-sampled data. This is a sign that there has been a change in model
behavior after using balancing methods. We would like to draw attention to the
need for careful consideration of this situation, and start a discussion on
this subject in the literature based on our findings.

## Acknowledgement

The work on this paper is financially supported by the NCN Sonata Bis-9 grant
2019/34/E/ST6/00052.

## Supplemental Materials

The R codes needed to reproduce the results can be found in the following
GitHub repository: https://github.com/mcavs/Explainable_xG_model_paper.

## References

  * [1] S. Green.  Assessing the performance of premier league goalscorer.  OptaPro Blog, 2012.
  * [2] F. Cardoso, S. González-Víllora, J. Guilherme, and I. Teoldo.  Young Soccer Players With Higher Tactical Knowledge Display Lower Cognitive Effort.  Percept Mot Skills, vol. 126, pp. 499–514, 2019.
  * [3] A. Rathke.  An examination of expected goals and shot efficiency in soccer.  J Hum Sport Exerc, vol. 12, pp. 514-–529, 2017.
  * [4] T. Tippana.  How accurately does the expected goals model reflect goalscoring and success in football?,  Bachelor’s Thesis, Aalto University, 2020.
  * [5] M. Pardo.  Creating a model for expected goals in football using qualitative player information.  Master’s thesis, Universitat Politecnica de Catalunya, 2020.
  * [6] C. Herbinet.  Predicting football results using machine learning techniques.  MEng thesis, Imperial College London, 2018.
  * [7] E. Wheatcroft and E. Sienkiewicz.  A probabilistic model for predicting shot success in football.  arXiv preprint, arXiv:2101.02104, 2021.
  * [8] L. Bransen and J. Davis.  Women’s football analyzed: interpretable expected goals models for women.  In AI for Sports Analytics (AISA) Workshop at IJCAI 2021.
  * [9] S. Sarkar and S. Kamath.  Does luck play a role in the determination of the rank positions in football leagues? A study of Europe’s big five.  Ann Oper Res, 2021.
  * [10] H. Eggels, R. Van Elk, and M. Pechenizkiy.  Explaining soccer match outcomes with goal scoring opportunities predictive analytics.  In proceedings of the Workshop on Machine Learning and Data Mining for Sports Analytics 2016 co-located with the 2016 European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases, 2016.
  * [11] G. Anzer and P. Bauer.  A goal scoring probability model for shots based on synchronized positional and event data in football (soccer).  Front. Sports Act. Living, vol. 3, pp. 1–15, 2021.
  * [12] T. Kharrat, I. G. MacHale, and J. L. Pena.  Plus-minus player ratings for soccer.  European Journal of Operational Research, vol. 283, pp. 726–736, 2020.
  * [13] W. Spearman.  Beyond expected goals.  MIT Sloan Sports Analytics Conference, 2018.
  * [14] M. Brechot and R. Flepp.  Dealing with randomness in match outcomes: how to rethink performance evaluation in European club football using expected goals.  Journal of Sports Economics, vol. 21, pp. 335–362, 2020.
  * [15] A. Fairchild, K. Pelechrinis, and M. Kokkodis,  Spatial analysis of shots in MLS: A model for expected goals and fractal dimensionality.  Journal of Sports Analytics, vol. 4, pp. 165–174, 2018.
  * [16] P. Biecek, and T. Burzykowski.  Explanatory Model Analysis.  Chapman and Hall/CRC, New York, 2021.
  * [17] J. H. Friedman.  “Greedy Function Approximation: A Gradient Boosting Machine”.  Annals of Statistics, vol. 29, pp. 1189–1232, 2000.
  * [18] P. Robberechts and J. Davis.  How Data Availability Affects the Ability to Learn Good xG Models.  in Brefeld, U., Davis, J., Van Haaren, J., Zimmermann, A. (eds) Machine Learning and Data Mining for Sports Analytics, MLSA 2020, Communications in Computer and Information Science, vol. 1324, Springer, Cham, 2020.
  * [19] J. Zivkovic, and T. ElHabr.  worldfootballR: Functions to Extract and Clean World Football (Soccer) Data. https://github.com/JaseZiv/worldfootballR, 2022.
  * [20] H. T. Ly, S. Szmajdzinski, and A. Kozak,  forester: Automated Machine Learning Model Solver. https://github.com/ModelOriented/forester, 2022.
  * [21] T. Chen and C. Guestrin.  XGBoost: A Scalable Tree Boosting System.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 785–794, 2016.
  * [22] L. Breiman.  Random forests.  Machine Learning, vol. 45, pp. 5–32, 2001.
  * [23] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T. Liu.  LightGBM: A Highly Efficient Gradient Boosting Decision Tree.  Advances in Neural Information Processing Systems, pp. 3149–3157, 2017.
  * [24] A. V. Dorogush, V. Ershov, and A. Gulin.  CatBoost: gradient boosting with categorical features support.  Workshop on ML Systems at NIPS 2017.
  * [25] P. Biecek.  DALEX: Explainers for Complex Predictive Models in R.  Journal of Machine Learning Research, vol. 19, pp. 1–5, 2018.
  * [26] X. Guo, Y. Yin, C. Dong, G. Yang, and G. Zhou.  On the class imbalance problem.,  In ICNC, pp. 192–-201, 2008.
  * [27] N. Lunardon, G. Menardi, and N. Torelli.  ROSE: A Package for Binary Imbalanced Learning.  The R Journal, vol. 6, pp. 79–89, 2014.
  * [28] G. Menardi and N. Torelli,  Training and assessing classification rules with imbalanced data.  Data Mining and Knowledge Discovery, vol. 28, pp. 92–122, 2014.
  * [29] J. V. Haaren.  Why would I trust your numbers? On the explainability of expected values in soccer.  arXiv preprint, arXiv:2105.13778, 2021.
  * [30] I. Umami, D. H. Gutama, and H. R. Hatta.  Implementing the expected goal (xG) model to predict scores in soccer matches.  International Journal of Informatics and Information Systems, vol. 4, pp. 38–54, 2021.
  * [31] J. Fernandez, L. Bornn, and D. Cervone.  A framework for the fne-grained evaluation of the instantaneous expected value of soccer possessions.  Machine Learning, vol. 110, pp. 1389-–1427, 2021.
  * [32] A. Patil, A. Framewala, and F. Kazi.  Explainability of SMOTE Based Oversampling for Imbalanced Dataset Problems.  3rd International Conference on Information and Computer Technologies, pp. 41–45, 2020.