Main

A neutron can be bound either in an atomic nucleus or in a neutron star. The free neutron has a lifetime of just under 15 min and decays into a proton, electron and antineutrino. The system made of two neutrons, the dineutron, is known to be unbound by only about 100 keV. Whether multi-neutron systems can exist as weakly bound states or very short-lived unbound resonant states has been a long-standing question1. The next simplest system of three neutrons is less likely to exist owing to the odd number of nucleons and therefore weaker binding; yet, a recent calculation has suggested its existence5. Following these considerations, the four-neutron system, the tetraneutron, is an appropriate candidate to address this question. An overview of previous experiments and theoretical approaches is given in ref. 1. 

Numerous attempts have been made to find a hint for the existence of the tetraneutron as a bound or resonant state. Among these attempts, experiments were performed searching for possible bound tetraneutrons produced in uranium fission reactions (see, for example, ref. 6). Other attempts, sensitive to both bound and resonant states, used pion-induced double-charge-exchange (DCX) reactions, mainly the \({}^{4}{\rm{H}}{\rm{e}}({{\rm{\pi }}}^{-},\,{{\rm{\pi }}}^{+})\) reaction (see, for example, ref. 7), as well as transfer reactions such as \({}^{8}{\rm{H}}{\rm{e}}({\rm{d}},{}^{6}{\rm{L}}{\rm{i}})\) (ref. 8). None of the experiments yielded a positive signal.

Most of the past experiments were performed with stable nuclei. Towards the twenty-first century, with the development of radioactive-ion beam facilities, it became possible to use extremely neutron-rich nuclei in which one can expect an enhanced formation of a tetraneutron system. The first indication for a possible bound tetraneutron was reported in 20022 from a break-up reaction of 14Be into 10Be + 4n. The result stimulated several theoretical studies, all agreeing on the same conclusion: a bound tetraneutron state cannot be obtained theoretically without significantly changing our understanding of the nuclear forces9,10,11. However, the possibility of the four-neutron system existing as a resonant quasi-bound state with a very short lifetime on the order of a few 10−22 s, before decaying, has remained an open and challenging question. It was later found that the result reported in ref. 2 is also consistent with such a resonant state with the limit on its energy \({E}_{{\rm{r}}}\lesssim 2\,{\rm{MeV}}\) (ref. 3).

A decade later, in 2016, an indication of a tetraneutron resonance was reported4. A DCX reaction was used, but in contrast to previous attempts, this time the reaction was induced by a high-energy 8He radioactive beam. 8He is the most neutron-rich bound isotope, and the 8He(4He, 8Be) reaction channel was investigated. The advantage of using a radioactive beam is the freedom of selecting the reaction partner in a so-called recoil-less production (without momentum transfer) of the four-neutron system. The energy of the state was found to be Er = 0.8 ± 1.4 MeV, and an upper limit on its width was estimated as Γ ≤  2.6 MeV. However, owing to the large experimental uncertainty, the possibility of a bound state could not be excluded by this experiment.

In this work, we used the quasi-elastic knockout of an α-particle (4He nucleus) from a high-energy 8He projectile induced by a proton target to populate a possible tetraneutron state. The inverse-kinematics knockout reaction \({}^{8}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{He}})\) at large momentum transfer is well suited because the 8He nucleus has the pronounced cluster structure of an α-core (4He) and four valence neutrons with small 4n centre-of-mass motion, such that after the sudden removal of the α-particle, a rather localized four-neutron system with small relative energy between the neutrons is produced, which may have a large overlap with a tetraneutron state12,13. The chosen kinematics at large momentum transfer between the proton and the α-particle ensures that the four-neutron system will recoil only with the intrinsic momentum of the 4He core in the 8He rest frame, without any further momentum transfer, thus allowing the recoil-less production. Furthermore, final-state interactions between the four neutrons and the charged particles are also minimized owing to the large momentum transfer, separating charged reaction partners from the neutron spectators in momentum space (Fig. 1).

Fig. 1: Schematic illustration of the quasi-elastic reaction investigated in this work.
figure 1

Top: quasi-elastic scattering of the 4He core from a 8He projectile off a proton target in the laboratory frame. The length of the arrows represents the momentum per nucleon (the velocity) of the incoming and outgoing particles. Zbeam is the beam axis. Bottom: the equivalent p–4He elastic scattering in their centre-of-mass frame, where we consider reactions at backward angles close to 180°, θc.m. ≳ 160°. In this frame, the momentum of the proton balances that of the 4He, \({{\bf{P}}}_{{\rm{p}}}=-{{\bf{P}}}_{{}^{4}{\rm{He}}}\), that is, the proton is four times faster than the 4He.

Full size image

The experiment took place at the Radioactive Ion Beam Factory operated by the RIKEN Nishina Center and the Center for Nuclear Study, University of Tokyo, using the Superconducting Analyzer for Multi-particles from Radio Isotope Beams (SAMURAI)14. A primary beam of 18O was directed onto a beryllium production target producing a cocktail of radioactive nuclei from fragmentation. The secondary 8He beam was separated using the BigRIPS fragment separator and transported with an energy of 156 MeV per nucleon to a 5-cm-thick liquid-hydrogen target15 located at the SAMURAI spectrometer (Fig. 2).

Fig. 2: Experimental set-up and charged fragments momenta.
figure 2

Left: schematic view of the experimental set-up. The 8He secondary beam at 156 MeV per nucleon is transported from the BigRIPS (Big RIKEN projectile-fragment separator) into the SAMURAI set-up, where it hits a liquid-hydrogen (LH2) target. In a quasi-elastic \(({\rm{p}},\,{{\rm{p}}}^{4}{\rm{He}})\) reaction, the 4He core is knocked out from the 8He projectile. Scintillator detectors and drift chambers are used for beam identification and tracking. The trajectories of the outgoing fragments are tracked by three silicon (Si) planes and bent afterwards through the SAMURAI spectrometer towards the focal-plane detectors. Two neutron-detector arrays were placed at a forward angle behind the SAMURAI. An additional scintillator wall was placed at smaller bending angle to detect the unreacted 8He beam. Right: measured momenta of the knocked-out 4He and the scattered proton after the quasi-elastic scattering (symbols). The momentum distribution of the incoming 8He beam is shown for comparison. The solid curves are the results from the simulation. The cyan (magenta) dotted line represents the upper (lower) limit of the 4He (proton) momentum expected from the simulation assuming a quasi-elastic scattering, and the orange line indicates the central beam momentum. 

Source data

Full size image

The incoming beam was measured upstream of the target on an event-by-event basis using scintillators for charge identification as well as momentum measurement, and two drift chambers for tracking (Extended Data Fig. 1).

The outgoing charged fragments (α-particle and proton) emerging from the quasi-elastic scattering were detected using a combination of detectors downstream of the target. Three planes of silicon-strip detectors, where each plane consists of two orthogonal layers enabling position measurements in both horizontal and vertical directions, served for tracking, energy-deposition measurement and reconstruction of the reaction vertex inside the target (Extended Data Figs. 2 and 3).

Behind the silicon planes, both charged fragments were bent through the magnetic field of the SAMURAI spectrometer, which was operated at a nominal magnetic field of 1.25 T in the centre of the magnet. The experiment was designed to detect an α-particle and a proton that emerge from quasi-elastic scattering close to 180° in the centre-of-mass frame (Fig. 1). Under these kinematical conditions, their resulting outgoing momenta are very different from each other in the laboratory frame, as shown in Figs. 1 and 2. The knocked-out α-particle is slowed down from its initial momentum, that is, with the incoming beam velocity, to a momentum of about 330 MeV/c per nucleon after the reaction (where c is the speed of light). In contrast, the proton, which was at rest in the initial state, becomes the fastest particle in the reaction, gaining a typical momentum of about 860 MeV/c. At the focal plane, a drift chamber is used for tracking of the fragments after the magnet, and two scintillator walls located side by side, which cover a wide momentum range, are used for energy-deposition and time-of-flight measurements. The α-particle and proton are identified from a combination of their measured energy deposition, each in a different scintillator wall, and their position in the drift chamber (Extended Data Fig. 4). Their momenta are determined precisely from their reconstructed trajectories through the SAMURAI spectrometer.

As no additional momentum is transferred to the neutrons in the reaction, they continue moving with nearly beam velocity and can be detected, in principle, by the neutron detectors placed at a forward angle behind the SAMURAI spectrometer. The detection efficiency for neutrons is much lower than that for charged particles and decreases quickly as a function of the number of detected neutrons. The small p–4He elastic cross-section at backwards centre-of-mass angles of less than 1 microbarn (ref. 16) resulted in the relatively low statistics of 422 events obtained for the \({}^{8}{\rm{H}}{\rm{e}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{H}}{\rm{e}})\) reaction. These factors made it impossible to detect more than two neutrons in coincidence with the charged particles. Therefore, the neutron detection is not a part of the current study, aside from a consistency check (provided in Supplementary Information) of the near recoil-less production of the free neutrons.

The combined selection of event-by-event identification of incoming 8He-beam particles in coincidence with the knocked-out α-particle and the scattered proton defines the \({}^{8}{\rm{H}}{\rm{e}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{H}}{\rm{e}})\) channel. From a precise measurement of the momenta of the charged particles, the energy spectrum of the 4n system is reconstructed assuming energy and momentum conservation through the missing mass:

$${E}_{4{\rm{n}}}=\sqrt{{E}_{{\rm{miss}}}^{2}-{{\bf{P}}}_{{\rm{miss}}}^{2}}-4{m}_{{\rm{n}}},$$
(1)

where Emiss (Pmiss) is the energy (momentum) component of the missing-momentum four-vector, and mn is the neutron mass. Using this notation, a bound 4n system will appear at E4n < 0 whereas a resonant state will appear at E4n > 0. The missing momentum in equation (1) is defined by \({\bar{P}}_{{\rm{miss}}}={\bar{P}}_{{}^{8}{\rm{He}}}+{\bar{P}}_{{\rm{p}}({\rm{tgt}})}-{\bar{P}}_{{}^{4}{\rm{He}}}-{\bar{P}}_{{\rm{p}}}\), where the four-momenta \(\bar{P}\) on the right-hand side of the equation are those of the incoming beam, target proton, knocked-out α-particle and scattered proton, respectively.

The \({}^{6}{\rm{H}}{\rm{e}}({\rm{p}},{{\rm{p}}}^{4}{\rm{H}}{\rm{e}})\) knockout reaction was measured with almost exactly the same experimental conditions as for 8He, except for some small differences in the energy of the incoming beam and the beam profile (Supplementary Table 2), and served as a benchmark for verifying the analysis and calibration procedures. In the case of 6He, the 2n system is produced by the sudden removal of the 4He core. The two-neutron relative-energy spectrum is expected to be well described by theory taking into account both the well established ground-state wavefunction and the final-state scattering wave of the two neutrons, predicting a low-energy peak around 100 keV. Similarly to the 8He case, we define the missing mass (\({\bar{P}}_{{}^{8}{\rm{He}}}\to {\bar{P}}_{{}^{6}{\rm{He}}}\) and \(4{m}_{{\rm{n}}}\to 2{m}_{{\rm{n}}}\)).The measured missing-mass spectrum for 6He is shown in the right panel of Fig. 3 together with the theoretical calculation17 convoluted with the experimental acceptance and resolution (blue curve). The energy range shown represents the one covered by the experimental acceptance. The calculation is compared with the data by implementing it into an event generator for the quasi-elastic reaction, which uses the measured p–4He differential elastic cross-section16 as an input, as well as the measured internal momentum distribution of the α-particle in 6He (ref. 18). The generated events are transported through the experimental set-up in Geant4 simulations to account for the experimental acceptance and detector resolutions. The excellent agreement of the simulated theoretical distribution with the measured spectrum confirms the analysis and the calibration for determining the missing mass. The missing-mass resolution obtained in the measurement is approximately 1 MeV sigma, and is almost constant over the measured energy range. The systematic uncertainty for the determination of the absolute energy was estimated from this measurement to be 0.4 MeV and that of the energy width to be 0.27 MeV (Methods). Also shown in the right panel of Fig. 3 (green curve) is a possible small background contribution coming from two-step process where 4He is produced in a first step (see Methods and following discussion for 8He). This background was estimated from the measured cross-section to contribute 1% of the total number of measured events.

Fig. 3: Missing-mass spectra.
figure 3

Left: missing-mass spectrum of the four-neutron system extracted from the \({}^{8}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{He}})\) reaction. The different curves represent a Breit–Wigner resonance (pink), a non-resonant continuum (dashed blue), the background from two-step processes (green) and the total sum (solid blue). Right: missing-mass spectrum of the two-neutron system extracted from the \({}^{6}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{He}})\) reaction. The blue curve represents the theoretical calculation17 convoluted with the experimental acceptance and resolution, and the green curve represents the background from the two-step reaction.

Source data

Full size image

The measured missing-mass spectrum of the four-neutron system from the \({}^{8}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{H}}{\rm{e}})\) reaction is shown in the left panel of Fig. 3. Two components are observed: a well pronounced peak in the low-energy region with an energy around 2 MeV and a broad distribution at higher energies attributed to a non-resonant continuum response13, a direct four-body decay.

The shape of the non-resonant continuum spectrum of the four neutrons has been studied theoretically for the case of the four-neutron structure formed after the sudden removal of the α-core from 8He (ref. 13). The creation of the system is investigated by introducing into the Schrödinger equation a source term that accounts for the reaction mechanism producing the four-body system, and that depends explicitly on the internal structure of the parent nucleus 8He. The 8He ground-state wavefunction (without final-sate interaction) was treated using the five-body (\({}^{4}{\rm{He}}+4{\rm{n}}\)) cluster orbital shell model approximation (COSMA)12. The exact shape of the non-resonant continuum is sensitive to the hyperradius of the source, ρsour an internal radius of the 4n system, described in the hyperspherical harmonics basis. A hyperradius of 5.6 fm is considered by the theory as the most realistic, as it reproduces the correct experimental radius of 8He in the COSMA model. This results in a broad distribution centred around 30 MeV, in good agreement with the observed experimental spectrum.

We model the spectrum as follows:

$$f({E}_{4{\rm{n}}})=a{f}_{{\rm{res}}}({E}_{4{\rm{n}}})+b{f}_{{\rm{con}}}({E}_{4{\rm{n}}})+c{f}_{{\rm{bkg}}}({E}_{4{\rm{n}}}),$$
(2)

where a, b, and c are constants, fres is a Breit–Wigner function representing the possible resonance structure, and fcon is the non-resonant continuum part presented above with the hyperradius as a parameter. The last term in equation (2), fbkg, represents possible background events coming from competing processes. Several processes were investigated and quantified (Methods), where the only non-negligible contribution found is from a two-step process involving 6He (4He) production: proton-induced break-up of 8He into 6He (4He) followed by a p–4He quasi-elastic scattering. The resulting energy distribution is broadened and shifted to lower energies compared with the pure 6He case (right panel of Fig. 3) owing to the two-step process, which has been taken into account in the simulation of fbkg. This background was estimated from measured cross-sections to contribute 2.6% to the total number of measured events (Methods), which has been used to determine the normalization constant c.

The experimental spectrum was then fitted with the energy-dependent function given in equation (2), where the fit function was convoluted with the experimental response, taking into account acceptance and detector resolutions. The experimental acceptance is not constant over the measured energy range. It is maximal in the region \(20\,{\rm{MeV}} < {E}_{4{\rm{n}}} < 40\,{\rm{MeV}}\) (Extended Data Fig. 5).

The result of the χ2 minimization is presented by the solid blue curve in the left panel of Fig. 3, together with the individual contributions. The statistical significance of the peak structure is well beyond the 5σ level (Methods).

The probability of populating a four-neutron system in a resonant state after the sudden removal of the α-core in 8He is determined by the overlap of the 4n wavefunction in the final state and the relative motion of the four neutrons in the 8He initial state. This overlap integral defines the ratio of cross-sections for the population of the resonance and the non-resonant continuum. Unconvoluting with the acceptance of the set-up, following the energy dependence of equation (2), we extract a probability of Pr = 18.7 ± 2.3%. For comparison, the relative motion of the four neutrons studied in the COSMA model12,13 yields a probability of about 30%. This value is obtained by considering the hyperradius of 5.6 fm, whereas the resulting value from the fit to the experimental data is 5.0 ± 0.1 fm, which would yield a smaller probability to populate the resonant state.

Assuming a resonant state, its energy and width as determined from the fit are

$$\begin{array}{c}{E}_{{\rm{r}}}=2.37\pm 0.38({\rm{stat}}.)\pm 0.44({\rm{sys.}})\,{\rm{MeV}},\\ \varGamma =1.75\pm 0.22({\rm{stat}}.)\pm 0.30({\rm{sys.}})\,{\rm{MeV}}.\end{array}$$

For comparison, Fig. 4 shows our result (full red symbol) together with the previous experimental result obtained from the DCX measurement4 (open red symbol). The energy of the resonance is in agreement within the uncertainty, despite the fact that different reactions were used to probe the 4n system, and is also in agreement with the upper limit given in ref. 3.

Fig. 4: Comparison of experimental results with theory predictions.
figure 4

Energy versus width of a tetraneutron resonance. Experimental data are shown in red: this work (full symbol), and the result from the DCX measurement4 (open symbol), where the red arrow indicates that the measured width is an upper limit. Theory predictions are shown in blue based on: NCSM19 and ref. 20 cited in ref. 20 (full stars), NCGSM (open star)20 (cross)21, where the blue arrow indicates that the width is predicted to be larger than 3.7 MeV, and QMC calculations5 (band). Whether this observation of a low-energy peak is attributed to a four-neutron resonant state or to other correlations between the neutrons in the final state, needs to be clarified by ab initio theories.

Source data

Full size image

From the theory side, there is no consensus among the different theories and their predictions are partly contradictory, although, there is a general agreement that a bound 4n does not exist. In 2003, Pieper10 studied this possibility using Green’s function Monte Carlo calculations. His conclusion was that the existence of a bound 4n state has to be excluded, unless nuclear forces are drastically modified. However, his calculation suggested that a possible resonance might exist near 2 MeV, but in such a case it must be very broad.

Using a similar approach, the quantum Monte Carlo (QMC) framework based on two-body and three-body chiral interactions was used to calculate the energy of a 4n resonance5. The result supports the existence of a resonant state with an energy of 2.1(2) MeV, while no prediction has been made for its width (blue band). An extended no-core shell model (NCSM) approach using a harmonic-oscillator representation predicts different resonant states (full stars) including their corresponding widths19 (see also ref. 20 cited in ref. 20). Calculations have been performed also in the framework of the no-core Gamow shell model (NCGSM)21. These resulted in Er ≈ 7 MeV and Γ ≈ 3.7 MeV (cross), where the conclusion was in fact that the energy of a 4n resonance might be compatible with the experimental value of ref. 4, albeit with a significantly larger width. As pointed out in a later study20, these calculations were incomplete, as they were performed only in truncated model spaces or with unphysically overbinding interactions. The authors of this work20 concluded that both the energy and the width of such a resonance are comparable with the experimental data (open star). At the same time, other calculations claim that to generate such a resonance, nuclear forces have to be significantly changed22,23,24,25,26,27, which would not be consistent with the present understanding. We note that some theories26,27 predict a non-resonant low-energy enhancement of the density of states in the four-neutron spectrum. Whether such a prediction is consistent with our observed resonance-like feature cannot be currently ascertained, as the energy spectrum of the four-neutron system is not given. The drastically different predictions resulting from different theoretical approaches highlight the importance of the current firm experimental observation.

In conclusion, we have presented the experimental observation of a resonance-like structure consistent with a tetraneutron state near threshold after 60 years of experimental attempts to clarify the existence of this state. The use of a high-energy knockout reaction in inverse kinematics allowed a precise measurement. The use of a radioactive 8He beam as the parent system and a direct, large momentum-transfer reaction opened up the opportunity to create the 4n system in a well defined one-step process and in a recoil-less undisturbed way. The optimized detection system enabled a precise determination of the final state and a high-resolution measurement. A well developed peak structure has been observed at an energy of 2.37 ± 0.38(stat.) ± 0.44(sys.) MeV with a striking statistical level. This is in agreement with the result of ref. 4 and the upper limit given in ref. 3. Both the energy and the extracted width of Γ = 1.75 ± 0.22(stat.) ± 0.30(sys.) MeV are consistent with a tetraneutron state that is unbound with a corresponding lifetime of (3.8 ± 0.8) × 10−22 s. Next-generation experiments using different reaction mechanisms and possibly detecting the four neutrons in coincidence will reveal more insights into the properties of the four-neutron system, including correlations among the neutrons. More elaborated ab initio nuclear theories accounting fully for the effect of the continuum are necessary to understand the observed low-energy peak and its origin.

Methods

Energy E 2n of the 2n system

The theoretical calculation used as an input for the E2n distribution includes both the ground-state wavefunction and the final-state interaction between the two neutrons17. These were obtained from a three-body (\({}^{4}{\rm{He}}+2{\rm{n}}\)) cluster model, which contains a phenomenological three-body force as well as the following local angular momentum l-dependent two-body interactions: a central s-wave interaction in the nn system, and central s- p- and d-wave interactions in the n–α system, as well as n–α spin–orbit interactions on the p- and d-waves. The ground-state wavefunction is computed in a three-body model using the computer code FaCE28, an acronym for Faddeev with Core Excitation, whereas the final-state interaction is taken into account by using the nn transition t-matrix approach. The calculated distribution is shown in Extended Data Fig. 6 (purple). For comparison, we also show the calculated distribution without the nn final-state interaction (green), which reflects the nn motion in the source, the 6He ground state. The distributions are normalized to the same maximum, such that only their shape is compared. It can be seen that the effect of the nn final-state interaction is very large. The ground-state distribution peaks at a larger energy with a far extending tail towards higher energies.

Systematic uncertainty of the missing-mass determination

As the theoretical calculation used as an input for the energy of the two-neutron system in the quasi-elastic \({}^{6}{\rm{He}}({\rm{p}},\,{\rm{p}}{}^{4}{\rm{He}})\) reaction is considered to be accurate, the comparison of it with the measured spectrum is used to estimate the systematic uncertainty of the missing-mass measurement. Different offsets covering the range of ±0.5 MeV are added to the generated theoretical distribution and for each offset the goodness of the agreement with the data is quantified by calculating the χ2 value. The generated distributions are convoluted with the energy response of the set-up, to account for the experimental acceptance and the detector resolutions. The minimum χ2 obtained is located very close to zero, that is, very close to the original distribution (Extended Data Fig. 7). The range of deviations in the energy scale still showing reasonable agreement with the data (\({\chi }_{{\rm{\min }}}^{2}\pm 1\) range) is interpreted as the systematic uncertainty for the absolute missing-mass determination. This results in an uncertainty with an average value of ±0.4 MeV. To determine the systematic uncertainty for the missing-mass width, the original theoretical distribution is used (that is, offset of zero). To check the sensitivity to the resolution, different resolutions are applied in the range \(\sigma =1\pm 0.5\,{\rm{MeV}}\,.\) The generated distributions are convoluted with the acceptance of the set-up and are smeared using a Gaussian function with a width σ. Similarly as described above, the χ2 values are calculated and the \({\chi }_{{\rm{\min }}}^{2}\pm 1\) range is interpreted as the systematic uncertainty (Extended Data Fig. 7). This results in an uncertainty with an average value of ±0.27 MeV. We note that for both absolute value and resolution, the minimum χ2 values fall very close to the expected values (which could be statistically by accident). As the given systematic uncertainties are determined solely by the statistics of the 6He data, we consider the given uncertainties as conservative estimates.

An additional systematic uncertainty from the background function in equation (2), on both the energy and the width, is added quadratically. The fitting procedure is repeated three times, with the green curve presented in Fig. 3 and with the 1σ lower and upper limits (Supplementary Information). The values for the energy and the width are taken as the average between the different fits, and the systematic uncertainty is taken as the difference between them. This results in an uncertainty of 0.18 MeV and 0.14 MeV for the energy and width, respectively.

Energy E 4n of the 4n system

The relative-energy spectrum of four neutrons produced by sudden removal of the α-particle from 8He as a source was studied theoretically in ref. 13. The final-state energy spectrum depends on the intrinsic 4n relative motion in 8He, that is, the 8He ground-state wavefunction as the source (without final-state interaction), and the final scattering state of the four isolated neutrons, that is, 4n final-state interaction. The wavefunction used in ref. 13 for 8He is based on the COSMA model12, and the source term contains the Fourier transform of the overlap between 8He and the α-particle wavefunctions. The method of hyperspherical harmonics29 was used to solve the equations of the model. This method is based on the link between the hyperspherical-function method and the oscillator NCSM and uses a Slater determinant representation of the hyperspherical harmonics functions. The resulting source function is described by the internal variables of the four-neutron system, hyperradius, hyperangle, hyperangular momentum and the 4He–4n relative motion. Calculations showed that only hyperangular momentum K = 2 significantly populates the 4n non-resonant continuum, reflecting the motion of the four neutrons in 8He. The resulting energy spectrum in ref. 13 is presented only up to 20 MeV. It was also used in ref. 4 to estimate the non-resonant background and can be modelled analytically as

$${f}_{{\rm{cont}}}({E}_{4{\rm{n}}})={E}_{4{\rm{n}}}^{\alpha }\times \exp (-{E}_{4{\rm{n}}}/{\epsilon }_{{\rm{a}}}),$$
(3)

where α = 7/2 + K with K = 2 as stated above. K = 0 and α = 7/2 corresponds to the free four-body phase space. \({\epsilon }_{{\rm{a}}}\approx \frac{3.3{\hbar }^{2}}{{m}_{{\rm{N}}}{\rho }_{{\rm{sour}}}^{2}},\) where mN is the nucleon mass, ρsour is the hyperradius, ħ is the reduced Planck constant and the factor 3.3 was estimated by matching to the calculated distribution in ref. 13.

Background events from competing reactions

Events with a proton and an α-particle in the final state might result from other competing reactions. However, owing to the unique kinematics of quasi-elastic p–4He scattering at large centre-of-mass angles, other direct reactions are excluded. Examples are one-neutron knockout, \({}^{8}{\rm{He}}{({\rm{p}},{\rm{p}}{\rm{n}})}^{7}{\rm{He}},\) or 6He knockout, \({}^{8}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{6}{{\rm{He}}}^{* })2{\rm{n}}.\) In the former case, the momentum of the proton is too small to be selected by the experiment, where in the latter, it is too high at the angular range covered by the set-up, owing to the larger mass of 6He relative to 4He Therefore, the only possible contribution of background events to our measured p–4He events can come from secondary reactions. Below we list the possible reactions and their expected contributions.

\({}^{8}{\rm{He}}{({\rm{p}},{\rm{pn}})}^{7}{\rm{He}}\). A first interaction is single-neutron knockout from 8He leading to break-up into 6He and two neutrons. This can be followed by a second interaction, quasi-elastic p–4He scattering at backward centre-of-mass angles. These processes will result in events with a missing-mass distribution similar to that of the one-step \({}^{6}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{4}{\rm{He}})\) reaction but broadened and shifted to smaller energies. This is due to the difference in the separation energies of two and four neutrons in 6He versus 8He of about 2 MeV, and due to the additional recoil (Supplementary Information). The single-neutron knockout cross-sections from 6,8He were measured in an experiment performed within the same experimental campaign as the experiment presented here, at a slightly higher beam energy of 185 MeV per nucleon using a series of targets with different nuclear charges Z (ref. 30). For 6He, the cross-section was extracted for the hydrogen target, whereas for 8He, the lowest-Z target used was carbon. We therefore scale the measured cross-section for 6He with hydrogen, \({\sigma }_{{}^{6}{\rm{He}},{\rm{H}}}\), by the carbon measurement \(({\sigma }_{{}^{8}{\rm{He}},{\rm{C}}}/{\sigma }_{{}^{6}{\rm{He}},{\rm{C}}})=1.26\pm 0.16.\). Using \({\sigma }_{{}^{6}{\rm{He}},{\rm{H}}}=47\pm 4\,{\rm{mb}}\) results in \({\sigma }_{1{\rm{n}}}=59\pm 9\,{\rm{mb}}\).

\({}^{8}{\rm{He}}{({\rm{p}},{\rm{p}})}^{8}{{\rm{He}}}^{* }\). A first interaction is inelastic excitation of 8He leading to a break-up into 6He and two neutrons. These cross-sections were measured as well in ref. 30, with scaling ratio of \({\sigma }_{{}^{8}{\rm{He}},{\rm{C}}}/{\sigma }_{{}^{6}{\rm{He}},{\rm{C}}}=1.02\pm 0.13.\) Using \({\sigma }_{{}^{6}{\rm{He}},{\rm{H}}}=11.0\pm 0.9\,{\rm{mb}}\) gives \({\sigma }_{{\rm{inel}}}=11.2\pm 1.7\,{\rm{mb}}\).

In total, the cross-section evaluated for the two processes described above equals \({\sigma }_{{\rm{tot}}}=71.1\pm 9.2\,{\rm{mb}}\,.\) Evaluating the relative contribution to the number of one-step p–4He quasi-elastic reactions Nreac at half of the target thickness (t) gives \({N}_{{\rm{bkg}}}/{N}_{{\rm{reac}}}\,=\sigma ({\rm{b}})\times t({\rm{g}}\,{{\rm{cm}}}^{-2})\times \) \(0.6({\rm{Avogadro}})/{\rm{A}}=0.071\times 0.175\times 0.6\)\(=\,0.75\pm 0.10 \% \), where A is the mass number of the target. These types of processes can also occur from 8He break-up along beamline materials before the target, plastic scintillators, kapton windows and air. These contributions were estimated starting from the last plastic scintillators at the BigRIPS fragment separator up to the kapton window at the entrance to the target cell, using the measured cross-sections for hydrogen and carbon targets30. The produced 6He has an angular spread in the transverse direction according to its intrinsic momentum, which has a width of \({\sigma }_{{}^{6}{\rm{He}}}=35\,{\rm{MeV}}/c\) (ref. 18). Assuming a central beam energy of 156 AMeV, its corresponding momentum is \({P}_{{}^{6}{\rm{He}}}=3.37\,{\rm{GeV}}/c\), leading to an angular θ width of \({\sigma }_{\theta }={\sigma }_{{}^{6}{\rm{He}}}/{P}_{{}^{6}{\rm{He}}}=0.035/3.37=10\,{\rm{mrad}}.\) The spread in the x and y directions evaluated at the entrance to the target equals \({\sigma }_{x,y}=d\cdot {\sigma }_{\theta }=d({\rm{mm}})\times 0.01({\rm{rad}})\,,\) where d is the distance from the material at which 6He was produced to the target entrance. The target enclosure radius equals 20 mm. Therefore, depending on the distance of production point and target, some fraction of the produced 6He is expected to not reach the target, which reduces its intensity. Using the measured distances to the target entrance, the estimated fraction of 6He produced before the target and reaching the effective target volume is 0.36 ± 0.01%. Overall, the estimated contribution from these two processes is 1.11 ± 0.10%.

\({}^{8}{\rm{He}}({\rm{p}},\,{{\rm{p}}}^{6}{\rm{He}})2{\rm{n}}\). A first interaction is 6He knockout to its ground state. This process can contribute only when the 6He is produced along the target, as, owing to its angular spread, it will not reach the target region in case of production in beamline materials. The cross-section for this process is not well known, and was studied previously only for a much higher energy18. However, it was shown to be compatible with the p–4He elastic scattering cross-section. We therefore adapt the measured p–4He total cross-section, which was measured at 156 MeV (ref. 16), σ = 91.8 mb, and use it as an upper limit. The relative contribution of these events is evaluated as \(0.0918\times 0.175\times 0.6=0.96 \% \).

\({}^{8}{\rm{He}}+{\rm{p}}{\to }^{4}{\rm{He}}\) . In the first interaction, 4He is produced from fragmentation of 8He. This can be, for example, a result of elastic p–4He scattering or inelastic excitation of 8He. The inclusive 4He cross-section was measured at higher energy for \({}^{8}{\rm{He}}+{}^{12}{\rm{C}}\to {}^{4}{\rm{He}}\) to be 95 ± 5 mb (ref. 31). We scale it by a factor of 1/2 to estimate the cross-section for the hydrogen target, leading to σ = 48 ± 5 mb. The relative contribution of these events is evaluated as \(0.048\times 0.175\times 0.6=0.50\pm 0.05 \% .\) This process will result in events with a negative missing-mass distribution below −3.1 MeV, the binding energy of 4He in 8He that will extend to more negative values owing to the additional recoil (Supplementary Information).

\({}^{8}{\rm{He}}{({\rm{p}},{\rm{p}})}^{8}{\rm{He}}\). A first interaction is p–8He elastic scattering at backwards centre-of-mass angles. This will lead to a fast proton in the final state. In a second interaction, 4He can be produced. However, in this case, the momentum of the proton is even larger than that resulting from direct 6He knockout, owing to the larger mass; therefore, this reaction channel is excluded.

For the benchmark measurement with 6He, only one reaction can contribute to a background, where 4He is produced in the first interaction, \({}^{6}{\rm{He}}+{\rm{p}}{\to }^{4}{\rm{He}}\). Similar to 8He, the inclusive 4He cross-section was measured at higher energy for \({}^{6}{\rm{He}}+{}^{12}{\rm{C}}\to {}^{4}{\rm{He}}\) to be 189 ± 14 mb (ref. 31). We scale it by a factor of 1/2 to estimate the cross-section for the hydrogen target, leading to σ = 95 ± 7 mb. The relative contribution of these events is evaluated as \(0.095\times 0.175\times 0.6=1.00\pm 0.07 \% .\) This process will result in events with a negative missing-mass distribution below −0.975 MeV, the binding energy of 4He in 6He that will extend to more negative values owing to the additional recoil (Supplementary Information).

Summary of background contributions

From the reactions listed above, we conclude that the only contributions from background processes come from reactions in which 6He or 4He were produced in a first step, and the overall contribution is evaluated as 2.6%. In the fitting of the missing-mass spectrum, fbgk is taken as the simulated distribution for the two-step processes described above. Each process is simulated individually, and adds weight to the total distribution according to its expected contribution. In addition, as the measured cross-sections adopted to estimate the background contributions were not measured directly for the reactions of interest, we consider the uncertainty as a variation by a factor of two on each one of the scaled cross-sections (Supplementary Information).

It is noted that the resulting background distribution can also explain the two events observed at the energy region of \(-10\,{\rm{MeV}} < {E}_{4{\rm{n}}} < -\,8\,{\rm{MeV}}.\) These two events cannot result from a direct reaction, as the lower limit of a bound tetraneutron is −3.1 MeV, corresponding to the binding energy of 8He against decay into \({}^{4}{\rm{He}}+4{\rm{n}}\). Therefore, these events can be attributed only to two-step reactions. Estimating the background function at −9 MeV, we expect 0.27 events in that energy region, leading to a 2σ significance of the measured two events, which is not considered as a statistically significant deviation from the expected background contribution in that energy region.

Finally, a major difference between one-step and two-step processes is expected in the reaction-vertex distribution of the proton and the α-particle (Extended Data Fig. 8). For a two-step process, an exponentially increasing yield of reaction products along the target is expected owing to the need to produce the 6He (4He) first. Such a distribution is not observed in the data, corroborating that the background contribution is indeed small. The estimated background contribution will be considered to evaluate the significance level of the observed peak.

To quantify the statistical significance of the resonance observed in Fig. 3, we evaluate the number of background events in the peak region defined as \({E}_{{\rm{r}}}\pm 2\varGamma \) and \(-2\,{\rm{MeV}} < {E}_{4{\rm{n}}} < 6\,{\rm{MeV}}.\) The number of measured events in the peak region amounts to N = 41 events. The number of background events equals the integral of the background and continuum functions (fbkg and fcont in equation (2)) in that region, resulting in an average number of eight events from the background and two events from the continuum, such that the total expected background amounts to ten events. Even though we allowed a variation of a factor of two for the background estimation, this does not change the result significantly. Hence, we conclude that the resonance structure is observed with high statistical significance, well above the 5σ observation threshold.