EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH

10 July 2023

Revised version:
30 August 2023

Search for $K^+$ decays into the $\pi^+e^+e^-e^+e^-$ final state

The NA62 Collaboration

The first search for ultra-rare $K^+$ decays into the $\pi^+e^+e^-e^+e^-$ final state is reported, using a dataset collected by the NA62 experiment at CERN in 2017–2018. An upper limit of  at 90% CL is obtained for the branching ratio of the $K^+\to\pi^+e^+e^-e^+e^-$ decay, predicted in the Standard Model to be $(7.2\pm0.7)\times 10^{-11}$. Upper limits at 90% CL are obtained at the level of $10^{-9}$ for the branching ratios of two prompt decay chains involving pair-production of hidden-sector mediators: $K^+\to\pi^+aa$, $a\to e^+e^-$ and $K^+\to\pi^+S$, $S\to A^\prime A^\prime$, $A^\prime\to e^+e^-$.

Accepted for publication in Physics Letters B

§ INTRODUCTION

Dark-sector models provide plausible dark-matter candidates, and represent a compelling new-physics direction to explore in its own right [1, 2]. So far, searches for the production of dark-sector mediators in meson decays have been focused on the production of a single particle, which is either invisible or decays into lepton or photon pairs. However, the parameter space of minimal dark sectors not yet excluded by experiment accommodates detectable pair-production of dark states [3, 4]. In the kaon sector, the process of particular interest is $K\to\pi XX$ followed by prompt $X\to e^+e^-$ decays, leading to characteristic multi-electron final states. Since this process has not been studied experimentally so far, ${\cal O}(10^{-6})$ sensitivity to its branching ratio is sufficient to improve existing constraints on dark-sector models.

The principal $K^+$ decay channel leading to the $\pi^+e^+e^-e^+e^-$ final state is the $K^+\to\pi^+\pi^0$ decay followed by the double Dalitz decay $\pi^0_{\rm DD}\to e^+e^-e^+e^-$. This decay chain, denoted $K_{2\pi {\rm DD}}$, with a measured branching ratio of $(6.9\pm0.3)\times 10^{-6}$ [5], is used for normalisation in this analysis. The non-resonant part of the  $K^+\to\pi^+e^+e^-e^+e^-$ decay, denoted $K_{\pi4e}$, is expected to occur via one-photon ($K^+\to\pi^+\gamma^*$) and two-photon ($K^+\to\pi^+\gamma^*\gamma^*$) exchange with an expected Standard Model (SM) branching ratio of $(7.2\pm0.7)\times 10^{-11}$ [6]. The resulting final state is of two-fold interest in the context of dark sectors:

* A short-lived QCD axion ($a$) decaying into an $e^+e^-$ pair is not completely ruled out by experiment [7, 8, 9], and provides a plausible explanation for the “17 MeV anomaly” in the mass spectra of the $e^+e^-$ pairs produced in the de-excitation of $^8$Be [10], $^4$He [11] and $^{12}$C [12] nuclei. Assuming an axion mass of 17 MeV/$c^2$, lower bounds of ${\cal B}(K_L\to\pi^0 aa)>10^{-7}$ and ${\cal B}(K^+\to\pi^+aa)>2\times 10^{-8}$ (differing by the $K_L/K^+$ lifetime ratio) are predicted [3], allowing a test of the QCD axion explanation for the 17 MeV anomaly.

* A scenario involving a dark scalar ($S$), and a dark photon ($A^\prime$) with masses satisfying the condition $m_S\ge 2m_{A^\prime}$ leads to a prompt cascade process $K^+\to\pi^+S$, $S\to A^\prime A^\prime$ [3]. The case of prompt $A^\prime\to e^+e^-$ decay leads to the same final state as the $K_{\pi4e}$ decay.

The NA62 experiment at CERN collected a dataset of $K^+$ decays into final states containing lepton pairs in 2017–2018. The first search, performed using this dataset, is reported here for the $K_{\pi4e}$ decay, and for pair-production of hidden-sector mediators in the prompt decay chains $K^+\to\pi^+aa$, $a\to e^+e^-$ (denoted $K_{\pi aa}$) and $K^+\to\pi^+S$, $S\to A^\prime A^\prime$, $A^\prime\to e^+e^-$ (denoted $K_{\pi S}$).

§ BEAM, DETECTOR AND DATA SAMPLE

The layout of the NA62 beamline and detector [13] is shown schematically in Fig. <ref>. An unseparated secondary beam of $\pi^+$ (70%), protons (23%) and $K^+$ (6%) is created by directing 400 GeV/$c$ protons extracted from the CERN SPS onto a beryllium target in spills of 3 s effective duration. The central beam momentum is 75 GeV/$c$, with a momentum spread of 1% (rms).

Beam kaons are tagged with a 70 ps time resolution by a differential Cherenkov counter (KTAG) that uses nitrogen gas at 1.75 bar pressure contained in a 5 m long vessel as the radiator. Beam particle positions, momenta and times

are measured by a silicon pixel spectrometer consisting of three stations (GTK1,2,3) and four dipole magnets forming an achromat; a toroidal muon sweeper (scraper, SCR) is installed between GTK1 and GTK2. A 1.2 m thick steel collimator (COL) with a $76\times40$ mm$^2$ central aperture and $1.7\times1.8$ m$^2$ outer dimensions is placed upstream of GTK3 to absorb hadrons from upstream $K^+$ decays; a variable-aperture collimator of $0.15\times0.15$ m$^2$ outer dimensions was used up to early 2018. A dipole magnet (TRIM5) providing a 90 MeV/c horizontal momentum kick is located in front of GTK3. Inelastic interactions of beam particles in GTK3 are detected by an array of scintillator hodoscopes (CHANTI). The beam is delivered into a vacuum tank evacuated to $10^{-6}$ mbar, which contains a 75 m long fiducial volume (FV) starting 2.6 m downstream of GTK3. The beam angular spread at the FV entrance is 0.11 mrad (rms) in both horizontal and vertical planes. Downstream of the FV, undecayed beam particles continue their path in vacuum.

Momenta of charged particles produced in $K^+$ decays in the FV are measured by a magnetic spectrometer (STRAW) located in the vacuum tank downstream of the FV. The spectrometer consists of four tracking chambers made of straw tubes, and a dipole magnet (M) located between the second and third chambers that provides a horizontal momentum kick of 270 MeV/$c$ in a direction opposite to that produced by TRIM5. The momentum resolution is $\sigma_p/p = (0.30\oplus 0.005\cdot p)\%$, with the momentum $p$ expressed in GeV/$c$.

A ring-imaging Cherenkov detector (RICH) consisting of a 17.5 m long vessel filled with neon at atmospheric pressure provides particle identification, charged particle time measurements with a typical resolution of 70 ps, and the trigger time. Two scintillator hodoscopes (CHOD), which include a matrix of tiles and two planes of slabs arranged in four quadrants located downstream of the RICH, provide trigger signals and time measurements with 200 ps precision.

Schematic side view of the NA62 beamline and detector used in 2018.

A $27X_0$ thick quasi-homogeneous liquid-krypton (LKr) electromagnetic calorimeter is used for particle identification and photon detection. The calorimeter has an active volume of 7 m$^3$, segmented in the transverse direction into 13248 projective cells of $2\times 2$ cm$^2$ size, and provides an energy resolution $\sigma_E/E=(4.8/\sqrt{E}\oplus11/E\oplus0.9)\%$, with $E$ expressed in GeV. To achieve hermetic acceptance for photons emitted in $K^+$ decays in the FV at angles up to 50 mrad from the beam axis, the LKr calorimeter is complemented by annular lead glass detectors (LAV) installed in 12 positions inside and downstream of the vacuum tank, and two lead/scintillator sampling calorimeters (IRC, SAC) located close to the beam axis. An iron/scintillator sampling hadronic calorimeter formed of two modules (MUV1,2) and a muon detector consisting of 148 scintillator tiles located behind an 80 cm thick iron wall (MUV3) are used for particle identification.

The data sample analysed is obtained from $8.3\times 10^5$ SPS spills recorded in 2017–2018, with the typical beam intensity increasing over time from  to  protons per spill. The latter value corresponds to a 500 MHz mean beam particle rate at the FV entrance, and a 3.7 MHz mean $K^+$ decay rate in the FV. The main trigger line is designed for the measurement of the ultra-rare $K^+\to\pi^+\nu\bar\nu$ decay [14]. The analysis presented here
is based on dedicated electron multi-track ($e$MT) and multi-track (MT) trigger lines operating concurrently with the main trigger line [15, 16], typically downscaled by factors of 8 and 100, respectively. Both trigger lines involve a low-level stage and a high-level stage. The low-level (L0) hardware trigger is based on RICH signal multiplicity and coincidence of signals in two diagonally opposite CHOD quadrants, with an additional requirement of at least 20 GeV deposited in the LKr calorimeter in the $e$MT line. The high-level (L1) software trigger involves beam $K^+$ identification by the KTAG and reconstruction of a negatively-charged STRAW track.

Monte Carlo (MC) simulation of particle interactions with the detector and its response is performed using a software package based on the toolkit [17]. Accidental activity and trigger line response are included in the simulation.

§ EVENT SELECTION

The $K_{2\pi {\rm DD}}$ events collected concurrently with signal candidates with the same trigger lines are used for normalisation. Signal and normalisation decay modes have the same set of particles in the final state, leading to a first-order cancellation of detector and trigger inefficiencies.

Five-track decays often produce final states with multiple tracks not fully contained in the geometric acceptance of the downstream detectors. The effect is particularly pronounced for the $K_{\pi 4e}$ and $K_{2\pi {\rm DD}}$ decays which typically produce soft, almost collinear electrons. For example, requiring a fully reconstructed five-track vertex leads to an acceptance of 1.2% for the $K_{2\pi {\rm DD}}$ decay, and requiring track separation in the STRAW1 plane decreases the acceptance by another order of magnitude, which determines the overall acceptance level.

Only 11% of the $K_{2\pi {\rm DD}}$ decays with a fully reconstructed five-track vertex have all tracks in the CHOD's geometric acceptance. To avoid an order-of-magnitude loss of acceptance, the analysis relies on the STRAW information only, and the STRAW tracks are not required to be in the geometric acceptance of any downstream detectors. The reduction in selectivity caused by the lack of information from the downstream detectors is mitigated by the exploitation

of kinematic constraints on the five-track final state. Main selection criteria are described below.

Five-track vertices compatible with $K_{\pi4e}$ and $K_{2\pi\rm DD}$ decays are selected as follows:

* Vertices are reconstructed by extrapolating STRAW tracks backward, taking into account the measured residual magnetic field in the vacuum tank. Exactly one five-track vertex should be present in the event within the FV. The longitudinal position of the vertex along the beam axis, $z_{\rm vtx}$, should be at least 5 m downstream of the start of the FV.

The total electric charge of the five tracks should be $q=+1$, their momenta should be in the range 5–45 GeV/$c$, and their trajectories through the STRAW chambers should be within the respective geometric acceptances. To suppress photon conversions and fake tracks, each pair of tracks should be separated by at least 15 mm in each STRAW chamber plane.

* The momentum excess defined as $\Delta p=p_{\rm vtx}-p_{\rm beam}$, where $p_{\rm vtx}$ is the total momentum of the five tracks and $p_{\rm beam}$ is the central beam momentum, should satisfy $|\Delta p|<2~{\rm GeV}/c$. The total transverse momentum of the tracks with respect to the beam axis should be $p_T<25~{\rm MeV}/c$. The quantity $p_{\rm beam}$ and the beam axis direction are monitored throughout the data taking, typically every few hours, with $K^+\to\pi^+\pi^+\pi^-$ decays.

* Track times are evaluated using the STRAW information, with a resolution of 6 ns. Each track is required to be within 20 ns of the trigger time.

Particle identification is based on event kinematics, and is applied as follows:

* The $K^+\to\pi^+e^+e^-e^+e^-$ final state involves three positively charged tracks (one $\pi^+$ and two $e^+$), and two negatively charged tracks (two $e^-$). Three assignments of the $\pi^+$ mass to one of the positively charged tracks are considered. In each mass assignment, the five-track mass, $m_{\pi 4e}$, the four-electron mass, $m_{4e}$, and the squared missing mass, $m_{\rm miss}^2=(P_K-P_\pi)^2$, are evaluated, where $P_K$ is the kaon four-momentum computed using the central beam momentum, and $P_\pi$ is the reconstructed $\pi^+$ four-momentum.

* The mass assignment corresponding to the minimal value of $|m_{\pi 4e}-m_K|$ is chosen, where $m_K$ is the $K^+$ mass, and it is required that $484~{\rm MeV}/c^2<m_{\pi4e}<504~{\rm MeV}/c^2$. The resolution of $m_{\pi4e}$ is 1.7 MeV/$c^2$. The probability of incorrectly assigning the mass is determined by simulation and found to be 4.5% for $K_{\pi 4e}$ decays after the full selection, and ${\cal O}(10^{-4})$ for $K_{2\pi {\rm DD}}$ decays after the full selection which involves a $\pi^0$ mass constraint.

The following condition is specific to $K_{2\pi\rm DD}$ selection:

* For the mass assignment chosen, it is required that $|m_{4e}-m_{\pi^0}|<10~{\rm MeV}/c^2$, where $m_{\pi^0}$ is the $\pi^0$ mass. The resolution of $m_{4e}$ is 0.9 MeV/$c^2$.

The following criteria are specific to $K_{\pi 4e}$ selection:

* To reject events consistent with a $K_{2\pi{\rm DD}}$ decay in any mass assignment, the requirements $|m_{4e}-m_{\pi^0}|>10~{\rm MeV}/c^2$, $m_{\rm miss}^2>0$ and $|m_{\rm miss}-m_{\pi^0}|>40~{\rm MeV}/c^2$ are applied for all three mass assignments. This reduces the $K_{\pi4e}$ acceptance by a factor of 0.74.

* To suppress the $K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm D}$ background, which predominantly enters by incorrect mass assignment (with a soft $e^+$ from a $\pi^0_{\rm D}\to\gamma e^+e^-$ decay considered as the $\pi^+$ candidate), the reconstructed $\pi^+$ momentum is required to be $p_\pi>10~{\rm GeV}/c$. This reduces the $K_{\pi4e}$ acceptance by a factor of 0.93.

Selection of the $K_{\pi aa}$ and $K_{\pi S}$ decay chains consists of the $K_{\pi 4e}$ selection with the following additional criteria:

* Two $X\to e^+e^-$ decay hypotheses are examined for the four $e^\pm$ candidates in the mass assignment chosen. In each hypothesis, the masses of the two $e^+e^-$ pairs, $m_{ee1}$ and $m_{ee2}$, are evaluated, and a discriminant is computed to quantify the consistency of the two mass values:

D=(m_ee1-m_ee2)^2   /   (4.9×10^-3·m_ee)^2.

Here $m_{ee} = (m_{ee1}+m_{ee2})/2$ is the reconstructed mass of the intermediate particle $X$ which is pair-produced, and the denominator represents the squared resolution of the quantity $m_{ee1}-m_{ee2}$, evaluated with simulation and found to be proportional to $m_{ee}$. The $X\to e^+e^-$ hypothesis with the smaller ${\cal D}$ value is considered, and it is required that ${\cal D}<9$ in this hypothesis.

* For each $X$ mass hypothesis ($m_X$), it is required that $|m_{ee}-m_X|<0.02\cdot m_X$.

The signal region in the data is kept masked until the completion and validation of the background evaluation.

§ NUMBER OF KAON DECAYS IN THE FIDUCIAL VOLUME

The reconstructed $m_{\pi 4e}$ distributions of $K_{2\pi{\rm DD}}$ candidates from both data and simulation are displayed in Fig. <ref>. The simulation reproduces the gaussian part of the $m_{\pi 4e}$ spectrum but does not fully describe the radiative tail; this is due to a soft-photon bremsstrahlung approximation [18] used to model the radiative corrections to the $\pi^0_{\rm DD}$ decay.

The number of $K^+$ decays in the FV is computed as

N_K = f ·N_DD/B(K_2π) ·B(π^0_DD) ·A_DD
= (8.58 ±0.19_stat ±0.07_MC ±0.41_ext)×10^11,

where $N_{\rm DD}=2023$ is the number of $K_{2\pi{\rm DD}}$ candidates in the data;
${\cal B}(K_{2\pi})=(20.67\pm0.08)\times 10^{-2}$
and ${\cal B}(\pi^0_{\rm DD})=(3.34\pm0.16)\times 10^{-5}$ are the $K^+\to\pi^+\pi^0$ and $\pi^0_{\rm DD}$ decay branching ratios [5]; $A_{\rm DD}=(3.41 \pm 0.02_{\rm stat}\pm 0.02_{\rm syst})\times 10^{-4}$ is the acceptance of the $K_{2\pi{\rm DD}}$ selection evaluated with simulation (the systematic error is evaluated by variation of the $m_{\pi4e}$ selection conditions); $f=0.9987$ is the purity of the $K_{2\pi{\rm DD}}$ sample evaluated with simulation  and accounting for the $K^+\to\pi^0_{\rm DD}\ell^+\nu$ background ($\ell=e,\mu$). The three uncertainties quoted in $N_K$ arise from the limited size of the data sample, the total uncertainty in the acceptance $A_{\rm DD}$, and the external uncertainties in the measured values of ${\cal B}(K_{2\pi})$ and ${\cal B}(\pi^0_{\rm DD})$.

The acceptance $A_{\rm DD}$ includes the effects of accidental activity and trigger efficiencies. The efficiencies of the CHOD and LKr trigger conditions for the $K_{2\pi\rm DD}$ sample are found with simulation to be 74% and 81% respectively, and are relatively low due to the absence of the CHOD and LKr geometric acceptance requirements for the tracks. The efficiencies of the RICH, KTAG and STRAW trigger conditions exceed 99%.

Reconstructed $m_{\pi4e}$ spectra for the data and a simulated $K_{2\pi \rm DD}$ sample obtained with the $K_{2\pi \rm DD}$ selection. The signal mass region is indicated by vertical arrows.

§ BACKGROUND ESTIMATION

Background from five-track and seven-track $K^+$ decays is evaluated with simulation; all $K^+$ decays with branching ratios exceeding $10^{-8}$ are considered. Background from coincidences of pairs of three-track $K^+$ decays in time (with one track not reconstructed) is evaluated with dedicated simulation of pairs of decays occurring within 25 ns; coincidences of all types of three-track $K^+$ decays with products of branching ratios exceeding $2\times 10^{-6}$ are considered.

The probability of two decays occurring within this 25 ns time window is dependent upon the instantaneous beam intensity, and is accounted for by normalising the total estimated background to the data in the momentum excess region .

Contributions from all possible triple coincidences of one or more $K^+\to\pi^+\pi^+\pi^-$ decays with $K^+\to\mu^+\nu$ decays and beam halo particles, and contributions from coincidences of three-track decays with inelastic interactions of beam particles in GTK3, are found to be much smaller than those from coincidences of pairs of three-track decays, and are neglected.

The reconstructed momentum excess ($\Delta p$) and five-track mass ($m_{\pi4e}$) spectra for data and simulated events after the $K_{\pi 4e}$ selection are shown in Fig. <ref>. The region of large $\Delta p$ values is dominated by backgrounds from coincidences of three-track decays. The integrals of the data and simulated spectra agree in the region $\Delta p>20~{\rm GeV}/c$ by construction. The region $\Delta p<0$ is dominated by backgrounds from single five-track decays which, with the exception of the $K_{2\pi\rm DD}$ decay, lead to photons and neutrinos in the final state.

The background model is validated against the data in the control region $\Delta p<-2~{\rm GeV}/c$ using the standard $K_{\pi4e}$ selection (Fig. <ref>, left) and a loose $K_{\pi4e}$ selection with the condition  removed. The results are shown in Table <ref>. The loose selection leads to a sharp increase in the $K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm D}$ background in the control region due to the soft momentum spectrum of the positrons from $\pi^0_{\rm D}$ decays misidentified as $\pi^+$ (Section <ref>).

Reconstructed momentum excess ($\Delta p$, left) and five-track mass ($m_{\pi 4e}$, right) spectra for data and simulated samples after the $K_{\pi 4e}$ selection obtained by removing the $\Delta p$ condition (left panel) or the $m_{\pi 4e}$ condition (right panel).

Contributions from three-track $K^+$ decays (the last six entries in the legend) are due to coincidences with a $K^+\to\pi^+\pi^+\pi^-$ decay. The sum of these contributions is normalised to the data in the region $\Delta p>20~{\rm GeV}/c$. The signal region is shown with pairs of vertical arrows in both plots.

Backgrounds to the $K_{\pi 4e}$ decay, their branching ratios (or products of branching ratios for coincidences of two decays), and estimated backgrounds with their statistical uncertainties in the control ($\Delta p<-2~{\rm GeV}/c$) and signal ($|\Delta p|<2~{\rm GeV}/c$) momentum excess regions shown in Fig. <ref> (left). Estimated backgrounds in the control region for a loose $K_{\pi4e}$ selection, with the $p_\pi>10~{\rm GeV}/c$ condition removed, are also shown. The $K^+\to\pi^+\pi^0_{\rm D}e^+e^-$ decay, as well as the $K^+\to\pi^+\pi^-e^+\nu$ decay in coincidence with a $K^+\to\pi^+\pi^+\pi^-$ decay, and coincidences of any of the two decays, $K^+\to\pi^+\pi^+\pi^-$ and $K^+\to\pi^+\pi^0_{\rm D}$, lead to negligible backgrounds in the control and signal regions, and are not listed.

Source     Branching ratio     Control region     Control region,     Signal region

(or their product)         loose selection

4lSingle five-track and seven-track decays

$K_{2\pi{\rm DD}}$     $6.9\times10^{-6}$     $0.06\pm0.06$     $0.06\pm0.06$     –

$K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm D}$     $2.4\times10^{-6}$     $0.30\pm0.06$
$2.47\pm0.16$     $0.04\pm0.02$

$K^+\to\pi^0_{\rm DD}e^+\nu$     $1.7\times10^{-6}$     $0.10\pm0.05$     $0.10\pm0.05$     –

$K^+\to\pi^+\pi^0\pi^0_{\rm DD}$     $1.2\times10^{-6}$     $0.03\pm0.03$     $0.03\pm0.03$     –

$K^+\to\pi^0_{\rm DD}\mu^+\nu$     $1.1\times10^{-6}$     $0.02\pm0.02$     $0.03\pm0.02$     –

$K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm DD}$     $1.4\times10^{-8}$     $0.05\pm0.02$
$0.10\pm0.02$     $0.01\pm0.01$

4lCoincidences of three-track decays with a $K^+\to\pi^+\pi^+\pi^-$ decay

$K^+\to\pi^0_{\rm D}e^+\nu$     $3.3\times10^{-5}$     $0.15\pm0.07$     $0.15\pm0.07$     $0.08\pm0.05$

$K^+\to\pi^+\pi^0\pi^0_{\rm D}$     $2.3\times10^{-5}$     $0.03\pm0.03$     $0.08\pm0.05$     –

$K^+\to\pi^0_{\rm D}\mu^+\nu$     $2.2\times10^{-5}$     $0.03\pm0.02$     $0.04\pm0.02$     $0.05\pm0.02$

Total        $0.77\pm0.13$     $3.06\pm0.21$     $0.18\pm0.06$

Data        1     4     0

The background model is further validated using the $m_{\pi4e}$ spectrum outside the signal region. Using the standard $K_{\pi4e}$ selection, three data events are observed with an expected background of $2.83\pm0.64_{\rm stat}$ (Fig. <ref>, right); using the loose $K_{\pi4e}$ selection, five data events are observed with an expected background of $3.99\pm0.74_{\rm stat}$. In all cases, the data are in agreement with simulation within statistical fluctuations, and the largest background comes from the $K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm D}$ decay.

The background contributions in the signal region after the $K_{\pi4e}$ selection are summarised in Table <ref>. The estimate of the total background, validated by agreement of the data with simulation in the control regions, is

N_B = 0.18±0.06_stat ±0.13_syst = 0.18 ±0.14,

where a 100% relative systematic uncertainty is assigned to the contribution from coincidences of three-track decays, based on the level of agreement of data and simulation in the momentum excess region $2~{\rm GeV}/c<\Delta p<20~{\rm GeV}/c$.

The requirement that the two reconstructed $e^+e^-$ masses ($m_{ee1}$, $m_{ee2}$) are consistent with each other, applied in the selection of the $K_{\pi aa}$ and $K_{\pi S}$ decay chains, leads to a reduction of the background with respect to the $K_{\pi4e}$ selection. The probability of the reconstructed $m_{ee1}$ and $m_{ee2}$ values being consistent is found to be 2% in simulated samples of $K^+\to\pi^+\pi^0_{\rm D}\pi^0_{\rm D}$ and $K^+\to\pi^+\pi^0_{\rm DD}$ decays. Furthermore, the probability of $|m_{ee}-m_X|<0.02\cdot m_X$ for the above samples does not exceed 10% in any $m_X$ hypothesis considered (where $X=a,A'$) as established with simulations. The background to the $K_{\pi aa}$ and $K_{\pi S}$ processes in each $m_X$ hypothesis is estimated as

N_B^'= (0.4±0.4)×10^-3,

where a 100% relative uncertainty is conservatively assigned.

§ SEARCH FOR THE SIGNAL DECAYS

The $K_{\pi4e}$ decay is modelled according to the SM description [6]. The non-resonant part of the decay amplitude is considered, including the dominant one-photon exchange diagram ($K^+\to\pi^+\gamma^*$), the two-photon exchange diagram ($K^+\to\pi^+\gamma^*\gamma^*$), and their interference. The contribution from the $\pi^0$ pole in the two-photon exchange is neglected; it does not interfere with the other amplitudes, peaks at $m_{4e}\approx m_{\pi^0}$, and accounts for 3% of the total decay rate in the region $|m_{4e}-m_{\pi^0}|>10~{\rm MeV}/c^2$ used in this analysis. The acceptance of the $K_{\pi 4e}$ selection for the $K_{\pi4e}$ signal is found to be $A_{\pi 4e}=(1.85\pm0.06_{\rm stat})\times 10^{-4}$.

The prompt decay chain $K^+\to\pi^+aa$, $a\to e^+e^-$, allowed kinematically for axion masses in the range $2m_e \le m_a \le (m_K-m_\pi)/2$, is simulated with a uniform $K^+\to\pi^+aa$ phase-space distribution [3] and isotropic $a\to e^+e^-$ decays, for 33 equally-spaced $m_a$ hypotheses in the range 10–170 MeV/$c^2$. The acceptance of the $K_{\pi aa}$ selection reaches a maximum of $7.1\times 10^{-3}$ for $m_a=155~{\rm MeV}/c^2$.

The prompt decay chain $K^+\to\pi^+S$, $S\to A^\prime A^\prime$, $A^\prime\to e^+e^-$, allowed kinematically in the mass region $2m_e \le m_{A^\prime} \le m_S/2 \le (m_K-m_\pi)/2$, is simulated as a series of isotropic decays for $m_{A^\prime}$ hypotheses in the range 10–170 MeV/$c^2$ and $m_S$ hypotheses in the range 20–340 MeV/$c^2$, with a 5 MeV/$c^2$ step in both $m_{A^\prime}$ and $m_S$. The region $120~{\rm MeV}/c^2<m_S<165~{\rm MeV}/c^2$ is excluded from the search as the signal acceptance is suppressed in this region by the selection conditions that reduce the $K_{2\pi\rm DD}$ background. The acceptance of the $K_{\pi S}$ selection reaches a maximum of $1.1\times 10^{-2}$ at the kinematic boundary for $m_{A^\prime}=150~{\rm MeV}/c^2$ and $m_S=300~{\rm MeV}/c^2$.

The $K_{\pi aa}$ and $K_{\pi S}$ decays typically produce energetic, well-separated electrons in the final state, which leads to significantly larger acceptances than for $K_{\pi 4e}$ and $K_{2\pi\rm DD}$ decays.

§ RESULTS

No data events are observed in the signal region after unmasking. In this case, upper limits at 90% CL on the signal branching ratios computed using the CL$_{\rm S}$ method [19] do not depend on the expected background and are equal to $2.3\,/(A\,N_K)$, where $A$ denotes the signal acceptance. The uncertainties in $N_K$ and signal acceptances are taken into account in the CL$_{\rm S}$ procedure.
An upper limit of ${\cal B}(K_{\pi4e}) < 1.4\times 10^{-8}$ is established at 90% CL for the branching ratio of the non-resonant part of the $K_{\pi4e}$ decay, which is a factor 200 larger than the SM expectation [6].

Upper limits obtained at 90% CL for the branching ratios of the $K_{\pi aa}$ and $K_{\pi S}$ decay chains are shown as functions of the assumed masses of the dark-sector mediators in Fig. <ref>. For the product ${\cal B}(K^+\to\pi^+aa)\times\left[{\cal B}(a\to e^+e^-)\right]^2$, the most stringent upper limit of $3.7\times 10^{-10}$ is obtained for $m_a=155~{\rm MeV}/c^2$. The upper limit of $2.1\times 10^{-9}$ obtained for the above quantity for $m_a=17~{\rm MeV}/c^2$ excludes the QCD axion explanation of the “17 MeV anomaly” in the mass spectra of the $e^+e^-$ pairs produced in nuclear de-excitation [10, 11, 12], according to the decay rate estimate reported in Ref. [3]. For the product ${\cal B}(K^+\to\pi^+S)\times{\cal B}(S\to A^\prime A^\prime)\times\left[{\cal B}(A^\prime\to e^+e^-)\right]^2$, the most stringent upper limit of $2.5\times 10^{-10}$ is obtained for $m_{A^\prime}=150~{\rm MeV}/c^2$ and $m_S=300~{\rm MeV}/c^2$.

Left: upper limit at 90% CL of the branching ratio of the prompt decay chain $K^+\to\pi^+aa$, $a\to e^+e^-$ as a function of the assumed axion mass. Right: upper limit at 90% CL of the branching ratio of the prompt decay chain $K^+\to\pi^+S$, $S\to A^\prime A^\prime$, $A^\prime\to e^+e^-$ as a function of the assumed dark-photon and dark-scalar masses. Discrete sets of points obtained with a step size of 5 MeV/$c^2$ are approximated by polynomial functions in each case.

§ SUMMARY

The first search for ultra-rare $K^+$ decays into the $\pi^+e^+e^-e^+e^-$ final state is reported, using a dataset collected by the NA62 experiment at CERN in 2017–2018. An upper limit of $1.4\times 10^{-8}$ at 90% CL is obtained for the branching ratio of the $K^+\to\pi^+e^+e^-e^+e^-$ decay, with the expected SM branching ratio of $(7.2\pm0.7)\times 10^{-11}$. Upper limits at 90% CL are obtained at the level of $10^{-9}$ for the branching ratios of two prompt decay chains involving pair-production of hidden-sector mediators: $K^+\to\pi^+aa$, $a\to e^+e^-$ and $K^+\to\pi^+S$, $S\to A^\prime A^\prime$, $A^\prime\to e^+e^-$. The QCD axion is excluded as a possible explanation of the “17 MeV anomaly” according to the decay rate estimate reported in Ref. [3].

§ ACKNOWLEDGEMENTS

We are grateful to Matheus Hostert and Maxim Pospelov for suggesting this study, and for the discussions of the phenomenology of pair-production of dark states.
It is a pleasure to express our appreciation to the staff of the CERN laboratory and the technical
staff of the participating laboratories and universities for their efforts in the operation of the
experiment and data processing.

The cost of the experiment and its auxiliary systems was supported by the funding agencies of
the Collaboration Institutes. We are particularly indebted to:
F.R.S.-FNRS (Fonds de la Recherche Scientifique - FNRS), under Grants No. 4.4512.10, 1.B.258.20, Belgium;
CECI (Consortium des Equipements de Calcul Intensif), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region, Belgium;

NSERC (Natural Sciences and Engineering Research Council), funding SAPPJ-2018-0017,  Canada;

MEYS (Ministry of Education, Youth and Sports) funding LM 2018104, Czech Republic;
BMBF (Bundesministerium für Bildung und Forschung) contracts 05H12UM5, 05H15UMCNA and 05H18UMCNA, Germany;
INFN  (Istituto Nazionale di Fisica Nucleare),  Italy;
MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca),  Italy;
CONACyT  (Consejo Nacional de Ciencia y Tecnología),  Mexico;
IFA (Institute of Atomic Physics) Romanian

CERN-RO No. 1/16.03.2016

and Nucleus Programme PN 19 06 01 04,  Romania;

MESRS  (Ministry of Education, Science, Research and Sport), Slovakia;
CERN (European Organization for Nuclear Research), Switzerland;
STFC (Science and Technology Facilities Council), United Kingdom;
NSF (National Science Foundation) Award Numbers 1506088 and 1806430,  U.S.A.;
ERC (European Research Council)  “UniversaLepto” advanced grant 268062, “KaonLepton” starting grant 336581, Europe.

Individuals have received support from:
Charles University Research Center (UNCE/SCI/013), Czech Republic;
Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR  “Futuro in ricerca 2012”  grant RBFR12JF2Z, Project GAP), Italy;

the Royal Society  (grants UF100308, UF0758946), United Kingdom;
STFC (Rutherford fellowships ST/J00412X/1, ST/M005798/1), United Kingdom;
ERC (grants 268062,  336581 and  starting grant 802836 “AxScale”);
EU Horizon 2020 (Marie Skłodowska-Curie grants 701386, 754496, 842407, 893101, 101023808).

[1]
J. Beacham et al., J. Phys. G47 (2020) 010501.

[2]
P. Agrawal et al., Eur. Phys. J. C81 (2021) 11.

[3]
M. Hostert and M. Pospelov, Phys. Rev. D105 (2022) 015017.

[4]
E. Goudzovski et al., Rep. Prog. Phys. 86 (2023) 016201.

[5]
R.L. Workman et al., Prog. Theor. Exp. Phys. 2022 (2022) 083C01.

[6]
T. Husek, Phys. Rev. D106 (2022) L071301.

[7]
D. Alves and N. Weiner, JHEP 07 (2018) 092.

[8]
D. Alves, Phys. Rev. D103 (2021) 055018.

[9]
D. Alves et al., Eur. Phys. J. C83 (2023) 230.

[10]
A.J. Krasznahorkay et al., Phys. Rev. Lett. 116 (2016) 042501.

[11]
A.J. Krasznahorkay et al., Phys. Rev. C104 (2021) 044003.

[12]
A.J. Krasznahorkay et al., Phys. Rev. C106 (2022) L061601.

[13]
E. Cortina Gil et al., JINST 12 (2017) P05025.

[14]
E. Cortina Gil et al., JHEP 06 (2021) 93.

[15]
R. Ammendola et al., Nucl. Instrum. Methods A929 (2019) 1.

[16]
E. Cortina Gil et al., JHEP 03 (2023) 122.

[17]
J. Allison et al., Nucl. Instrum. Methods A835 (2016) 186.

[18]
K. Kampf, J. Novotný and P. Sanchez-Puertas, Phys. Rev. D97 (2018) 056010.

[19]
A.L. Read, J. Phys. G28 (2002) 2693.

The NA62 Collaboration

Université Catholique de Louvain, Louvain-La-Neuve, Belgium

E. Cortina Gil0000-0001-9627-699X,
A. Kleimenova$\,$[1]0000-0002-9129-4985,
E. Minucci$\,$[2]0000-0002-3972-6824,
S. Padolski0000-0002-6795-7670,
P. Petrov,
A. Shaikhiev$\,$[3]0000-0003-2921-8743,
R. Volpe$\,$[4]0000-0003-1782-2978

TRIUMF, Vancouver, British Columbia, Canada

T. Numao0000-0001-5232-6190,
Y. Petrov0000-0003-2643-8740,
B. Velghe0000-0002-0797-8381,
V. W. S. Wong0000-0001-5975-8164

University of British Columbia, Vancouver, British Columbia, Canada

D. Bryman$\,$[5]0000-0002-9691-0775,
J. Fu

Charles University, Prague, Czech Republic

Z. Hives0000-0002-5025-993X,
T. Husek0000-0002-7208-9150,
J. Jerhot$\,$[6]0000-0002-3236-1471,
K. Kampf0000-0003-1096-667X,
M. Zamkovsky$\,$[7]0000-0002-5067-4789

Aix Marseille University, CNRS/IN2P3, CPPM, Marseille, France

B. De Martino0000-0003-2028-9326,
M. Perrin-Terrin0000-0002-3568-1956

Institut für Physik and PRISMA Cluster of Excellence, Universität Mainz, Mainz, Germany

A.T. Akmete0000-0002-5580-5477,
R. Aliberti$\,$[8]0000-0003-3500-4012,
G. Khoriauli$\,$[9]0000-0002-6353-8452,
J. Kunze,
D. Lomidze$\,$[10]0000-0003-3936-6942,
L. Peruzzo0000-0002-4752-6160,
M. Vormstein,
R. Wanke0000-0002-3636-360X

Dipartimento di Fisica e Scienze della Terra dell'Università e INFN, Sezione di Ferrara, Ferrara, Italy

P. Dalpiaz,
M. Fiorini0000-0001-6559-2084,
A. Mazzolari0000-0003-0804-6778,
I. Neri0000-0002-9669-1058,
A. Norton$\,$[11]0000-0001-5959-5879,
F. Petrucci0000-0002-7220-6919,
M. Soldani0000-0003-4902-943X,
H. Wahl$\,$[12]0000-0003-0354-2465

INFN, Sezione di Ferrara, Ferrara, Italy

L. Bandiera0000-0002-5537-9674,
A. Cotta Ramusino0000-0003-1727-2478,
A. Gianoli0000-0002-2456-8667,
M. Romagnoni0000-0002-2775-6903,
A. Sytov0000-0001-8789-2440

Dipartimento di Fisica e Astronomia dell'Università e INFN, Sezione di Firenze, Sesto Fiorentino, Italy

E. Iacopini0000-0002-5605-2497,
G. Latino0000-0002-4098-3502,
M. Lenti0000-0002-2765-3955,
P. Lo Chiatto0000-0002-4177-557X,
I. Panichi0000-0001-7749-7914,
A. Parenti0000-0002-6132-5680

INFN, Sezione di Firenze, Sesto Fiorentino, Italy

A. Bizzeti$\,$[13]0000-0001-5729-5530,
F. Bucci0000-0003-1726-3838

Laboratori Nazionali di Frascati, Frascati, Italy

A. Antonelli0000-0001-7671-7890,
G. Georgiev$\,$[14]0000-0001-6884-3942,
V. Kozhuharov$\,$[14]0000-0002-0669-7799,
G. Lanfranchi0000-0002-9467-8001,
S. Martellotti0000-0002-4363-7816,
M. Moulson0000-0002-3951-4389,
T. Spadaro0000-0002-7101-2389,
G. Tinti0000-0003-1364-844X

Dipartimento di Fisica “Ettore Pancini” e INFN, Sezione di Napoli, Napoli, Italy

F. Ambrosino0000-0001-5577-1820,
T. Capussela,
M. Corvino0000-0002-2401-412X,
M. D'Errico0000-0001-5326-1106,
D. Di Filippo0000-0003-1567-6786,
R. Fiorenza$\,$[15]0000-0003-4965-7073,
R. Giordano0000-0002-5496-7247,
P. Massarotti0000-0002-9335-9690,
M. Mirra0000-0002-1190-2961,
M. Napolitano0000-0003-1074-9552,
I. Rosa0009-0002-7564-182,
G. Saracino0000-0002-0714-5777

Dipartimento di Fisica e Geologia dell'Università e INFN, Sezione di Perugia, Perugia, Italy

G. Anzivino0000-0002-5967-0952,
F. Brizioli$\,$[7]0000-0002-2047-441X,
E. Imbergamo,
R. Lollini0000-0003-3898-7464,
R. Piandani$\,$[16]0000-0003-2226-8924,
C. Santoni0000-0001-7023-7116

INFN, Sezione di Perugia, Perugia, Italy

M. Barbanera0000-0002-3616-3341,
P. Cenci0000-0001-6149-2676,
B. Checcucci0000-0002-6464-1099,
P. Lubrano0000-0003-0221-4806,
M. Lupi$\,$[17]0000-0001-9770-6197,
M. Pepe0000-0001-5624-4010,
M. Piccini0000-0001-8659-4409

Dipartimento di Fisica dell'Università e INFN, Sezione di Pisa, Pisa, Italy

F. Costantini0000-0002-2974-0067,
L. Di Lella$\,$[12]0000-0003-3697-1098,
N. Doble$\,$[12]0000-0002-0174-5608,
M. Giorgi0000-0001-9571-6260,
S. Giudici0000-0003-3423-7981,
G. Lamanna0000-0001-7452-8498,
E. Lari0000-0003-3303-0524,
E. Pedreschi0000-0001-7631-3933,
M. Sozzi0000-0002-2923-1465

INFN, Sezione di Pisa, Pisa, Italy

C. Cerri,
R. Fantechi0000-0002-6243-5726,
L. Pontisso$\,$[18]0000-0001-7137-5254,
F. Spinella0000-0002-9607-7920

Scuola Normale Superiore e INFN, Sezione di Pisa, Pisa, Italy

I. Mannelli0000-0003-0445-7422

Dipartimento di Fisica, Sapienza Università di Roma e INFN, Sezione di Roma I, Roma, Italy

G. D'Agostini0000-0002-6245-875X,
M. Raggi0000-0002-7448-9481

INFN, Sezione di Roma I, Roma, Italy

A. Biagioni0000-0001-5820-1209,
P. Cretaro0000-0002-2229-149X,
O. Frezza0000-0001-8277-1877,
E. Leonardi0000-0001-8728-7582,
A. Lonardo0000-0002-5909-6508,
M. Turisini0000-0002-5422-1891,
P. Valente0000-0002-5413-0068,
P. Vicini0000-0002-4379-4563

INFN, Sezione di Roma Tor Vergata, Roma, Italy

R. Ammendola0000-0003-4501-3289,
V. Bonaiuto$\,$[19]0000-0002-2328-4793,
A. Fucci,
A. Salamon0000-0002-8438-8983,
F. Sargeni$\,$[20]0000-0002-0131-236X

Dipartimento di Fisica dell'Università e INFN, Sezione di Torino, Torino, Italy

R. Arcidiacono$\,$[21]0000-0001-5904-142X,
B. Bloch-Devaux0000-0002-2463-1232,
M. Boretto$\,$[7]0000-0001-5012-4480,
E. Menichetti0000-0001-7143-8200,
E. Migliore0000-0002-2271-5192,
D. Soldi0000-0001-9059-4831

INFN, Sezione di Torino, Torino, Italy

C. Biino0000-0002-1397-7246,
A. Filippi0000-0003-4715-8748,
F. Marchetto0000-0002-5623-8494

Instituto de Física, Universidad Autónoma de San Luis Potosí, San Luis Potosí, Mexico

A. Briano Olvera0000-0001-6121-3905,
J. Engelfried0000-0001-5478-0602,
N. Estrada-Tristan$\,$[22]0000-0003-2977-9380,
M. A. Reyes Santos$\,$[22]0000-0003-1347-2579

Horia Hulubei National Institute for R&D in Physics and Nuclear Engineering, Bucharest-Magurele, Romania

P. Boboc0000-0001-5532-4887,
A. M. Bragadireanu,
S. A. Ghinescu0000-0003-3716-9857,
O. E. Hutanu

Faculty of Mathematics, Physics and Informatics, Comenius University, Bratislava, Slovakia

L. Bician$\,$[23]0000-0001-9318-0116,
T. Blazek0000-0002-2645-0283,
V. Cerny0000-0003-1998-3441,
Z. Kucerova$\,$[7]0000-0001-8906-3902

CERN, European Organization for Nuclear Research, Geneva, Switzerland

J. Bernhard0000-0001-9256-971X,
A. Ceccucci0000-0002-9506-866X,
M. Ceoletta0000-0002-2532-0217,
H. Danielsson0000-0002-1016-5576,
N. De Simone$\,$[24],
F. Duval,
B. Döbrich$\,$[25]0000-0002-6008-8601,
L. Federici0000-0002-3401-9522,
E. Gamberini0000-0002-6040-4985,
L. Gatignon$\,$[3]0000-0001-6439-2945,
R. Guida,
F. Hahn$\,$[2],
E. B. Holzer0000-0003-2622-6844,
B. Jenninger,
M. Koval$\,$[23]0000-0002-6027-317X,
P. Laycock$\,$[26]0000-0002-8572-5339,
G. Lehmann Miotto0000-0001-9045-7853,
P. Lichard0000-0003-2223-9373,
A. Mapelli0000-0002-4128-1019,
R. Marchevski$\,$[1]0000-0003-3410-0918,
K. Massri0000-0001-7533-6295,
M. Noy,
V. Palladino0000-0002-9786-9620,
J. Pinzino$\,$[27]0000-0002-7418-0636,
V. Ryjov,
S. Schuchmann0000-0002-8088-4226,
S. Venditti

School of Physics and Astronomy, University of Birmingham, Birmingham, United Kingdom

T. Bache0000-0003-4520-830X,
M. B. Brunetti$\,$[28]0000-0003-1639-3577,
V. Duk$\,$[4]0000-0001-6440-0087,
V. Fascianelli$\,$[29],
J. R. Fry0000-0002-3680-361X,
F. Gonnella0000-0003-0885-1654,
E. Goudzovski$\,$[1]0000-0001-9398-4237,
J. Henshaw0000-0001-7059-421X,
L. Iacobuzio,
C. Kenworthy0009-0002-8815-0048,
C. Lazzeroni0000-0003-4074-4787,
N. Lurkin$\,$[6]0000-0002-9440-5927,
F. Newson,
C. Parkinson0000-0003-0344-7361,
A. Romano0000-0003-1779-9122,
J. Sanders0000-0003-1014-094X,
A. Sergi$\,$[30]0000-0001-9495-6115,
A. Sturgess0000-0002-8104-5571,
J. Swallow$\,$[7]0000-0002-1521-0911,
A. Tomczak0000-0001-5635-3567

School of Physics, University of Bristol, Bristol, United Kingdom

H. Heath0000-0001-6576-9740,
R. Page,
S. Trilov0000-0003-0267-6402

School of Physics and Astronomy, University of Glasgow, Glasgow, United Kingdom

B. Angelucci,
D. Britton0000-0001-9998-4342,
C. Graham0000-0001-9121-460X,
D. Protopopescu0000-0002-3964-3930

Faculty of Science and Technology, University of Lancaster, Lancaster, United Kingdom

J. Carmignani$\,$[31]0000-0002-1705-1061,
J. B. Dainton,
R. W. L. Jones0000-0002-6427-3513,
G. Ruggiero$\,$[32]0000-0001-6605-4739

School of Physical Sciences, University of Liverpool, Liverpool, United Kingdom

L. Fulton,
D. Hutchcroft0000-0002-4174-6509,
E. Maurice$\,$[33]0000-0002-7366-4364,
B. Wrona0000-0002-1555-0262

Physics and Astronomy Department, George Mason University, Fairfax, Virginia, USA

A. Conovaloff,
P. Cooper,
D. Coward$\,$[34]0000-0001-7588-1779,
P. Rubin0000-0001-6678-4985

Authors affiliated with an Institute or an international laboratory covered by a cooperation agreement with CERN

A. Baeva,
D. Baigarashev$\,$[35]0000-0001-6101-317X,
D. Emelyanov,
T. Enik0000-0002-2761-9730,
V. Falaleev$\,$[4]0000-0003-3150-2196,
S. Fedotov,
K. Gorshanov0000-0001-7912-5962,
E. Gushchin0000-0001-8857-1665,
V. Kekelidze0000-0001-8122-5065,
D. Kereibay,
S. Kholodenko$\,$[27]0000-0002-0260-6570,
A. Khotyantsev,
A. Korotkova,
Y. Kudenko0000-0003-3204-9426,
V. Kurochka,
V. Kurshetsov0000-0003-0174-7336,
L. Litov$\,$[14]0000-0002-8511-6883,
D. Madigozhin0000-0001-8524-3455,
M. Medvedeva,
A. Mefodev,
M. Misheva$\,$[36],
N. Molokanova,
S. Movchan,
V. Obraztsov0000-0002-0994-3641,
A. Okhotnikov0000-0003-1404-3522,
A. Ostankov$\,$[2],
I. Polenkevich,
Yu. Potrebenikov0000-0003-1437-4129,
A. Sadovskiy0000-0002-4448-6845,
V. Semenov$\,$[2],
S. Shkarovskiy,
V. Sugonyaev0000-0003-4449-9993,
O. Yushchenko0000-0003-4236-5115,
A. Zinchenko$\,$[2]

$^{\footnotemark[1]}$Corresponding author: E. Goudzovski, email<EMAIL_ADDRESS>
$^{1}$Present address: Ecole Polytechnique Fédérale Lausanne, CH-1015 Lausanne, Switzerland

$^{2}$Present address: Syracuse University, Syracuse, NY 13244, USA

$^{3}$Present address: Faculty of Science and Technology, University of Lancaster, Lancaster, LA1 4YW, UK

$^{4}$Present address: INFN, Sezione di Perugia, I-06100 Perugia, Italy

$^{5}$Also at TRIUMF, Vancouver, British Columbia, V6T 2A3, Canada

$^{6}$Present address: Université Catholique de Louvain, B-1348 Louvain-La-Neuve, Belgium

$^{7}$Present address: CERN, European Organization for Nuclear Research, CH-1211 Geneva 23, Switzerland

$^{8}$Present address: Institut für Kernphysik and Helmholtz Institute Mainz, Universität Mainz, Mainz, D-55099, Germany

$^{9}$Present address: Universität Würzburg, D-97070 Würzburg, Germany

$^{10}$Present address: European XFEL GmbH, D-22869 Schenefeld, Germany

$^{11}$Present address: School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK

$^{12}$Present address: Institut für Physik and PRISMA Cluster of Excellence, Universität Mainz, D-55099 Mainz, Germany

$^{13}$Also at Dipartimento di Scienze Fisiche, Informatiche e Matematiche, Università di Modena e Reggio Emilia, I-41125 Modena, Italy

$^{14}$Also at Faculty of Physics, University of Sofia, BG-1164 Sofia, Bulgaria

$^{15}$Present address: Scuola Superiore Meridionale e INFN, Sezione di Napoli, I-80138 Napoli, Italy

$^{16}$Present address: Instituto de Física, Universidad Autónoma de San Luis Potosí, 78240 San Luis Potosí, Mexico

$^{17}$Present address: Institut am Fachbereich Informatik und Mathematik, Goethe Universität, D-60323 Frankfurt am Main, Germany

$^{18}$Present address: INFN, Sezione di Roma I, I-00185 Roma, Italy

$^{19}$Also at Department of Industrial Engineering, University of Roma Tor Vergata, I-00173 Roma, Italy

$^{20}$Also at Department of Electronic Engineering, University of Roma Tor Vergata, I-00173 Roma, Italy

$^{21}$Also at Università degli Studi del Piemonte Orientale, I-13100 Vercelli, Italy

$^{22}$Also at Universidad de Guanajuato, 36000 Guanajuato, Mexico

$^{23}$Present address: Charles University, 116 36 Prague 1, Czech Republic

$^{24}$Present address: DESY, D-15738 Zeuthen, Germany

$^{25}$Present address: Max-Planck-Institut für Physik (Werner-Heisenberg-Institut), München, D-80805, Germany

$^{26}$Present address: Brookhaven National Laboratory, Upton, NY 11973, USA

$^{27}$Present address: INFN, Sezione di Pisa, I-56100 Pisa, Italy

$^{28}$Present address: Department of Physics, University of Warwick, Coventry, CV4 7AL, UK

$^{29}$Present address: Center for theoretical neuroscience, Columbia University, New York, NY 10027, USA

$^{30}$Present address: Dipartimento di Fisica dell'Università e INFN, Sezione di Genova, I-16146 Genova, Italy

$^{31}$Present address: School of Physical Sciences, University of Liverpool, Liverpool, L69 7ZE, UK

$^{32}$Present address: Dipartimento di Fisica e Astronomia dell'Università e INFN, Sezione di Firenze, I-50019 Sesto Fiorentino, Italy

$^{33}$Present address: Laboratoire Leprince Ringuet, F-91120 Palaiseau, France

$^{34}$Also at SLAC National Accelerator Laboratory, Stanford University, Menlo Park, CA 94025, USA

$^{35}$Also at L.N. Gumilyov Eurasian National University, 010000 Nur-Sultan, Kazakhstan

$^{36}$Present address: Institute of Nuclear Research and Nuclear Energy of Bulgarian Academy of Science (INRNE-BAS), BG-1784 Sofia, Bulgaria


# Quantifying high-order interdependencies on individual patterns via
the local O-information: theory and applications to music analysis

Tomas Scagliarini  Dipartimento Interateneo di Fisica, Universitá degli Studi
Aldo Moro, Bari and INFN, Italy  Daniele Marinazzo  Department of Data
Analysis, Ghent University, Belgium  Yike Guo  Department of Computer Science,
Hong Kong Baptist University, Hong Kong  Data Science Institute, Imperial
College London, United Kingdom  Sebastiano Stramaglia  Dipartimento
Interateneo di Fisica, Universitá degli Studi Aldo Moro, Bari and INFN, Italy
Center of Innovative Technologies for Signal Detection and Processing (TIRES),
Universitá degli Studi Aldo Moro, Italy  Fernando E. Rosas  Data Science
Institute, Imperial College London, United Kingdom  Centre for Psychedelic
Research, Department of Brain Science, Imperial College London, United Kingdom
Centre for Complexity Science, Imperial College London, United Kingdom

###### Abstract

High-order, beyond-pairwise interdependencies are at the core of biological,
economic, and social complex systems, and their adequate analysis is paramount
to understand, engineer, and control such systems. This paper presents a
framework to measure high-order interdependence that disentangles their effect
on each individual pattern exhibited by a multivariate system. The approach is
centred on the _local O-information_ , a new measure that assesses the balance
between synergistic and redundant interdependencies at each pattern. To
illustrate the potential of this framework, we present a detailed analysis of
music scores from J.S. Bach, which reveals how high-order interdependence is
deeply connected with highly non-trivial aspects of the musical discourse. Our
results place the local O-information as a promising tool of wide
applicability, which opens new perspectives for analysing high-order
relationships in the patterns exhibited by complex systems.

##  I Introduction

The analysis of interdependence is crucial for understanding the staggering
complexity of structures and behaviours manifested in biological, economic and
social systems. The unprecedented amount of data available for scientific
scrutiny provides unique opportunities to deepen our understanding of
multivariate co-evolving complex systems, including the orchestrated activity
of multiple brain areas, the interactions between different genes, and the
relationship between various econometric indices. Importantly, what allows
these systems to be more than the sum of their parts is not to be found in the
material nature of their parts, but in the fine structure of their
interdependencies [1].

Information theory provides an ideal framework to study interdependencies in
multivariate system, which establishes the notion of _information_ as a common
currency under which diverse systems can be measured and compared [2]. A
particularly promising approach for analysing the structure of
interdependencies is the partial information decomposition (PID), which
distinguishes different ‘modes’ of information that multiple predictors convey
about a target variable [3, 4, 5]. Two paradigmatic examples of such modes are
synergy and redundancy [6, 7, 8, 9, 10]: redundancy corresponds to information
which can be retrieved independently from more than one source, while synergy
correspond to statistical relationships that exist in the whole but cannot be
seen in the parts — this being rooted in the elementary fact that variables
can be pairwise independent while being globally correlated.

Despite continuous efforts to develop PID, the precise way in which synergies
and redundancies should be calculated is still being revised [11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24] . One way to circumvent this challenge
is to avoid computing the full decomposition, and study mixtures of PID modes
that can be captured by linear combinations of Shannon measures. One such
measure is the O-information [25], which has been shown to effectively capture
the overall balance between redundant and synergistic modes. The effectiveness
of the O-information in practical analyses has been verified by recent
applications on populations of spiking neurons [26], and the relationship
between neural patterns and ageing [27].

An important limitation of the O-information is that it characterises a
multivariate system with a single number, which summarises to the aggregated
effect of various patterns. Building on the rich literature of pointwise
information measures [28, 29, 30], in this paper we introduce the _local
O-information_ , which evaluates each pattern separately — such that its
ensemble average recovers the O-information. More specifically, the local
O-information constitutes an overall measure that characterise the high-order
interdependencies between the parts of a multivariate system at each possible
pattern of activity. Put differently, the local O-information evaluates the
‘statistical quality’ of each pattern, providing a signed scalar that assesses
the balance between redundancies an synergies at each individual pattern.

This paper presents the theory behind the local O-information, and then
illustrates its rich capabilities by analysing the scores of the chorales of
J.S. Bach. Our results show how the local O-information is capable of
revealing subtle musical relationships, including properties of different
intervals, chord dispositions, harmonic depth, and the relationship between
music and text. Thanks to its ability to uncover such highly non-trivial
relationships, the local O-information is a valuable addition to the toolkit
of data analysts interested in the study of complex systems.

The rest of the paper is organised as follows. Section II provides background
information about the O-information, introduces the new local O-information,
and then illustrates its basic properties on small spin systems. Then, Section
III presents a detailed analysis of the local O-information on the chorales of
J.S Bach, and finally Section VI summarises our main conclusions.

##  II A local measure of information quality

Let us consider a scientist who is interested in studying a given complex
system, whose state can be appropriately described by the vector
$\bm{X}^{n}=(X_{1},\dots,X_{n})$. We focus on scenarios where the scientist
has enough data to build a reliable statistical description of the statistics
of $\bm{X}^{n}$, which is denoted by $p(\bm{X}^{n})$. A question of interest
is how to leverage the statistics encoded in $p$ in order to deepen our
understanding of the structure of interdependencies that characterize
$\bm{X}^{n}$. Such understanding can lead either to the building of
statistical markers to classify different systems or different states of the
same system, or to establish parallels between seemingly heterogeneous systems
based on the similarity of their relational structure.

Through this section, random variables are denoted by capital letters (e.g.
$X,Y$) and their realisations by lower case letters (e.g. $x,y$). Random
vectors and their realisations are denoted by capital and lower case boldface
letters, respectively.

###  II.1 O-Information

Shannon’s mutual information is a popular metric of interdependency, which
overcomes the limitations of correlation metrics such as Pearson’s in that it
captures both linear and non-linear relationships and being applicable to
ordinal data. However, the mutual information can only asses the relationships
between two (sets of) variables, being unable to fully explore the rich
interplay that can take place within triple or higher-order interactions.

Two multivariate extensions of the mutual information are the Total
Correlation (TC) [31] and the Dual Total Correlation (DTC) [32], which are
defined as

$\displaystyle\text{TC}(\bm{X}^{n})$
$\displaystyle:=\sum_{i=1}^{n}H(X_{i})-H(\bm{X}^{n}),$
$\displaystyle\text{DTC}(\bm{X}^{n})$
$\displaystyle:=H(\bm{X}^{n})-\sum_{i=0}^{N}H(X_{i}\mid\bm{X}^{n}_{-i}).$

Above, $H(X_{i})=-\sum_{x_{i}}p(x_{i})\log p(x_{i})$ corresponds to the
Shannon entropy, $H(X_{i}|X_{j})=H(X_{i},X_{j})-H(X_{j})$ is the conditional
Shannon entropy, and $\bm{X}^{n}_{-i}$ is the vector of all variables except
$X_{i}$ (i.e.,
$\bm{X}^{n}_{-i}=\left(X_{1},\ldots,X_{i-1},X_{i+1},\ldots,X_{n}\right)$);
hence, the term $H(X_{i}\mid\bm{X}^{n}_{-i})$ quantifies how $X_{i}$ is
independent from the other $n-1$ variables. As the mutual information, both TC
and DTC are non-negative quantities which are zero if and only if all
variables $X_{0},\dots,X_{N}$ are jointly statistically independent — i.e. if
$p(\bm{X}^{n})=\prod_{i=1}^{n}p(X_{i})$.

Despite the similarities between TC and DTC, these two metrics provide
distinct but complementary views on the strength of the interdependencies in a
multivariate system. On the one hand, the TC accounts for _collective
constraints_ , which refer to regions of the phase space that the system
explore less; while on the other hand, the DTC measures the degree of _shared
randomness_ between the variables, i.e. the amount of information that can be
collected in one variable that also refers to the activity of another [25]. An
attractive way to exploit these complementary views is by considering their
difference,

$\displaystyle\Omega_{n}(\bm{X}^{n})$
$\displaystyle=\text{TC}(\bm{X}^{n})-\text{DTC}(\bm{X}^{n})~{}.$ (1)

which is known as the _O-information_ [25]. The O-information can be seen as a
revision of the measure of neural complexity proposed by Tononi, Sporns and
Edelman in [33], which provides a mathematical construction that is closer to
their original desiderata 111For a comparison between the original measure
proposed in Ref. [33] and the O-information, please see Ref. [25].. In effect,
the O-information is a signed metric that captures the balance between high-
and low-order statistical constraints 222While low-order constraints impose
strong restrictions on the system and allow little amount of shared
information between variables, high-order constraints impose collective
restrictions that enable large amounts of shared randomness.. By construction,
$\Omega(\bm{X}^{n})<0$ implies a predominance of high-order constraints within
the system $\bm{X}^{n}$, a condition that is usually referred to as
statistical synergy. Conversely, $\Omega(\bm{X}^{n})>0$ implies that the
system $\bm{X}^{n}$ is dominated by low-order constraints, which imply
redundancy of information. This nomenclature is further supported by the
following key properties:

  * (1)

It captures genuine high-order effects, at it is zero for systems with only
pairwise interdependencies: if the joint distribution of $\bm{X}^{n-1}$ (for
$n$ odd) can be factorised as
$p_{\bm{X}^{n}}(\bm{x}^{n})=\prod_{k=1}^{n/2}p_{X_{2k},X_{2k+1}}(x_{2k},x_{2k+1})$,
then $\Omega(\bm{X}^{n})=0$.

  * (2)

The O-information is maximised by redundant distributions where the same
information is copied in multiple variables, and is minimised by synergistic
(“xor-like”) distributions: e.g. for binary variables, $\Omega$ is maximised
by the “n-bit copy” where $X_{1}$ is a Bernoulli r.v. with parameter $p=1/2$
and $X_{0}=\dots=X_{N-1}$, and is minimised when $X_{0},\dots,X_{N-2}$ are
i.i.d. fair coins and $X_{N-1}=\sum_{j=0}^{N-2}X_{j}\ (\mathrm{mod}\ 2)$.

  * (3)

The O-information characterises the dominant tendency, being additive over
non-interactive subsystems: if the system can be factorised as $p_{\bf
X^{N}}(\bm{x}^{N})=p_{X_{1},\dots,X_{m}}(x_{1},\dots,x_{m})\times
p_{X_{m+1},\dots,X_{N}}(x_{m+1},\dots,x_{N})$, then
$\Omega(\bm{X}^{n})=\Omega(X_{1},\dots,X_{m-1})+\Omega(X_{m},\dots,X_{n})$.

For more details related to the O-information, we refer the reader to Ref.
[25].

###  II.2 Local O-Information

Building on the properties of the O-information reviewed in the previous
section, one can design a measure that can capture these effects on a state-
by-state basis. In effect, the O-information provides a single scalar that
characterises the interdependencies of a system _on average_. However, in many
systems of interest this average represents a mid-point between highly
heterogeneous states, which could potentially be of limited value to
understand the role of individual patterns. In this section we introduce a
_pointwise_ measure that allows to calculate the O-information on individual
states.

As a first step, let us introduce a _local total correlation_ and _local dual
total correlation_ which are given by

$\displaystyle\text{tc}(\bm{x}^{n}):=$
$\displaystyle\sum_{j=1}^{n}h(x_{j})-h(\bm{x}^{n}),$ (2)
$\displaystyle\text{dtc}(\bm{x}^{n}):=$
$\displaystyle\>h(\bm{x}^{n})-\sum_{j=1}^{n}h(x_{j}|\bm{x}^{n}_{-j}),$ (3)

where $h(\bm{x}^{n})=-\log p(\bm{x}^{n})$ is the information content of the
state $\bm{x}^{n}$ [28]. These quantities capture how the strength of the
multivariate interdependencies vary with the state of the system, providing a
generalisation to the pointwise mutual information introduced in Ref [36].

Using these definitions, we can define the _local O-information_ as follows:

$\displaystyle\omega(\bm{x}^{n}):=$
$\displaystyle\;\text{tc}(\bm{x}^{n})-\text{dtc}(\bm{x}^{n})$ (4)
$\displaystyle=$
$\displaystyle\;(n-2)h(\bm{x}^{n})+\sum_{j=1}^{n}\Big{(}h(x_{j})-h(\bm{x}_{-j}^{n})\Big{)}.$
(5)

In contrast to $\Omega(\bm{X}^{n})$, which provides a single value for the
random variable, $\omega(\bm{x}^{n})$ assigns a number of each possible state
$\bm{x}^{n}$. In particular, the local O-information has the following useful
properties:

  * $\bullet$

$\Omega(\bm{X}^{n})=\mathbb{E}\\{\omega(\bm{x})\\}$.

  * $\bullet$

$\inf_{\bm{x}}\omega(\bm{x})\leq\Omega(\bm{X}^{n})\leq\sup_{\bm{x}}\omega(\bm{x})$.

Generally speaking, the local O-information provides more information than its
global counterpart as, technically, $\omega(\bm{X}^{n})$ is a random variable
whose mean value is $\Omega(\bm{X}^{n})$. Therefore, the whole range of values
of $\omega(\bm{X}^{n})$ can naturally provide a more fine-grained description
of the system than its mere first moment.

Additionally, note that Lemma 3 of Ref. [25] provides upper and lower bounds
for $\Omega(\bf{X})$ when the variables take values in a finite alphabet; in
particular, if $X_{j}\in\mathscr{X}$ for all $j$ then

$-(n-2)\log|\mathscr{X}|\leq\Omega(\bm{X}^{n})\leq(n-2)\log|\mathscr{X}|~{}.$
(6)

Note that these bounds don’t apply to the local O-information — they apply to
its average value, but extreme values can be larger than it. Nonetheless, the
quantity $(n-2)\log|\mathscr{X}|$ establishes a natural bound which that, when
surpassed, values of $\omega$ can be considered to be particularly large. This
provides a useful rule of thumb to interpret the ranges of values obtained by
evaluations of $\omega$.

###  II.3 Proof of concept

An useful way of employing $\omega$ is to classify states among different
_types_. In particular, building on the properties of the O-information, we
say that a state $\bm{x}^{n}$ for which $\omega(\bm{x}^{n})>0$ is redundancy-
dominated, while if $\omega(\bm{x}^{n})<0$ we say the state is synergy-
dominated. In this section we illustrate this capability of the local
O-information by using it to analyse a small Ising system.

Let’s consider three coupled spins denoted by $(S_{1},S_{2},S_{3})$, whose
joint probability distributions follows a Boltzmann-Gibbs distribution of the
form:

$p(s_{1},s_{2},s_{3})=\frac{e^{J(s_{1}s_{2}+s_{1}s_{3}+s_{2}s_{3})}}{Z},$ (7)

with $Z=\sum_{s1,s2,s3}e^{J(s_{1}s_{2}+s_{1}s_{3}+s_{2}s_{3})}$ is a
normalisation factor. For positive values of $J$, the configurations with all
the spins in agreement (i.e. $\uparrow\uparrow\uparrow$ and
$\downarrow\downarrow\downarrow$) satisfy all bonds, and hence the system may
be seen as a small ferromagnet. In contrast, for negative $J$ the system is
said to be “frustrated” as there is not a configuration satisfying all bonds
simultaneously.

Figure 1: Concerning the Ising toy model of three spins, $\Omega$ is plotted
versus the coupling $J$ (top panel). The values of $\omega$ for the $8$
possible configurations of the three spins are shown for the unfrustrated case
when $J=1$ (middle panel) and for the frustrated case, when $J=-1$ (bottom
panel).

By analysing Eq. (7) via the O-information, one finds that ferromagnetic
behaviour are redundancy-dominated while frustrated systems are synergy-
dominated (see Figure 1). Interestingly, the local O-information shows that
configurations of spin agreement are redundancy-dominated states, whilst the
six configurations with disagreement are synergy-dominated — this for both
positive and negative values of $J$. This let us conclude that what makes the
system redundancy- or synergy-dominated for different values of $J$ is the
different frequency with which either redundancy- or synergy-dominated
configurations are visited.

This toy example shows how the local O-information can reveal different
qualities of various states — in this case, either states with agreement or
disagreement. Furthermore, this example also reveals an intriguing connection
between synergy and frustration in spin models, which will be further
investigated in a future publication.

##  III Case study: high-order relationships in Bach’s chorales

To illustrate the usefulness of the local O-information for practical data
analysis, this section presents a study of the multivariate statistics of
musical scores from the Baroque period. Ref. [25] provided an analysis of
music scores based on the global O-information; however, such analyses could
not provide information about individual chords, and hence could not explore
further musical aspects related to harmony and tonality. Here we show how the
local O-information can greatly expand this type of analyses, revealing subtle
aspects of the musical discourse that are reflected in the high-order
interactions.

In the following, Section III.1 describes the procedure to obtain and analyse
the data, and Section III.2 discusses our main findings.

###  III.1 Processing pipeline

####  III.1.1 Data

Our analysis focuses on the chorales for four voices (soprano, alto, tenor,
and bass) written by Johann Sebastian Bach (1685-1750). These works are
characterised by an elaborate counterpoint between the melodic lines that
leads to rich harmonic progressions, which in turn results into a broad range
of chords displayed along the repertoire. An additional point of interest of
these pieces is that, as typical in the Baroque period (approx. 1600–1750),
they display a balance in the interest and richness of each of the four
voices. This contrasts with the subsequent Classic (1730–1820) and Romantic
(1780–1910) periods, where higher voices tend to take the lead while the lower
voices provide mere support.

Figure 2: Processing pipeline Top: each chorale was selected among those
written in major mode and then transposed in C major. Middle: each melodic
line was transformed into a time series of $13$ possible values. Bottom: for
each 4-note state the local O-information was calculated using Eq. (5).

Our analysis is based on the electronic scores publicly available at
http://kern.ccarh.org, a website that hosts professionally curated digital
scores [37]. Our pre-processing pipeline is the same as in Ref. [25], which we
describe here for completeness. The scores were pre-processed in Python using
the Music21 package (http://web.mit.edu/music21), which allowed us to select
only the pieces written in Major mode. Each chorale was transposed to C Major,
and each melodic line was transformed into time series of 13 possible values
(one for each note plus one for the silence), using a small rhythmic duration
as common time unit. This resulted in a total of 172 chorales, which gave
$\approx 4\times 10^{4}$ four-note chords.

With this data, the joint distribution of the values for the four-note chords
was estimated using their empirical frequency 333Please note that
regularisation methods — such as Laplace smoothing — can have significant
effects on the results. We decided not to use such methods, as some chords
(e.g. C-C$\sharp$-D-D$\sharp$) are not representative of the Baroque
repertoire.. This leads to a probability assigned to each four-note chord,
which assess simply the odds of picking that chord when randomly selecting one
out of the whole repertoire. One can express this probability as the
multivariate statistic $p(x_{1},x_{2},x_{3},x_{4})$, with each variable
corresponding to the different voices.

Finally, we used $p(x_{1},x_{2},x_{3},x_{4})$ to calculate the local
O-information $\omega(\bm{x})$ for each chord $\bm{x}$ using Eq. (5), which
determines the dominant statistical behaviour (in terms of synergy and
redundancy) associated with each chord. The overall pipeline starting from the
music score and arriving to the local O-information is illustrated in Figure
2.

####  III.1.2 Research questions and tools

We studied the multivariate properties of each of the possible four-note
chords of Bach’s chorales. Our analysis focuses exclusively on harmony and
chords, leaving melodic and rhythmic properties to future studies. We focus on
the question of what harmonic properties of the music tend to give rise to
synergistic or redundant high-order relationships between the four voices.

Let us denote by $\bm{X}=(X_{1},X_{2},X_{3},X_{4})$ the random vector that
follow the statistics encapsulated by $p(x_{1},x_{2},x_{3},x_{4})$. Following
standard musical practice, we follow the convention that the variables go from
lower to higher range, so that $X_{1}$ corresponds to the bass and $X_{4}$ is
the soprano. Moreover, we use the shorthand notation CEGE when referring to
the chord
$(x_{1},x_{2},x_{3},x_{4})=(\texttt{C},\texttt{E},\texttt{G},\texttt{E})$.

Note that $\bm{X}$ can adopt $13^{4}=28561$ possible values, and that $p$ is
generally not invariant under changes of ordering between the four voices.
Since $X_{1},\dots,X_{4}$ take values among alphabets of cardinality
$|\mathscr{X}|=13$, we do all calculations using logarithms to base $13$, so
that $H(X_{k})\leq 1$ for all $k\in\\{1,\dots,4\\}$ — we call this unit a
_mut_ , for _musical bit_. Eq. (6) implies that $-2<\Omega(\bm{X})<2$, and
hence most values of $\omega$ are expected to have absolute value less than 2
muts — which gives a sense of how to interpret the magnitude of local
O-information values.

We expected to find a correspondence between tonality and O-information
values. In particular, we hypothesise that the principal tonal chords (C, F,
and G major) would be associated with redundant behavior, while chords that
are farther away from the tonal centre (i.e. involve many sharp or flat
alterations) would be related to synergistic events. Additionally, we expect
dissonance to be associated with less redundancy, as it involves more complex
combinations of notes.

###  III.2 Results

####  III.2.1 Analysis of the extreme values of the local O-information

Out of the $13^{4}$ possible chords, we found that only $1715$ of them were
observed at least once in the chorales, corresponding to only $6\%$ of the
possibilities — reflecting the specificity of the chord choices used in Bach’s
chorales. A weak correlation is observed between frequency and local
O-information: more frequent chords tend to have a higher $\omega$ — which
suggests that more visited chords tend to be made by more redundant parts (see
Figure 3). The most frequently encountered chords are shown in the Appendix
(Table A).

Figure 3: Scatter plot of the logarithm of the probability of each state
versus $\omega$. Almost all the states have $\omega$ between $-2$ and $+2$ —
i.e. the bounds for $\Omega$. The outlier in redudancy corresponds to the
chord RRRR.

Some interesting observations can be made by observing the most positive
(redundant) and the most negative (synergistic) states in terms of $\omega$,
which are presented in Table III.2.1. On the one hand, the most redundant
states tend to contain few alterations (sharp notes, denoted in the table with
the symbol $\sharp$) and mostly consonant intervals 444For a description of
consonant and dissonant intervals, please see next section.. In contrast,
synergistic chords tend to contain more alterations and dissonant intervals,
which in the Western culture are typically associated with harshness and
unpleasantness. For example, the most synergistic chord contains a major
seconds (D-E), while the second most synergistic has one minor second
(F$\sharp$-G) and one major second (E-F$\sharp$). The “chord” with highest
local O-information is found to be RRRR, where the redundancy can be
interpreted as a consequence of the voices doing the same thing — not signing.

[/csv/respect sharp=true, filter=¡20, table head= Redundancy Synergy
Chord $\omega$ Chord $\omega$
]Table1.csv

Table 1: Chords with the highest (redundance) and lowest (synergy) local
O-information. Letters refers to the standard music nomenclature (plus R is
for silence), and the ordering of the voices is Bass-Tenor-Alto-Soprano — from
left to right.

####  III.2.2 The role of intervals

The results reported in the previous section imply a link between the musical
(harmonic) properties of a chord and the type of the statistical
interdependencies within it constituent notes. In particular, results suggest
that synergistic interdependencies may be linked to the presence of
dissonances, while redundancy may be related to consonance. A consonant
interval occurs when the ratio of the frequencies between two notes is very
simple, like $(1\\!:\\!2)$ for the octave, $(2\\!:\\!3)$ for the perfect fifth
or $(4\\!:\\!5)$ for the major third. In the Western culture, consonance is
typically associated by listeners with pleasantness and acceptability [40]. In
Westerner music theory, dissonant intervals typically include the major second
($8\\!:\\!9$) and minor second ($15\\!:\\!16$), major seventh ($8\\!:\\!15$)
and minor seventh ($9\\!:\\!16$), and the augmented forth (so-called ‘tritone’
or diabolus in musica).

To further explore the relationship between high-order statistics and harmony,
we studied how the local O-information depends on the number of dissonant
intervals (either seconds/sevenths or augmented fourths) a chord possesses.
The results are depicted in Figure 4. Analysis of variance with Bonferroni
correction revealed a significant dependency between number of dissonances and
$\omega$ for all number of chords — except for the contrast between 3 and 4
dissonances, arguably due to the small number of chords with 4 dissonances.
The results of all the comparison are shown in Table III.2.2.

[ table head=Dissonances $p-$value Cohen’s $d$
]Table2.csv The chords that contain 3 or more dissonant intervals are
presented in the Appendix (Table A), most of which exhibits negative values of
$\omega$. Remarkably, the only two chords of those with $\omega>0$ can be
identified as part of G major with added 7-th, which is the most frequent
dissonant chord in classical harmony.

Table 2:

Figure 4:  Dissonance vs local O-information. Each box represents a state with
a different number of dissonant interval. Each category is statistically
different from the one with no dissonance (pure consonant intervals).

A question that raises from the results shown in Figure 4 is why purely
consonant chords can be dominantly synergistic, as shown by the variance of
the values of $\omega$ for zero dissonance. For example, the chord
$\omega(\texttt{CEGC})=0.42$ is redundant while $\omega(\texttt{EGCC})=-0.45$,
being both C major chords but having a different pitch in the bass (the
fundamental note in the first case, and the third in the second). Leveraging
music theory, an explanation from this comes from the notion of ‘chord
inversion:’ a triad chord is in _first inversion_ if the third (either major
or minor) is in the bass, it is in _second inversion_ if the fifth is in the
bass, and it is in _root position_ if the first/fundamental note is in the
bass. In Western classical music each inversion tends to be associated with
specific sensations – the first inversion gives a sense of lightness, while
the second inversion and root position are typically associated with
instability and stability, respectively. By considering the values of $\omega$
corresponding to different inversions, t-test shows a tendency towards lower
values of $\omega$ in chords in first (Cohen’s $d\simeq 0.36$) and second
(Cohen’s $d\simeq 0.31$) inversion when compared to chords in root position,
as shown in Figure 5.

Figure 5: The effect of chords inversions on the local O-information. Root
form occurs when the root note of the triad is in the bottom (e.g. CEG), first
inversion when the third of the triad is the bass note (e.g. EGC), and second
inversion when the fifth is in the bass (e.g. GCE).

Figure 6: Networks of redundant (orange) and synergy-dominated (blue)
relationships in Bach’s chorales. For each voice pair the links were created
by averaging the local O-information of all chords that contain that pair of
notes.

Finally, as a complementary way to study the role of intervals on the
O-information, we considered the average value of $\omega$ for given notes at
specific voices — averaging over all possible notes adopted by the other two
voices. Results are shown in Figure 6. It was found that redundancy (i.e. the
most positive values of O-information) is ‘localised’ in few intervals, taking
place mainly between notes involving the tonic (C major, CEG) or dominant (G
major, GBD) chords, or between silences. Also, most redundancy in the bass is
associated to the fundamental note of each chord — either C or G. In contrast,
synergy (i.e. the most negative values of O-information) are much more
widespread. Interestingly, the redundancies between the two extreme voices
(soprano and bass) are relatively weak (except between their silence), while
synergies between them are not. Please note that the extreme voices tend to
carry an important role in Bach chorales — the soprano carrying out the main
melody, and the bass leading the harmony.

####  III.2.3 Harmonic depth

While the previous section focused on the role of single intervals, now our
analyses focuses on harmonic considerations. Harmony is organised around a
_tonality_ (also called ‘key’), which plays the role of centre of gravity
around which music discourse revolves. The axial pitch of a tonality is called
the _root_ , which in turn gives name to the tonality — e.g. C is the root of
the tonality of C major. In classical Western music there are 12 different
major tonalities, one for each of each pitch. Also, each major tonality has an
associated minor tonality, which is located a minor third below (e.g. C major
is associated with A minor). Each of these 12 major tonalities are made of 7
distinct pitches, and are naturally ordered by a notion of proximity depending
on how many pitches do they have in common. This gives raises to the _circle
of fifths_ : major keys separated by a fifth have only one note different. For
example, C and G major are only distinguished by the note F, which is sharp
for the latter but natural for the former.

A simple way to explore the impact of harmony on the high-order statistics is
by analysing the dependency between $\omega$ and the number of alterations
(sharps or flats) that a chord has. In effect, please recall that all the
chorales analysed are in major mode, and have been shifted to C (see Section
III.1). Moreover, chords belonging to C major have no alterations, while
chords for more distant tonalities have progressively more alterations —
either sharps if going up the cycle of fifths, or flats otherwise. Therefore,
we ran statistical analyses (t-test corrected for multiple comparisons) on the
effect of the number of alterations on $\omega$, whose results are shown in
Figure 7. Results revealed significantly decreases of $\omega$ ($p<0.01$) for
states with one (Cohen’s $d\simeq 0.41$) or two (Cohen’s $d\simeq 0.44$)
alterations with respect to chords without any alterations.

Figure 7: The local O-information is related to the number of alterations
(sharps and flat notes) inside each state. States with 1 and 2 alterations
have $\omega$ significantly lower ($p<0.05$) than states with no alterations.
Figure 8: The effect of harmonic distance. The figure shows the mean values
(and confidence intervals) of $\omega$ for musical states in every different
major (top) and minor (bottom) chord. The blue line indicates purely consonant
triads, while the orange line corresponds to chords with one or more
dissonances.

In order to go deeper into the impact of harmony on the O-information we
introduced the notion of _harmonic depth_ , which corresponds to the smallest
number of steps (going clockwise or anti-clockwise) in the circle of fifths
are required to go from $C$ major to a given tonality — or from $A$ minor, in
case of minor tonality. For example, the chord $D$ major has an harmonic depth
of $+2$, while the chord F minor (denoted as Fm) has an harmonic depth of
$+4$.

We are interested to study the relationship between $\omega$ and harmonic
depth. For this purpose, we consider the values of $\omega$ for triads that
belong to a specific tonality, regardless to the arrangement of notes between
the voices For example, determining a chord is C just accounts for triads
containing only the pitches C,E,G, regardless how they are arranges among the
voices. Results are shown in Figure 8, and show that the value of $\omega$
decreases as soon as the tonality moves away from C major, being this
difference more pronounced in minor chords. This suggest that synergy may also
be associated with more complex harmonic explorations involving more
harmonically distant chords.

####  III.2.4 Lyric analysis

As a final step in our analysis, we investigated the relationship between the
values of $\omega$ due to the chords and the corresponding word that is sang
as part of the lyrics. For this purpose, we consider the different values of
$\omega$ that correspond to each time a given word is sang. In cases of
melismas (i.e. when many notes are sang under the same syllable), the values
of the whole progression were averaged and counted as one realisation of the
word.

As a first analysis, we calculated a word cloud where words associated with
negative or positive values of $\omega$ are represented in red and blue,
respectively. As shown in Figure 9, many of the most common words (like Gott,
Herr, Sohn) are redundant, with the exception of Jesu that is synergistic. As
the majority of the chords explored by Bach are synergistic (the average
O-information is negative, see Figure 3), this prevalence of redundant words
is highly non-trivial.

Figure 9: Word cloud of Bach’s chorales lyrics. Most common words found in
text, with their size representing their frequency and their colour their
local O-information sign (plus or minus).

Another insight that can be drawn from the word cloud is that words that are
not the subject of the phrase seems to be more synergistic. To verify this, we
evaluated the effect on $\omega$ of words being in root form (nominative case)
with respect to all others 555German language has four cases: nominative
(subject), accusative (direct object), dative (indirect object), and genitive
(possessive).. Results confirmed our conjecture, showing that words in root
form have a tendency towards higher values of $\omega$ (see Figure 10). As an
speculation, this may be interpreted by noting that root form words correspond
to the most important part of the sentence, and hence a redundant underlying
harmony might contribute to a easier comprehension. The most frequent cases of
words for which we found both the root and the non root word are shown in
Table III.2.4.

[/csv/respect sharp=true, table head=Root $\omega$ No Root $\omega$
]Table3.csv

Table 3: Most common words and their mean local O-information value.

Figure 10: Root vs non root words. Each box shows the distribution of local
O-information $\omega$ for each category. A statistical comparison between the
two populations has been carried out with a two sample $t$-test, which
rejected the null hypothesis with $p<0.01$; the effect size can be expressed
as Cohen’s $d\simeq 0.26$.

##  IV Conclusions

This paper introduces a new framework to study the high-order
interdependencies observed in complex multivariate systems, which is capable
of disentangling their effects on individual patterns of activity. The
approach is centred on the local O-information, a measure that quantifies the
balance between redundancy and synergy at each pattern. Because of its
information-theoretic nature, this measure is widely applicable — being
suitable to assess systems with categorical, discrete, and continuous
variables.

The capabilities of the proposed framework were showcased in an analysis of
the scores of the chorales of J.S. Bach, which illuminated the high-order
relationships that exist between the different voices. In particular, our
results synergy-dominated interdependencies tend to be associated with complex
musical elements, including dissonances, chord inversions, and harmonic
distance from the tonal centre. Taken together, our findings provide
converging evidence about the relationship between statistical synergy and the
complexity of the musical discourse.

These findings have interesting parallels with recent studies on the human
brain, which are revealing a close relationship between synergistic
interdependencies in neural activity and high cognitive functions.
Historically, the notion of synergistic information was originated in
theoretical neuroscience as an effort to characterise aspects of complex
neural activity [33, 42, 6, 7, 43, 44]. Moreover, recent empirical work has
shown that synergy characterises the interactions between multimodal,
association brain areas that have undergone more evolutionary expansion [45],
and is also associated with changes in brain function due to anaesthesia and
disorders of consciousness [46]. Under the light of these findings, the
results presented in this work let us speculate that elaborated musical
discourse may have key similarities with the type of neural activity that
underpins high brain functions. One possible commonality could be presence of
emergent phenomena, which have been recently characterised formally in terms
of statistical synergy [47]. This would not be the first time music is shown
to share some of the hallmark properties of complex systems; in effect,
properties of musical discourse has been show to be related to non-linear
fluctuations and self-organised criticallity [48, 49], and also to entropy
production and irreversibility [50].

Pointwise information measures, initially proposed in Ref. [29] w.r.t. the
local transfer entropy, are a promising set of techniques that can provide a
detailed description of information transfer mechanisms in complex systems.
Recently this paradigm has been applied to implement the local Granger
causality [30], which has shown interesting results on physiological and
neural data. The fine descriptions allowed by the formalism introduced in this
paper bring a new perspective over high-order interdependencies, which
complement existent pointwise information decomposition approaches (e.g. Ref.
[51]) by being applicable to larger systems — hence greatly extending their
domain of practical applicability. The extension of these ideas to dynamical
information decomposition, such as the integrated information decomposition
framework [52], constitutes a promising direction for future research.

## Acknowledgements

The authors thank Pablo Padilla and Alejandro Reyes for insightful
discussions. T.S. wants to thank M. Paolo Daniele for having transmitted to
him a bit of his passion for music. F.R. was supported by the Ad Astra
Chandaria foundation. S.S. was supported by MIUR project PRIN 2017WZFTZP
“Stochastic forecasting in complex systems”.

## Code availability

Code for the analysis reported here is publicly available at
https://github.com/tomscag/local_O_information

##  Appendix A Tables

[/csv/respect sharp=true, table head=Chord $\omega$ Frequency
]Table4.csv

Table 4: Most common chords.

[/csv/respect sharp=true, table head= Redundancy Synergy
Chord $\omega$ Chord $\omega$
]Table1.csv

Table 5: Table shows the first states with the highest (redundance) and lowest
(synergy) $\omega$ values. The letters refers to the standard nomenclature of
notes in music (R is for rest) and the order of the voices is Bass-Tenor-Alto-
Soprano, from left to right.

[/csv/respect sharp=true, table head=Chord $\omega$ Occurrence
]Table5.csv  [/csv/respect sharp=true, table head=Chord $\omega$ Occurrence
]Table6.csv  Table 7: Most common dissonant chords (only chords with 3 or more
dissonant intervals are shown).

Table 6: Most common tonal chords.

## References

  * [1] James P. Crutchfield.  The calculi of emergence: computation, dynamics and induction.  Physica D: Nonlinear Phenomena, 75(1-3):11–54, 1994.  Publisher: Elsevier.
  * [2] James P. Crutchfield and David P. Feldman.  Regularities unseen, randomness observed: Levels of entropy convergence.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 13(1):25–54, 2003.  Publisher: American Institute of Physics.
  * [3] Paul L. Williams and Randall D. Beer.  Nonnegative decomposition of multivariate information.  arXiv preprint arXiv:1004.2515, 2010.
  * [4] Virgil Griffith.  Quantifying synergistic information.  California Institute of Technology, 2014.
  * [5] Joseph T. Lizier, Nils Bertschinger, Jürgen Jost, and Michael Wibral.  Information decomposition of target effects from multi-source interactions: Perspectives on previous, current and future work, 2018.
  * [6] Itay Gat and Naftali Tishby.  Synergy and redundancy among brain cells of behaving monkeys.  Advances in neural information processing systems, pages 111–117, 1999.  Publisher: MIT; 1998.
  * [7] Elad Schneidman, William Bialek, and Michael J. Berry.  Synergy, redundancy, and independence in population codes.  Journal of Neuroscience, 23(37):11539–11553, 2003.  Publisher: Soc Neuroscience.
  * [8] Nicholas Timme, Wesley Alford, Benjamin Flecker, and John M. Beggs.  Synergy, redundancy, and multivariate information measures: an experimentalist’s perspective.  Journal of computational neuroscience, 36(2):119–140, 2014.  Publisher: Springer.
  * [9] Fernando Rosas, Vasilis Ntranos, Christopher J. Ellison, Sofie Pollin, and Marian Verhelst.  Understanding interdependency through complex information sharing.  Entropy, 18(2):38, 2016.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [10] Michael Wibral, Viola Priesemann, Jim W. Kay, Joseph T. Lizier, and William A. Phillips.  Partial information decomposition as a unified approach to the specification of neural goal functions.  Brain and cognition, 112:25–38, 2017.  Publisher: Elsevier.
  * [11] Malte Harder, Christoph Salge, and Daniel Polani.  Bivariate measure of redundant information.  Physical Review E, 87(1):012130, 2013.  Publisher: APS.
  * [12] Virgil Griffith, Edwin KP Chong, Ryan G. James, Christopher J. Ellison, and James P. Crutchfield.  Intersection information based on common randomness.  Entropy, 16(4):1985–2000, 2014.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [13] Adam B. Barrett.  Exploration of synergistic and redundant information sharing in static and dynamical Gaussian systems.  Physical Review E, 91(5):052802, 2015.  Publisher: APS.
  * [14] Eckehard Olbrich, Nils Bertschinger, and Johannes Rauh.  Information decomposition and synergy.  Entropy, 17(5):3501–3517, 2015.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [15] Robin AA Ince.  The Partial Entropy Decomposition: Decomposing multivariate entropy and mutual information via pointwise common surprisal.  arXiv preprint arXiv:1702.01591, 2017.
  * [16] Conor Finn and Joseph T. Lizier.  Pointwise partial information decomposition using the specificity and ambiguity lattices.  Entropy, 20(4):297, 2018.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [17] Ryan G. James, Jeffrey Emenheiser, and James P. Crutchfield.  Unique information via dependency constraints.  Journal of Physics A: Mathematical and Theoretical, 52(1):014002, 2018.  Publisher: IOP Publishing.
  * [18] Ryan G. James, Jeffrey Emenheiser, and James P. Crutchfield.  Unique information and secret key agreement.  Entropy, 21(1):12, 2019.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [19] Nihat Ay, Daniel Polani, and Nathaniel Virgo.  Information decomposition based on cooperative game theory.  arXiv preprint arXiv:1910.05979, 2019.
  * [20] Fernando E. Rosas, Pedro AM Mediano, Borzoo Rassouli, and Adam B. Barrett.  An operational information decomposition via synergistic disclosure.  Journal of Physics A: Mathematical and Theoretical, 53(48):485001, 2020.  Publisher: IOP Publishing.
  * [21] Conor Finn and Joseph T. Lizier.  Generalised measures of multivariate information content.  Entropy, 22(2):216, 2020.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [22] Aaron J. Gutknecht, Michael Wibral, and Abdullah Makkeh.  Bits and pieces: Understanding information decomposition from part-whole relationships and formal logic.  Proceedings of the Royal Society A, 477(2251):20210110, 2021.  Publisher: The Royal Society Publishing.
  * [23] Kyle Schick-Poland, Abdullah Makkeh, Aaron J. Gutknecht, Patricia Wollstadt, Anja Sturm, and Michael Wibral.  A partial information decomposition for discrete and continuous variables.  arXiv preprint arXiv:2106.12393, 2021.
  * [24] Abdullah Makkeh, Aaron J. Gutknecht, and Michael Wibral.  Introducing a differentiable measure of pointwise shared information.  Physical Review E, 103(3):032149, 2021.  Publisher: APS.
  * [25] Fernando E. Rosas, Pedro AM Mediano, Michael Gastpar, and Henrik J. Jensen.  Quantifying high-order interdependencies via multivariate extensions of the mutual information.  Physical Review E, 100(3):032305, 2019.  Publisher: APS.
  * [26] Sebastiano Stramaglia, Tomas Scagliarini, Bryan C. Daniels, and Daniele Marinazzo.  Quantifying dynamical high-order interdependencies from the o-information: an application to neural spiking dynamics.  Frontiers in Physiology, 11:1784, 2021.  Publisher: Frontiers.
  * [27] Marilyn Gatica, Rodrigo Cofré, Pedro AM Mediano, Fernando E. Rosas, Patricio Orio, Ibai Diez, Stephan P. Swinnen, and Jesus M. Cortes.  High-order interdependencies in the aging brain.  Brain connectivity, 2021.  Publisher: Mary Ann Liebert, Inc., publishers 140 Huguenot Street, 3rd Floor New ….
  * [28] Terry Bossomaier, Lionel Barnett, Michael Harré, and Joseph T. Lizier.  Transfer entropy.  In An introduction to transfer entropy, pages 65–95. Springer, 2016\.
  * [29] Joseph T. Lizier, Mikhail Prokopenko, and Albert Y. Zomaya.  Local information transfer as a spatiotemporal filter for complex systems.  Physical Review E, 77(2):026110, 2008.  Publisher: APS.
  * [30] Sebastiano Stramaglia, Tomas Scagliarini, Yuri Antonacci, and Luca Faes.  Local granger causality.  Physical Review E, 103(2):L020102, 2021.  Publisher: APS.
  * [31] Satosi Watanabe.  Information theoretical analysis of multivariate correlation.  IBM Journal of research and development, 4(1):66–82, 1960.  Publisher: IBM.
  * [32] Han Te Sun.  Nonnegative entropy measures of multivariate symmetric correlations.  Information and Control, 36:133–156, 1978.  Publisher: Elsevier.
  * [33] Giulio Tononi, Olaf Sporns, and Gerald M. Edelman.  A measure for brain complexity: relating functional segregation and integration in the nervous system.  Proceedings of the National Academy of Sciences, 91(11):5033–5037, 1994.  Publisher: National Acad Sciences.
  * [34] For a comparison between the original measure proposed in Ref. [33] and the O-information, please see Ref. [25].
  * [35] While low-order constraints impose strong restrictions on the system and allow little amount of shared information between variables, high-order constraints impose collective restrictions that enable large amounts of shared randomness.
  * [36] Robert M Fano.  Transmission of information: a statistical theory of communications.  Mit Press, 1968.
  * [37] Craig Stuart Sapp.  Online database of scores in the humdrum file format.  In ISMIR, pages 664–665, 2005.
  * [38] Please note that regularisation methods — such as Laplace smoothing — can have significant effects on the results. We decided not to use such methods, as some chords (e.g. C-C$\sharp$-D-D$\sharp$) are not representative of the Baroque repertoire.
  * [39] For a description of consonant and dissonant intervals, please see next section.
  * [40] Imre Lahdelma and Tuomas Eerola.  Cultural familiarity and musical expertise impact the pleasantness of consonance/dissonance but not its perceived tension.  Scientific reports, 10(1):1–11, 2020.  Publisher: Nature Publishing Group.
  * [41] German language has four cases: nominative (subject), accusative (direct object), dative (indirect object), and genitive (possessive).
  * [42] Giulio Tononi, Gerald M. Edelman, and Olaf Sporns.  Complexity and coherency: integrating information in the brain.  Trends in cognitive sciences, 2(12):474–484, 1998.  Publisher: Elsevier.
  * [43] Peter E. Latham and Sheila Nirenberg.  Synergy, redundancy, and independence in population codes, revisited.  Journal of Neuroscience, 25(21):5195–5206, 2005.  Publisher: Soc Neuroscience.
  * [44] Elad Ganmor, Ronen Segev, and Elad Schneidman.  Sparse low-order interaction network underlies a highly correlated and learnable neural population code.  Proceedings of the National Academy of sciences, 108(23):9679–9684, 2011.  Publisher: National Acad Sciences.
  * [45] Andrea I. Luppi, Pedro AM Mediano, Fernando E. Rosas, Negin Holland, Tim D. Fryer, John T. O’Brien, James B. Rowe, David K. Menon, Daniel Bor, and Emmanuel A. Stamatakis.  A synergistic core for human brain evolution and cognition.  BioRxiv, 2020.  Publisher: Cold Spring Harbor Laboratory.
  * [46] Andrea I. Luppi, Pedro AM Mediano, Fernando E. Rosas, Judith Allanson, John D. Pickard, Robin L. Carhart-Harris, Guy B. Williams, Michael M. Craig, Paola Finoia, and Adrian M. Owen.  A synergistic workspace for human consciousness revealed by integrated information decomposition.  BioRxiv, 2020.  Publisher: Cold Spring Harbor Laboratory.
  * [47] Fernando E. Rosas, Pedro AM Mediano, Henrik J. Jensen, Anil K. Seth, Adam B. Barrett, Robin L. Carhart-Harris, and Daniel Bor.  Reconciling emergences: An information-theoretic approach to identify causal emergence in multivariate data.  PLOS Computational Biology, 16(12):e1008289, 2020.  Publisher: Public Library of Science San Francisco, CA USA.
  * [48] Luciano Telesca and Michele Lovallo.  Revealing competitive behaviours in music by means of the multifractal detrended fluctuation analysis: application to Bach’s Sinfonias.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2134):3022–3032, 2011.  Publisher: The Royal Society Publishing.
  * [49] Alfredo González-Espinoza, Hernán Larralde, Gustavo Martínez-Mekler, and Markus Müller.  Multiple scaling behaviour and nonlinear traits in music scores.  Royal Society open science, 4(12):171282, 2017.  Publisher: The Royal Society Publishing.
  * [50] Alfredo González-Espinoza, Gustavo Martínez-Mekler, and Lucas Lacasa.  Arrow of time across five centuries of classical music.  Physical Review Research, 2(3):033166, 2020.  Publisher: APS.
  * [51] Conor Finn and Joseph T. Lizier.  Pointwise partial information decomposition using the specificity and ambiguity lattices.  Entropy, 20(4):297, 2018.  Publisher: Multidisciplinary Digital Publishing Institute.
  * [52] Pedro AM Mediano, Fernando Rosas, Robin L. Carhart-Harris, Anil K. Seth, and Adam B. Barrett.  Beyond integrated information: A taxonomy of information dynamics phenomena.  arXiv preprint arXiv:1909.02297, 2019.