F Miguel Marqués

Jaume Carbonell

The quest for light multineutron systems

Keywords: . Trineutron, Tetraneutron, Multineutron, Few-Body Systems, Few-Nucleon Problem, ab initio

The long history of the research concerning the possible existence of bound or resonant states in light multineutron systems, essentially 3 n and 4 n, is reviewed. Both the experimental and the theoretical points of view have been considered, with the aim of showing a clear picture of all the different detection and calculation techniques that have been used, with particular emphasis in the issues that have been found. Finally, some aspects of the present and future research in this field are discussed.

Introduction

Understanding the structure of nuclei and how their properties emerge from the underlying forces between nucleons is a major goal of modern nuclear physics. Light nuclei have long played a fundamental role in this respect, and those exhibiting very asymmetric neutron-to-proton ratios have proven to be particularly sensitive to details of the two-and few-body forces used in nuclear models. Therefore, the question about the existence of multineutrons, the most extreme combinations one can find, raises many experimental and theoretical challenges. The discovery of such neutral systems as bound or resonant states would have far-reaching implications for many facets of nuclear physics, from the nature of the force itself up to the way it builds nuclei, as we will see in Sec. 3, and also for the modeling of neutron stars (see for example Ref. [1]).

The quest for neutral nuclei may be traced back to the early 1960s [START_REF] Ogloblin | Nuclei Far From Stability, Treatise on Heavy-Ion Science[END_REF]. However, the last two decades have witnessed a renewed and enhanced interest in studying light neutron systems. This is essentially due to two experimental results claiming the possible observation of bound or low-lying resonant tetraneutron states [START_REF] Marqués | [END_REF]4], which motivated new experiments but mostly new theoretical studies, leading to progress in the computation of the exact solutions of the few-nucleon system in the continuum (see for example Refs. [5,6,7,8,9,10]).

The main difficulty in the study of multineutron states is that there are no particle-stable substructures, unlike standard charged nuclei. If the dineutron were bound, one could more comfortably inquire about the existence of 3 n (as n+ 2 n) or 4 n (as 2 n+ 2 n), both experimentally and theoretically. Despite some recent speculations [11] this seems however totally excluded, making the progress in this field particularly non trivial. This intrinsic difficulty explains the lack of strong experimental evidence concerning the existence of such states, and leads on the theoretical side to contradictory results for the very same systems, even when described using the same nn interaction.

As the title says, we have restricted this review paper to small clusters of neutrons, basically 3 n, 4 n and 6 n, the ones that can be realistically accessed in experiments and that can be theoretically tackled with exact ab initio methods. We have deliberately omitted all the interesting studies concerning nuclear matter and neutron droplets, involving a large number of particles, for they require a very different approach from the experimental as well as the theoretical points of view. The interested reader can find helpful reviews in Refs. [12,13,14,15,16].

The paper is divided in two main sections, corresponding to the review of the experimental (Sec. 2) and theoretical (Sec. 3) works. On the experimental side, we have tried to class the long series of very different experiments, across a wide range of techniques, into a few categories for clarity, and in doing so we have restricted ourselves to openly accessible works. In the theoretical part we have focused our efforts mainly on the results obtained during the last twenty years, which also coincide with the main experimental signals. We will close the paper with some conclusions and perspectives in the field.

Multineutron experiments

Detecting neutral particles represents an experimental challenge in general. Charged particles crossing a detector material interact with the atomic electron clouds, leading to detection efficiencies of 100%. Neutral particles, however, must interact with a nucleus, a much less probable process, leading to typical detection efficiencies (for similarly sized detectors) of few per cent. Moreover, the detection of x neutrons becomes exponentially harder, since the efficiency will decrease roughly as ε xn ≈ (ε 1n ) x , as shown in Fig. 1. To make things worse, the latter is just an upper limit due to cross-talk effects [17] (a neutron may interact several times), that require the application of rejection algorithms with subsequent losses of total efficiency.

Even if one obviates the detection problem, the fact that neutrons are difficult to guide and unstable themselves does not allow us to build a system of several neutrons out of its components. Therefore, all the experimental approaches must face the challenge of overcoming those two issues: how to build multineutrons and then how to detect them. In reviewing all the existing works since the 1960s we have followed the chronology when possible, but at the same time we have tried to organize them in subsections according to the probe they used.

Paradoxically, even if the aim of these experiments is the 'observation' of systems of several neutrons, almost all of them have in common the absence of any neutron detection, due to the aforementioned very low ε xn efficiencies. We have identified basically two main categories of experiments according to their principle:

a ○------→ b ○ ⇓ ⇑ c ○------→ ( A n) missing mass a ○------→ b ○ ⇓ ⇑ ( A n) ------→ (X) two step
In both of them one must detect in the final state only one charged particle, b ○, either directly or indirectly. On the left panel, representing the missing-mass category, the multineutron is sought to be formed in a twobody collision a ○+ c ○ leading to a two-body final state b ○+ A n. Only in that case, the detection of b ○ can sign the x (detected neutrons) (ε 1n )

x

[%]

Fig. 1. The neutron efficiency to the power of the number of detected neutrons, as a function of the latter. From bottom to top, the lines/symbols correspond to neutron efficiencies ε1n of 10, 20, 30, 40 and 50%, respectively. The vertical scale goes from 1 to 60%. In practice, due to cross-talk effects the multineutron (x > 1) efficiency will be even lower, εxn < (ε1n) x .

population of states in the missing multineutron system through the unique constraints of two-body kinematics.

The main advantages (+) and issues (-) of the missingmass technique are:

+ It only requires the detection of one charged particle. + The multineutron mass number A is well defined. + Both bound and resonant A n states can be probed.

-The cross-sections that bring all the protons into b ○ without breaking it are generally (extremely) low.

-Any beam or target contaminant different from a ○ or c ○ would lead to a missing partner(s) of b ○ that is not a multineutron.

On the right panel above, representing the two-step category, a bound multineutron is supposed to be produced in a previous reaction, and then one seeks to sign its subsequent interaction with a nucleus a ○ that induces a transformation into b ○. The advantages and issues of the two-step category are: + It only requires the detection of one charged particle.

-Only bound A n states can lead to this second reaction.

-There is no sensitivity to the multineutron energy.

-Only a lower limit of the mass number A can be inferred from the transformation into b ○. -A contaminant different from a ○ could lead to b ○ without requiring a multineutron.

-The generally uncontrolled previous reaction may produce a huge background of many particles, that could be eventually responsible for the production of b ○.

Concerning the mass number of the multineutron candidate systems, although 3 n is the simplest one, the pairing effects observed along the neutron dripline suggest that 4,6,8 n could more likely exhibit bound or resonant states. However, the availability of beams and targets, and the significance of the cross-sections involved, have mostly limited the search to 3,4 n. This is particularly true in the missing-mass technique, since the charged partners define the multineutron mass number and A > 4 would require too unbalanced partners, and/or negligible crosssections. On the other hand, the ambiguity in mass number of two-step reactions has opened the (uncertain) way to the search for heavier multineutrons, as we will see.

Therefore, due to accessibility and pairing considerations, the tetraneutron has been the most natural aim of these searches. Concerning its eventual binding energy, an upper limit is provided by the minimum of the experimental 4n separation energy (S 4n ) of bound nuclei, since a higher tetraneutron binding would make them unbound. For most of the quest this value was S 4n ( 8 He) = 3.1 MeV, but recent mass measurements have decreased it to S 4n ( 19 B) = 1.5(4) MeV [18]. However, if the tetraneutron were bound by more than 1 MeV, α+ 4 n would be the first particle threshold in 8 He. As the breakup of 8 He is dominated by the 6 He channel [19], the tetraneutron, if bound, should be so by less than 1 MeV.

Before starting the detailed review of the experiments, note that we will only address the results that have been published in accessible international journals. Moreover, we will not cover experiments that tried to probe the existence of multineutrons indirectly through the search of analog states in neighboring nuclei, like 3,4 H and 3,4 He.

The pion probe

One of the cleanest probes is the double charge exchange (DCX) reaction (π -, π + ), which falls within the missingmass category. The incoming negative pion becomes positive, changing two protons in the target into neutrons. By measuring the momenta of the pions on a 3,4 He target, one can measure the missing mass of the 3,4 n system. Obviously, since no other helium isotopes can be used as targets, no other multineutrons are accessible.

To our knowledge, this axis started in 1965 when Gilly et al. searched for the 4 He(π -, π + ) 4 n reaction [20]. They found no evidence of the tetraneutron, that would manifest itself through peaks in the π + spectrum due to twobody kinematics. The spectrum could be explained by introducing a final-state interaction (FSI) between two of the neutrons, as can be seen in Fig. 2 (red curve).

In 1970 Sperinde et al. used a 3 He target and searched for the 3 He(π -, π + ) 3 n reaction [21]. No evidence for a bound trineutron was found, but an enhancement at low 3n energies was observed. The same enhancement was observed in 1974 [22], but it was finally explained by the FSI between two of the neutrons or between all of them in a two-step DCX [23].

In 1976 Bistirlich et al. used a variant of this technique by studying the single charge exchange (SCX) on hydrogen, 3 H(π -, γ)3n [24]. After the subtraction of the different contaminants present in the tritium cell, no evidence for a trineutron was found in the γ spectrum. In [n+n+n+n] Fig. 2. Experimental results for the reaction 4 He(π -, π + )4n at θ = 0 • . The curves represent the phase space for the final states n + n + n + n, (2n) + n + n, and (2n) + (2n), normalized at 232 MeV and including the experimental conditions. The small axis represents the energy of the neutrons around the 4n threshold. Adapted from Ref. [20].

1980, Miller et al. undertook the same study with better statistics and resolution [25], and confirmed that no trineutron was needed to describe the observed spectrum, only the 2n FSI was required.

In 1984 Ungar et al. revisited the tetraneutron with the 4 He(π -, π + ) 4 n reaction [26]. As a test, they also measured the 12 C(π -, π + ) 12 Be reaction, in which the formation of states in 12 Be appeared clearly as a peak in the π + spectrum (Fig. 3, triangles). In the 4n channel, however, no apparent peak was observed (Fig. 3, circles). A few events laid in the expected region of bound 4 n, although at a rate similar to the events beyond that region, which are kinematically forbidden. The latter background was due to imperfect rejection of the events in which the π + decayed inside the spectrometer. In search of a possible resonant behavior, the spectrum below the 4n threshold was found to exhibit a broad enhancement with respect to four-neutron phase space, but it was consistent with the FSI of two neutron pairs [26].

In 1986, Stetz et al. extended the search to the trineutron with the reactions 3,4 He(π -, π + ) 3,4 n [27], at several energies and angles, but again no evidence of trineutron or tetraneutron states was found.

In 1989 Gorringe et al. repeated Ungar's experiment at lower energy, to reduce the continuum DCX contribution, but at a higher angle, corresponding to a higher momentum transfer to the neutrons [28]. After background subtraction, 6 counts remained in the possible bound 4 n window, but those were consistent with the estimated continuum contribution. Only an upper limit of the production cross-section could be deduced.

In 1997, a systematic study of the 3 He(π -, π + )3n reaction by Yuly et al. [29] at 120-240 MeV did not show 4 He(π -,π + ) 4n 12 C(π -,π + ) 12 Be (÷2)

T(π -) = 165 MeV Fig. 3. Experimental results for the reactions 4 He(π -, π + )4n at θ = 0 • (circles) and 12 C(π -, π + ) 12 Be at θ = 8 • (triangles, divided by 2). The curves are fits to guide the eye, with a Woods-Saxon distribution only (red) plus an additional Gaussian function (blue). The peak in the 12 C channel corresponds to the formation of the 12 Be ground and two first excited states, and the range in yellow in the 4 He channel to the region expected for a bound tetraneutron. Adapted from Ref. [26].

evidence for a trineutron, and indicated that the DCX process proceeds as a sequential SCX one. In 1999, Gräter et al. studied the same reaction at lower energies, 65-120 MeV, and found no evidence either [30].

For the sake of completeness, let us mention two special, non-conventional uses of the negative-pion probe. Earlier, in 1979, Chultem et al. had proposed an original use of a π -beam in a two-step process [31]. They hoped to produce a bound tetraneutron in the 208 Pb(π -, π + ) 4 n reaction, by DCX on an α cluster inside lead, and then measure the tetraneutron absorption by another lead nucleus, leading to 212 Pb. However, they found no α particles from its decay chain into 212 Bi and then 212 Po. Finally, in 1991 Gornov et al. studied the 3n missing mass in the reactions 9 Be(π -, t 3 He) and 9 Be(π -, d 4 He) [32]. The missing-mass spectrum was tentatively described using a very broad trineutron resonance at 3 MeV, but the very limited resolution did not allow to draw firm conclusions.

In summary, after more than 30 years of DCX experiments with pion beams, only upper limits following the non-observation of multineutrons have been set, and the technique is not being presently used. However, the DCX technique has been recently revisited by Kisamori et al. using exotic nuclei [4], as we will see in Sec. 2.5.

Fission and other activation probes

The DCX pion reactions are a very clean and powerful probe, although the cross-sections involved are extremely low. From the start of the multineutron program a second avenue was taken, exploiting the two-step principle described previously. In the first step, the multineutron is supposed to be produced in a high-flux reaction, not necessarily well-characterized, like spallation or induced fission. Assuming that a bound multineutron was thus produced, together with many other particles, we may use it in a secondary two-body reaction that transforms a given sample in a unique way. Therefore, we need to demonstrate that the sample was transformed, but above all that it could not be transformed in any other alternative way.

Already in 1963, Schiffer et al. irradiated samples of nitrogen and aluminum in a nuclear reactor, searching for the reactions 14 N( 4 n, n) 17 N and 27 Al( 4 n, t) 28 Mg [33]. If bound tetraneutrons had been produced in the fission of the uranium fuel, they might have observed the β-delayed neutron decay of 17 N and the γ-rays from 28 Mg and its daughter 28 Al. However, none of them were significantly observed above the background.

In 1965 Cierjacks et al. carried out a similar experiment, with samples of nitrogen, oxygen, magnesium and other heavier elements surrounding a uranium target bombarded with deuterons [34]. For the three lighter samples, they searched for bound tetraneutrons emitted in uranium fission through the reactions 14 N( 4 n, n) 17 N, 16 O( 4 n, t) 17 N and 26 Mg( 4 n, 2n) 28 Mg. However, they were not able to observe the corresponding β-delayed neutrons or γ-rays above the background.

In 1968 Fujikawa et al. searched for bound trineutrons in the reaction n+ 7 Li, that would transform a barium sample through the secondary reaction 138 Ba( 3 n, n) 140 Ba [35]. Barium was in the form of a high-purity BaCO 3 powder, and the potential production of 140 Ba would have been identified through the γ-rays of its daughter, 140 La. After chemical separation of lanthanum from the powder, no γ-rays were observed, and an upper limit of the crosssection was provided.

In 1977 the multineutron search seemed to come to a happy end at CERN. Détraz obtained the first positive result after exposing a natural zinc sample to a tungsten block irradiated with a proton beam of 24 GeV [36] (Fig. 4). Block and sample were separated by an aluminum screen, supposed to be thick enough to stop any charged particle that could have been produced, and the reaction that was sought was 64,66,67,68,70 Zn( A n, xn) 72 Zn, i.e. the transformation of a natural zinc isotope into 72 Zn through the interaction with a neutral cluster. The relatively long half-life of 72 Zn (46.5 h) and its daughter 72 Ga (14.1 h), and the several γ-rays emitted by the latter, made it possible to remove the sample from the high-activity area and perform a clean measurement of those γ-rays.

After several-day exposure, the 100 g of natural zinc were sent by air freight to Orsay, where a Ge(Li) detector revealed the unambiguous production of 72 Zn through five prominent γ-rays of its daughter 72 Ga. This observation required up to four chemical separations of gallium, [36]. In the collision of protons and tungsten a neutral cluster was supposed to be produced, go through a screen, and induce ( A n, xn) reactions on natural zinc. The new 72 Zn isotope was detected through its decay into 72 Ga * with a germanium detector.

due to the high activity induced in the sample by other beam-induced reactions. The production of 72 Zn could have been the result of any (A > 1) bound multineutron, but since the previous searches for A = 2-4 had failed, Détraz concluded that most likely bound hexa-or octaneutrons (A = 6, 8) had been produced [36].

Only a few months later, Turkevich et al. tried to confirm Détraz's exciting results using a similar principle [37]. As the first step, they irradiated uranium with 700 MeV protons, and exposed a lead sample behind an aluminum screen, again supposed to stop any charged particle. For the second step, they used the one that was described in the previous section for Chultem [31], the search for the 208 Pb( A n, xn) 212 Pb reaction through the α particles from its decay into 212 Bi and 212 Po. No trace of "polyneutrons" was found, shedding doubt on Détraz's results.

In 1980, de Boer et al. followed suit. Since Turkevich concluded that their results would have not been sensitive to multineutrons with A < 6 [37], a possible solution to Détraz's puzzle was that he had observed only 4 n, and not 6,8 n as he had suggested [36]. Therefore, de Boer irradiated a block of tellurium with 3 He ions and searched for the interaction of bound tetraneutrons in the same block through the 130 Te( 4 n, 2n) 132 Te reaction [38]. No γ-rays from the decay of 132 Te were observed, and de Boer concluded, in agreement with Détraz, that the latter had in fact underestimated the transmission of tritons through his aluminum shield and that the reaction observed had been (nat) Zn(t, p) 72 Zn, without any multineutron [38].

This tree-year period of excitement ended up with some frustration, but the multineutron quest went on. However, following these results, the two-step or activation probe was abandoned for more than thirty years in favor of the cleaner missing-mass experiments.

Anyhow, in 2012 Novatsky et al. revisited the technique by inducing uranium fission with 62 MeV α particles, and then searched for the interaction of bound multineutrons within strontium and aluminum samples shielded by thin Kapton and beryllium layers [39,40]. The authors claimed they had observed the reaction 88 Sr( A n, xn) 92 Sr, through γ-rays from the daughter 92 Y [39], and the reaction 27 Al( A n, p xn) 28 Mg, through γ-rays from the decay into 28 Al and 28 Si [40]. They concluded that, due to the lack of evidence for a bound tetraneutron, their results should "certainly" correspond to heavier clusters (A 6). However, the experimental hall exhibited strong (and diverse) activity, the number of reaction channels following fission is extremely high, and it could be that the screens that were used were not thick enough to shield all the many charged particles produced, as it had been the case for Détraz in 1977.

To conclude, the so-called activation probe appeared as a powerful tool from the very first years, lead to some hope about bound multineutrons in the late 70s, but after the refutation of the claim it was left out and is not being followed nowadays in any large-scale facility.

The multinucleon-transfer probe

A third research axis was also opened in the early years. Like the pion DCX, it relied on the cleaner missing-mass principle. However, instead of changing two protons into two neutrons, which limited the search to 3,4 He targets, the multineutron system was probed in a transfer reaction between two stable nuclei, increasing the possible combinations. But as in the DCX reaction, the cross-sections were expected to be very low, since one needed to: transfer all the protons away from one nucleus; sometimes bring some neutrons back in the opposite sense; and in any event lead exclusively to a two-body final state.

In 1965 Ajdačić et al. started to use this probe with a simple (n, p) transfer reaction on the triton [41]. With a 14 MeV neutron beam they detected protons from the 3 H(n, p)3n reaction. Some events were observed at a missing mass of about 1 MeV below the 3n threshold, the first candidates of a (quite) bound trineutron. However, knowing that the very first experiments had failed to find a bound tetraneutron, more likely to exist due to pairing, they concluded that their result was "highly improbable". One year later, Thornton et al. repeated the same experiment with 21 MeV neutrons and better resolution [42], and found no evidence for a bound trineutron.

In 1968 Ohlsen et al. used a triton beam and searched already for a more complex transfer reaction, 3 H(t, 3 He)3n [43]. The missing mass reconstructed from 3 He lead to a deviation from four-body phase space, only at forward angles, that could be consistent with a low-energy trineutron resonance. They were not able to exclude, however, an effect from the reaction mechanism itself.

In 1974 Cerny et al. started to use heavier nuclei and searched for the tri-and tetraneutron in the reactions 7 Li( 7 Li, 11 C)3n and 7 Li( 7 Li, 10 C)4n [44]. Concerning the trineutron, the intense 11 C channel led to a very high statistics spectrum. Unlike the 7 Li( 7 Li, 11 B) 3 H channel, that exhibited several structures, the 3n missing-mass spectrum could be well described by four-body phase space, plus some small peaks from known target contaminants (that lead to 11 C partners different from 3n). 7Li, 10 C)4n reaction at θ = 7.4 • . Known contaminant reactions are indicated either explicitly or with an asterisk. The red curve corresponds to five-body phase space, and the range in yellow to the region expected for a bound tetraneutron. Adapted from Ref. [44].

In the tetraneutron channel, however, the low 10 C production led to a poor separation from the tail of the much stronger 11 C distribution. The resulting 4n missing-mass spectrum could be described by five-body phase space plus the known contaminants, as can be seen in Fig. 5. Although some events are visible in the possible region for bound 4 n, the signal was not significant with respect to the background level. They concluded that the purity was still much worse than in the DCX experiments, but that by using different beam-target combinations the way had been opened to heavier multineutrons.

In 1988 Belozyorov et al. improved on the main issues of Cerny's work, the target purity and fragment identification, and searched for the trineutron in the 7 Li( 11 B, 15 O)3n reaction and for the tetraneutron in the 7 Li( 11 B, 14 O)4n, 7 Li( 9 Be, 12 N)4n and 9 Be( 9 Be, 14 O)4n reactions [45]. In all the reactions on 7 Li, the missing-mass spectra above the 3n and 4n thresholds were well described by the corresponding four-and five-body phase space, showing no evidence for multineutron resonances. Although the two reactions leading to 14 O showed a few events below the 4n threshold, they were consistent with the background due to pulse pileup or to beryllium target impurities [45].

In 1995 Bohlen et al. used a 14 C beam in order to probe very neutron-rich missing masses, among them the trineutron in the reaction 2 H( 14 C, 13 N)3n on a CD 2 target [46]. The missing mass was fully described below the threshold by the carbon contribution and above it by the decay of a broad 15 N resonance. Finally, in 2005 Aleksandrov et al. repeated Cerny's experiment, with similar beam energy and target, and obtained the same negative results for both the tri-and tetraneutron [47].

This technique appeared as a good compromise between pion DCX and activation. It could access heavy multineutrons from a variety of beam-target combinations, and the potential signals were supposed to be unambiguous. However, the absence in practice of clear signals led towards a need for higher purities, and the technique has been put aside for the last fifteen years.

The GANIL 2002 result

The experiments performed in the XX century used mainly stable beams and targets. The beams could thus be very intense, but building a neutral system from balanced combinations of protons and neutrons required reactions with very low cross-sections. Moreover, the potential multineutron signal often shared parts of the spectra with background from contaminant species, and due to the low cross-sections used the background contributions became too important for a signal to be clearly established.

In 2002 Marqués et al. proposed at the GANIL facility a new technique that could solve those issues [START_REF] Marqués | [END_REF]. With the advent of radioactive secondary beams, the possible preformation of multineutrons inside very neutron-rich nuclei was considered, similar to the preformation of α particles in the process of α decay. Within this scenario, the until then complex formation step of multineutrons could be reduced to the breakup of one of those nuclei, with an increase in cross-section of several orders of magnitude (mb, compared to the nb or pb of the previous probes) due to the weak binding of these clusters.

The radioactive beam was 14 Be, in which the 10 Be+4n threshold is at only 5 MeV. Following the breakup on a carbon target at 35 MeV/N, the detection of the 10 Be fragment provided a clean signature of the channel. For the detection of a potentially liberated tetraneutron cluster, the principle was similar to that used by Chadwick in the discovery of the neutron [48]: deduce the mass of the neutral particle from the recoil induced by elastic scattering on charged particles. The recoil energy E p of a proton in the organic scintillator detector is related to the energy per nucleon E n of the incoming A n cluster, obtained from its time of flight: (E p /E n ) 4A 2 /(A + 1) 2 . Since the dineutron is unbound, and the trineutron should also be due to pairing, the measurement of proton recoils over the incoming neutron energy of 1.5-2.5 could only be attributed to a bound tetraneutron.

The method was also applied to data from 11 Li and 15 B beams, but only in the case of 14 Be some events were observed with characteristics consistent with the production and detection of a multineutron cluster. The 6 events are shown in Fig. 6 in yellow, with proton recoils 1.4-2.2 times higher than those expected for individual neutrons, and appear all in coincidence with the detection of a 10 Be fragment (note the absence of events in the most abundant 12 Be channel). Special care was taken to estimate the effects of pileup, i.e. the detection of several neutrons in the same detector, and it was found that it could at most account for some 10% of the observed signal. The conclusion for the most probable scenario was the formation of a bound tetraneutron in coincidence with 10 Be [START_REF] Marqués | [END_REF].

As we will see in Sec. 3, this result triggered several theoretical calculations that could not explain the possible

(E p /E n ) max = f(A) 1 2 4
Fig. 6. Scatter plot of the particle identification parameter PID vs the proton recoil in the neutron detector (normalized to the neutron energy) for the reaction ( 14 Be, X + n). The dotted lines show the region centered on the 10 Be peak and with Ep/En > 1.4, and the 6 events in yellow are candidates to the formation of a bound tetraneutron. The scale on the upper axis shows the maximum proton recoil as a function of the multineutron mass number. Adapted from Ref. [START_REF] Marqués | [END_REF].

binding of the tetraneutron. Moreover, another work questioned the probe itself [49], arguing that a weakly bound tetraneutron would rather undergo breakup than elastic scattering. Marqués et al. addressed both issues [50], finding that the signal observed could be generated also by a low-energy tetraneutron resonance (E 2 MeV) through an enhancement of pileup, or by the breakup of a bound tetraneutron in the scintillator followed by the detection of some of the neutrons. Interestingly, the latter scenario had been already suggested by Brill et al. in 1964 as a potential probe of bound multineutrons, but no abnormal signals were observed at that time in the bombardment of a 48 Ca target with 12 C and 3 He beams [51]. On the other hand, another experiment in 1971 by Koral et al. had already tried to detect bound trineutrons from abnormal induced recoils [52]. Following the bombardment of a lithium target with neutrons similar to the one in Ref. [35], they searched abnormal recoils of 4 He ions by comparing the time of flight and the pulse height signals in a helium scintillator, the same principle used at GANIL, but no evidence was found.

The decrease in beam intensities at GANIL and the aging of the neutron detector used did not allow an unambiguous confirmation of this result. Other experiments at GANIL using the missing-mass technique were not able to find positive signals either. After some theoretical works in the early 2000s, the field became quiet. The dashed (blue) curve represents only the latter, multiplied by a factor of 10 in order to make it visible. The 4 events at threshold are highlighted in yellow. Adapted from Ref. [4].

The RIKEN 2016 result

In 2016, Kisamori et al. proposed a new probe at RIKEN: 4 He( 8 He, 8 Be) 4 n, a DCX reaction using exotic nuclei [4]. Sending a very intense 8 He beam at 186 MeV/N onto a liquid 4 He target, the exit channel was selected through the detection in an spectrometer of the two α particles from the decay in flight of 8 Be. The large Q value of the ( 8 He, 8 Be) reaction almost compensated the binding energy of the α particle and allowed for the formation of a 4n system with small momentum transfer. The authors expected in this way to enhance the odds of a weakly interacting tetraneutron in the final state, that would remain in the target area (and would not be directly detected). The 4n missing-mass spectrum (Fig. 7) showed 4 events very close to threshold. The relative energy and angle between the two α particles were consistent with the formation of the 8 Be ground state. The background was estimated from the probability of having two 8 He beam particles from the same bunch breaking up, due to the high beam intensity, and leading to the detection of two independent α particles at small angle that could mimic a 8 Be decay. It was found to be uniformly spread over the whole range of missing mass (blue curve in Fig. 7), with an integrated value of about 2 events [4].

The red curve in Fig. 7 corresponds to a calculation of the direct decay of the 4n final state within a wave packet similar to the initial 4 He, including the interaction between neutrons and between neutron pairs [4]. It includes also the estimated background described above, not visible at the scale of the figure. This curve, without the hypothesis of a tetraneutron resonance state, clearly cannot explain the events observed around the 4n threshold. Note that there was only 1 event in the kinematically forbidden region (at about -20 MeV in Fig. 7), an additional indication of the low background.

The 4 events in the region 0 < E 4n < 2 MeV were found to be consistent with the formation of a tetraneutron resonance at E( 4 n) = 0.8 ± 1.3 MeV, with a width Γ < 2.6 MeV, taking into account both the statistic and systematic errors (the latter were dominant due to the missing-mass calibration procedure). It should be noted that, because of the missing-mass uncertainty, this result is also consistent with the formation of a bound tetraneutron. The cross-section corresponding to the candidate events was estimated at about 4 nb, compatible with a simple estimation of the two-step process assuming the Gamow-Teller and the spin-dipole transition for the projectile and the target system, respectively [4].

The neutron probe

As mentioned when introducing the experimental techniques, paradoxically almost none of the many experiments already reviewed has tried to detect the neutrons. With the recent improvements in beam intensities and neutron detection efficiencies in the world leading facilities, trying to detect the neutrons from the decay of a multineutron system seems a logical next step, which in addition would give access to eventual neutron correlations in their decay:

a ○------→ X ⇓ ⇑ c ○------→ ( A n) -→ A n ○ neutron detection
For the sake of completeness, we note that in 2016 Bystritsky et al. claimed the first direct detection of multineutrons with A = 6, 8 [53]. They used an array of 20 3 He counters, with a total estimated efficiency of ε 1n ∼ 20%, surrounding a 238 U sample, and searched for the neutrons that would stem from the multineutron cluster decay of uranium. Although they based their claim in the detection of few events with 5 neutron hits, they also concluded that the confirmation would require an increase of the statistics, and a reduction of the cosmic and natural backgrounds, by at least one order of magnitude [53].

The Radioactive Isotope Beam Factory (RIBF) of the RIKEN Nishina Center provides nowadays the highest intensities of light neutron-rich beams, together with high neutron-detection efficiencies. Since the 2016 result, two experiments have been undertaken at RIKEN aiming at the detection, for the first time, of all the neutrons emitted in a multineutron (x > 2) decay. They benefited both from the combined capabilities of the NEBULA [54] and Neu-LAND demonstrator [55] neutron detector arrays. With an average neutron efficiency of ε 1n ∼ 45%, the granularity of the arrays allowed for an unprecedented fourneutron efficiency of ε 4n ∼ 1-2% (see Fig. 1).

The first experiment searched for the ground state of 28 O, one of the grails of nuclear structure physics, using . The solid histogram represents the number of experiments reviewed that were searching for multineutrons (36 total, mainly for tri-and tetraneutrons) as a function of the year of publication of the results. In yellow are those searching specifically for the tetraneutron. The stars represent the three positive signals reported, the empty one that was refuted [36] and the two solid ones that have not been contested yet [START_REF] Marqués | [END_REF]4]. The pale Gaussians guide the eye through the recurring pattern.

the reaction H( 29 F, 24 O)4n [START_REF] Kondo | RIBF Proposal NP1312-SAMURAI21[END_REF]. The moderate intensity of the very exotic 29 F beam lead to some 100 complete 4n events. Although the analysis of the results is still in progress, it seems that the decay of the ground state of 28 O would proceed through the narrow 26 O ground state, i.e. it would be a sequential 2n-2n decay. As such, no information should be derived concerning multineutrons.

The second experiment used a similar setup and technique, the proton removal from a beam and the detection of a fragment plus four neutrons, but this time in search of the 7 H ground-state energy and its potential tetraneutron decay with the reaction H( 8 He, 3 H)4n [START_REF] Marqués | RIBF Proposal NP1512[END_REF]. With respect to the previous experiment, the much higher intensity of the 8 He beam will lead to several orders of magnitude increase in statistics, providing tens of thousands of 4n events, while the absence of low-lying 4,5,6 H resonances should allow for the observation of a direct 4n decay, and the eventual correlations within.

A sixty-year quest

In this section we have reviewed experiments taking place over a sixty-year period, all searching for evidence of multineutron existence. Their chronology (Fig. 8) exhibits several trends. Even if the techniques have been diverse and their sensitivity has increased with time, we can see a recurring pattern of 'bunches' during the first forty years, with experiments accumulating from mid to end of each decade. Some experiments at mid-decade triggered others, and then the overall negative results lead to a stop in the program, until someone else restarted it a few years later. Towards the end of the century the number of experiments in each bunch decreased, showing signs of exhaustion due to the lack of positive signals.

The number of experiments in the present century has been much lower, although two positive signals of a bound or low-lying resonant tetraneutron were obtained. As said in the introduction those signals have renewed the interest in the field, both experimentally and, and as we will see in the next section, theoretically. Therefore, for the next extension of Fig. 8 we expect an upcoming significant 'bunch' of results. Taking into account the increasing accuracy and sensitivity of these new experiments, the next few years will possibly see the end of this quest whatever the outcome, at least with respect to 3,4 n.

Moreover, besides the aforementioned experiments aiming at the detection of four neutrons already carried out [START_REF] Kondo | RIBF Proposal NP1312-SAMURAI21[END_REF][START_REF] Marqués | RIBF Proposal NP1512[END_REF], complementary missing-mass experiments are also being programmed at RIKEN, without neutron detection but with increased sensitivity. The tetraneutron has already been revisited using the same reaction as in 2016, 4 He( 8 He, 8 Be) 4 n, with several improvements in the experimental conditions [START_REF] Shimoura | RIBF Proposal NP1512-SHARAQ10[END_REF]. It has also been probed in the knockout of an α particle off 8 He at backward angles in quasi-free conditions, H( 8 He, pα) 4 n [START_REF] Paschalis | RIBF Proposal NP1406-SAMURAI19[END_REF], avoiding the FSI of the 4n with the other particles in the reaction. There are also new plans to probe the trineutron, with the reaction 3 H(t, 3 He) 3 n [START_REF] Miki | RIBF Proposal NP1712-SHARAQ11[END_REF].

Among those missing-mass experiments, some are already searching for the next heavier system, the hexaneutron. One has already been carried out, knocking out two α particles from 14 Be in the reactions H( 12,14 Be, pαα) 4,6 n [START_REF] Beaumel | RIBF Proposal NP1206-SAMURAI12[END_REF], and a second one is planned in a near future, knocking out an α particle and a proton from 11 Li in the reaction H( 11 Li, ppα) 6 n [START_REF] Nakamura | RIBF Proposal NP1812[END_REF]. Depending on their results, future experiments could be planned in order to study neutron correlations in the decay of hexaneutron states.

However, the hexaneutron seems to represent a mass frontier difficult to cross in the laboratory. In the next section we will see how theoretical calculations have been dealing with the lightest multineutrons, and how they could help us go beyond A = 6 in order to understand possible binding energy trends in these systems.

Theoretical calculations

The theoretical interest in bound or resonant multineutron states is as old as nuclear theory, and its development goes in parallel with the technical possibility to obtain accurate quantum mechanical solutions of few interacting particles using reliable interactions. This is an essential ingredient when dealing with systems very close to, or above, their dissociation threshold. The progress in this domain has been slow. The first rigorous formulation of the three-body problem in Quantum Mechanics dates from 1960 [63], and the first realistic solution took still several years [64]. As it has been the case in the experiments, the main efforts of the theory have been concentrated on 3,4 n, and to a lesser extent on the next heavier systems, 6,8 n.

Even if we can trace the origin of both theoretical and experimental programs back to the 1960s, it seems clear that the theoretical community has significantly increased its interest in 3n and 4n systems following each of the two experimental signals described in the previous section. We will thus review the many theoretical works from that perspective. We will describe the activity prior to GANIL results on the breakup of 14 Be into 10 Be+4n [START_REF] Marqués | [END_REF], then review the different approaches that were triggered by this result, and finally discuss the many calculations that followed the second signal in the DCX reaction 4 He( 8 He, 8 Be)4n [4]. But before moving to the specific works, let us start with some general considerations about the theoretical treatment of the few-neutron problem.

The main difficulties

The theoretical study of A n systems has faced two hardwearing difficulties. The first one is finding a reliable method to solve exactly the A-body problem (A = 3, 4...) for loosely-bound states as well as for resonant states, which could be embedded in the continuum and very far from the physical region. By "exactly" we mean a rigorous mathematical formulation of the quantum mechanical problem and a uniformly convergent numerical scheme which could lead to the required accuracy. This implies an ab initio solution, i.e. only in terms of the nn interaction. We anticipate that for such diluted systems the three-nucleon forces (3NF) should play a negligible role.

Since such systems, if they manifest at all in nature, would be a subtle balance between the attractive nn interaction and the Pauli 'repulsion' and live in the continuum, any uncontrolled approximation, in the antisymmetry or in the boundary conditions, could be fatal and lead to illusionary conclusions. This constraint delayed the first serious attempt to study the 3n system to the end of the 60s, by Mitra et al. [65], and this was still based on an ad-hoc method they had developed for solving the Schrodinger equation with a rank-1 separable interaction, which furthermore turned to be not very realistic.

The second difficulty is related to the fact that, until the contrary is proven, A n has no bound subsystems. Therefore, the usual scattering theory with two-cluster asymptotic channels cannot be applied to obtain the scattering amplitude and eventually the location of the corresponding poles. The only pertinent quantity, the amplitude for the 3 → 3 or 4 → 4 scattering processes at positive energy (with respect to the corresponding thresholds), is well defined but very difficult to compute, even for the simplest cases. Moreover, it is not clear that it could provide some information about possible resonant states when they are far enough from the real axis.

Last but not least, we would like to point out that the very concept of "resonance" is far from being clear, and in any event it is certainly not shared by the physics community (experimentalist and theorists from many fields), even at the two-body level. As we will show, it is possible to generate extremely sharp 'resonant' cross-sections without any dynamical near-threshold state counterpart. Therefore, associating a 'state' to a sharp structure in the cross-section could be a dangerous shortcut. Moreover, the quite widely adopted definition of resonance as a rapid increase of a phase-shift crossing π/2 is not very satisfactory. First, it is not clear what "rapid" means and, most importantly, a slight modification in the interaction could significantly lower the value of the phase-shift and make a resonance 'disappear', a bit surprising when dealing with a physical state. Not to talk about all the complexities appearing in a coupled-channel problem.

With all that in mind, it seems that the only unambiguous approach to conclude about the existence of a A n resonance is determining the location of the relevant S-matrix pole, or any equivalent quantity like the A-neutron Green function, in the complex-energy manifold (or complex-momentum plane). This allows to define the intrinsic parameters of the resonance in the usual way:

E = E R -i Γ 2 (1)
whatever its proximity to the physical axis and independently of any phase-shift it can create anywhere.

To this aim, it could be worth reminding the location of the singularities in the complex-momentum plane and in the corresponding energy Riemann surface. This is sketched in Fig. 9 for the simplest case of a two-cluster problem. In the left panel a bound (

k 0 = +i | k 0 |) and a virtual (k 1 = -i | k 1 |
) state appear as singularities in the imaginary k-axis, in the upper (bound) and lower (virtual) half planes. They are mapped onto the two Riemann sheets E I and E II of the E-surface, in the right panel, respectively at the points E 0 ∈ E I and E 1 ∈ E II . Note that the energies of the virtual states are negative, as are for bound states, despite being in the continuum. The upper half k-plane is mapped onto the first ("physical") Riemann sheet E I , and the lower one onto the second ("unphysical") Riemann sheet E II , which are glued to each other by the positive real axis E > 0.

Resonances manifest as poles in the lower half k-plane, k 2 = k R -ik I with k I > 0, and are mapped onto E 2 in the unphysical sheet E II . They appear always by pairs (k 2 , k ′ 2 = -k * 2 ) symmetric with respect to the imaginary k-axis. Note that Re(E 2 ) > 0 only if Im(k 2 ) < Re(k 2 ). The situation is more complicated in presence of open channels with open thresholds as well as dynamical cuts [66], but what is described above is enough for our purposes. The interested reader can find a more detailed description in the standard books of scattering theory [START_REF] Baz | Scattering, reactions and decay in nonrelativistic Quantum Mechanics[END_REF][START_REF] Taylor | Scattering Theory[END_REF][START_REF] Goldberger | Collision Theory[END_REF][START_REF] Newton | Scattering Theory of Waves and Particles[END_REF][START_REF] Kukulin | Theory of Resonances: principles and applications[END_REF].

Whether or not this A n singularity could manifest in any observable, and in which way, must be considered as a different (and difficult) problem, especially taking into account that such states can be produced in complex nuclear reactions that could mask, generally in an unknown way, the intrinsic parameters (E R , Γ ) of an eventual resonance.

Multineutron states in the continuum

The research for an eventual bound state is relatively easy. If a direct bound state calculation gives a negative eigenstate of the Hamiltonian, the system is bound. On the contrary, we conclude that there is no bound state... within the accuracy of the method, which provides then a lower limit of the binding energy. The method for computing bound states must be accurate enough to deal with loosely bound systems, i.e. quite extended objects in configuration space. We would like to point out from now that no serious calculation has ever predicted a 3n or 4n bound state.

Γ k 1 E E 0 2 E 1 0 k k 2 E=k E=k 2 2 k-plane E-surface E I E II γ 1 1 Γ 2 γ 2
However, it is worth noting that the aforementioned pioneering calculation of Mitra et al. [65] was quite optimistic. They realized that a pure S-wave nn interaction generates a strong repulsion in the 3n system due to the Pauli principle, but concluded that "a moderately attractive 3 P force ... is enough to produce a bound 3 n state", and added that the "range and strength of the force seem to be in good accord with the 3 P 0 phase-shift data of Bryan and Scott", one of the most realistic potentials available at that time. Their conclusion, certainly due to the simplicity of their dynamical approach and interactions, would be denied very soon.

In the absence of bound states, the optimal choice is to find the location of the nearest singularity in the continuum, i.e. the parameters (E R , Γ ) of a virtual or resonant state with all the caveats mentioned. This is however a difficult task when working in configuration space, mainly due to the exponentially increasing boundary conditions. The 'regalian' way to avoid the latter is the Complex Scaling Method (CSM) [START_REF] Nuttal | [END_REF]73,74]. In the two-body case it consists of an analytic continuation of the Schrödinger equation by means of a complex rotation of the coordinates (r → r e iθ ), which transforms the exponentially increasing boundary conditions of a resonance into square integrable ones, provided the "rotation angle" θ is larger than half the argument θ r of the complex resonance energy (1):

0 < θ - θ r 2 < π 2 tan θ r = Γ 2E R
Under these conditions, illustrated in Fig. 10, the resonant state can be computed by means of the standard bound-state problem techniques. This approach, supplemented with many variants, has been successfully extended to A = 3, 4, 5 [START_REF] Nuttal | [END_REF]75,76,77,78,79]. Note, however, that the numerical difficulty increases dramatically with the rota- tion angle, in turn limited by the analytical cuts of the two-body interaction, and this restricts the practical application of CSM to relatively narrow states [77].

To circumvent this problem, one possibility is to introduce an additional binding to the physical system under study until it gets artificially bound. One can then remove adiabatically the 'artifact' and follow the evolution of the system in the continuum until it reaches the physical state. A commonly used artifact consists of a "scaling factor" s in front of the nn interaction, which can eventually depend on each partial wave α ≡ {L, S, J}:

V α sα (r) = s α V α nn (r) (2)
and increase s until the multineutron system binds. Note that only moderate s values are required to bind the dineutron with all realistic and semi-realistic V nn models, as shown in Fig. 11. For instance, the critical value to get B nn > 0 is s = 1.080 using the Argonne AV18 potential [80] and s = 1.087 using Nijmegen Reid93 [81]. For the phenomenological CD-MT13 [82,83], it is only slightly larger, s = 1.101, probably due to its shorter range. Another possibility consists in placing the multineutron system inside an external potential, which acts as a confining trap. This is implemented by adding to the total potential energy V of the system a term of the form:

V = A i<j=1 V ij → V = A i<j=1 V ij + A i=1 V T (r i )
where V T can be an harmonic oscillator (HO) or a Woods-Saxon potential. The quantum mechanical solution of 2n and 3n systems in a Woods-Saxon trap is treated with some detail in the Appendix.

The greatest advantage of this approach is the preparation of the system in an a priori well-controlled bound state. Moreover, the size of the extra energy required to bind the system constitutes already a first and valuable indication of how far from the physical region a resonant state may be. The main drawback, however, is that the artifact can dramatically change the properties of the system and induce some parasitic phenomena that have no physical reality. This unpleasant situation may happen if one makes use of the two following procedures. The first one is the introduction of spurious processes that can modify the evolution of the system. For example, if one artificially binds 3 n by introducing a scaling factor s in the nn interaction, the dineutron will be quickly bound by several tens of MeV. For s ≈ 3, a typical value required to bind 3 n, B nn ≈ 100 MeV. Therefore, a multineutron state built in such a way will contain one or several subsystems with very peculiar properties (e.g. extremely compact), that will certainly affect its evolution into the continuum. What one believes to be a bound 3 n state will be in fact a system strongly decaying into n+ 2 n.

If the method used to solve the three-body problem does not allow to properly account for open scattering channels, like when using basis of square integrable functions, this fact will not be detected in the numerical calculations but will influence the properties and the evolution of a state bound by brute force in an uncontrolled way. In any event, the evolution of the multineutron system when decreasing s will be totally dominated by this process and one may attribute to 3 n resonances what corresponds to an unphysical decay 3 n → n+ 2 n. The same will happen when studying the tetraneutron with the corresponding spurious channels 4 n → n+ 3 n and 4 n → 2 n+ 2 n.

The appearance of these spurious thresholds can be avoided by adding an A-neutron force involving only the total number of neutrons of the studied system. For example, for 3 n we can introduce a 3n force that we can choose for simplicity of hyper-radial form:

W 3n (ρ) = -W 0 f (ρ) ρ 2 = x 2 + y 2 (3)
where f (ρ) is an asymptotically decreasing function of ρ (like e -ρ or Yukawa) and (x, y) are the intrinsic Jacobi coordinates. This artifact does not modify any subsystem threshold and the evolution W 0 → 0 can be done safely.

A similar distortion is also induced when 'confining' neutrons in a Woods-Saxon trap, like for example: a common practice in the research on neutron droplets [14,15,16]. For small multineutrons, the evolution of their binding energies as a function of the potential strength V 0 has been used to draw conclusions about eventual 3n and 4n resonances. However, this approach presents some instabilities, as we illustrate for the simplest case of 3n in the upper panel of Fig. 12. If B 2n > B 3n for a given set (V 0 , R, a) of the confining potential (4), the 3 n will decay into a dimer, staying in the well, plus an outgoing neutron with kinetic energy proportional to B 2n -B 3n . If the methods used to compute the bound state account for the full quantum dynamics of the system, the supposed confined multineutron state must decay. This fact may strongly impact the extrapolation pattern. Since the presence of thresholds is associated with algebraic branching points, the open decay channels will generate sharp structures in the extrapolation function and the results of any extrapolation procedure become unreliable. The binding energy of multineutrons in a trap is in fact only fictitious, as are the properties extrapolated from their evolution as a function of the trap parameters. The same will happen for a 4 n state, with respect to the open thresholds n+ 3 n and 2n+ 2 n (Fig. 12, lower panel).

V T (r) = V 0 1 + e (r-R)/a (4)
The second mechanism that may induce parasitic phenomena is related to the evolution from a bound state to a state in the continuum. It is not a mere extrapolation from negative to positive values in an energy axis. As already noted, the quantum mechanical energy continuum is not an axis, even not a plane, but a Riemann manifold with dynamical cuts and delicate analytic properties.

The evolution of a bound state into the continuum can be properly described by the so-called Analytic Continuation on the Coupling Constant Method (ACCCM) [START_REF] Kukulin | Theory of Resonances: principles and applications[END_REF]. It consists in adding to the Hamiltonian H, for which we aim to determine the parameters of a resonance, a perturbation λW . We can then make use of the analytic dependence of the spectrum of the 'perturbed' Hamiltonian H λ = H + λW on the coupling constant λ. If k is the wave number of a bound state with negative energy E of H λ (E = 2 k 2 /2µ), its dependence on λ is given by the Padé-like analytic expansion:

k(λ) = i P N (z) Q M (z) = i c 0 + c 1 z + c 2 z 2 + . . . c N z N 1 + d 1 z + d 2 z 2 + . . . d M z M (5)
where P N and Q M are polynomials on the complex variable z = √ λ -λ c . The critical value λ c corresponds to a zero-energy bound state that moves into the continuum and to an algebraic branching point of the complex function (5).

ACCCM is schematically illustrated in Fig. 13 and proceeds as follows. In a first step, a series of bound state values (λ n , k n ) are computed, with λ > λ c and k n purely imaginary. They are inserted in the formal expansion (5) to determine the N + M unknown coefficients (c i , d i ), which are real. The k(λ) dependence being fixed, one can decrease λ → 0 to determine the resonance position of H in the complex plane. Note that for λ < λ c , z is a complex number, as will be k and the corresponding E. This method has been successfully used to compute the position of the nearest 3n and 4n resonances [84,85,86], as well as of the 4,5 H isotopes [87]. It requires an accurate determination of the bound state (λ n , k n ) parameters and a careful stability study as a function of the degrees N and M of the Padé expansion (5).

After this overview of the general issues related to the calculation of few-neutron systems, in the following we will review the most relevant theoretical works that have been published prior, in between and after the two experimental signals observed at GANIL in 2002 and RIKEN in 2016.

Previous to GANIL 2002

Before 2002 we find in the literature several works devoted to 3n and 4n systems which, although far from being uninteresting, can be considered as preliminary or incomplete from the physical point of view, either due to the quantum mechanical treatment of the three-body problem and/or to a very simplistic nn interaction. As we noted in the introduction, the first result to our knowledge is from Mitra et al. back to 1966 [65], using a technique for solving the three-body Schrödinger equation with rank-1 separable S-and P-wave nn interactions. Their optimistic conclusion about a bound 3 n was immediately contested by Okamoto et al. [88] and Barbi [89], showing with variational approaches that based on their knowledge of the NN interaction the trineutron was unbound by at least 10 MeV, and that the depth of the required nn potential well to bound it should be unphysically large.

The first 'exact' 3n calculation solving the Faddeev equations in momentum space, although using a rank-1 separable Yamaguchi force, is due to Glöckle in 1978 [90]. It is also the first time that the trajectory of the complexenergy 3 n resonance as a function of a scaling factor s was presented, together with a careful discussion of the analytic continuation properties of the scattering amplitude. The critical value for binding 3n was found to be s = 4.2, and the possibility of a low-energy resonance was ruled out. Of course the latter was on the basis of a pure Swave interaction, but for an educated practitioner of the few-nucleon problem it was clear that if such a huge extra binding was needed, the P-wave would not be enough.

The severe dynamical limitations of the first Glöckle paper were partly solved one year later together with Offermann [91], using S+P waves and the Reid potential including a tensor force. They found that when adding attractive P-waves in V nn , the critical S-wave scaling factor to bind 3n was slightly lowered to s = 3.7, but still remained unacceptably large. They also determined that when using a realistic nn interaction the most favorable candidate was J π = 3/2 -, dominated by a 3 PF 2 force but with a scaling factor of s ≈ 4. The other candidates were J π = 1/2 -, 3/2 + , 1/2 -, by order and in contradiction with Mitra [65], who had found the positive-parity state the most favorable one. Despite all these refinements they reached, however, the very same conclusion: no chance for a 3n bound state and a "low energy resonance is extremely unlikely". They ended with a non-prophetical remark: "The [required modifications in V nn ] are not so small that one likes to push forward a very expensive experiment to measure a three-neutron resonance". Would it be enough to discourage any further experimental or theoretical research? Certainly not! We are aware of other works that appeared during the following decade [92,93,94,95,96], but we are not going to comment on them with further detail for they represent still a preliminary stage, due to the interactions and/or the techniques used to solve the three-body problem.

At the end of this period, two works devoted to 3n and 4n resonances were published, which deserve special attention for they were also based on exact solutions of the three-and four-body problem.

In 1997 Sofianos et al. [97], within the framework of the hyper-spherical formalism and CSM, computed the position of the 3n and 4n resonances in the complexmomentum plane, identified with the zeros of the Jost function or the S-matrix poles. The nn interaction was

Re(E) [MeV]

Im(E)

[MeV] limited to two different central S-wave potentials, one purely attractive Yukawa and the semi-realistic MT13 [82]. With the MT13 potential, the lowest-lying 3n resonance was E = -4.9 -6.9 i MeV with J π = 1/2 -(degenerated with 3/2 -within the considered potentials), and it was located in the third quadrant of the second Riemann sheet (Fig. 9)1 . The search for a 4n resonance with the hyperspherical method gave for the nearest state J π = 0 + and E = -2.64 -8.95 i MeV, i.e. also located in the third quadrant of the second Riemann sheet. Although limited by using the "minimal approximation" in the hyper-spherical expansion and a semi-realistic NN potential, the conclusion from this work was clear. First, any 3n and 4n possible bound state is very far, since one requires an unphysically large enhancement factor in the NN interaction (typically s ≈ 4 and 3.6 respectively for 3n and 4n with the MT13 potential), well beyond any possible uncertainties of the models used. Second, the resonant states, although well identified, are located in the unphysical region with no possibility to manifest in an scattering experiment.

In 1999 Witala and Glöckle [98,99] published a solution of the Faddeev equations in configuration space for the 3n system. The nn interaction was taken from the phenomenological GPdT model [100] and from the realistic Reid93 Nijmegen potential [81], which includes the 1 S 0 , 3 P 0 and 3 PF 2 channels, all attractive. They introduced a different scaling factor s α in each partial wave of the nn potential, except in 1 S 0 to keep the dineutron unbound. They first forced the 3n system into a bound state and decreased s α progressively following the resonant trajectory in the complex-momentum plane. The J π = 1/2 -, 3/2 -, 3/2 + resonance positions were obtained by using the so-called Smooth Exterior Complex Scaling Method [101].

Due to numerical instabilities they could not reach the physical case with all s = 1. Their conclusion was however similar to the previous works: the 3n resonance energies will have such a large imaginary part that one should not expect to see any trace of them in experiments.

From GANIL to RIKEN

In 2002, few months after GANIL publication of the 4 n signal [START_REF] Marqués | [END_REF], there was another paper from Glöckle et al. concerning the 3n system [102]. They used rank-2 separable nn forces ( 1 S 0 , 3 P 0 , 3 P 2 , 1 D 2 ) and solved the Faddeev equations, which were analytically continued into the unphysical energy sheet below the positive real-energy axis. The trajectories of the 3n S-matrix poles in the J = 1/2 ± and 3/2 ± states were traced as a function of the scaling factors of the nn forces (Fig. 14), but only in the 3 P 0 , 3 P 2 and 1 D 2 channels in order to keep 2n unbound. The final positions of the S-matrix poles, once the artificial factors removed, were found to be in all cases far from the positive real-energy axis, which provided a strong indication for the nonexistence of nearby 3n resonances.

This was the last of a series of papers which were pioneers on this topic [90,91,98]. The overall conclusion was that 3n resonances could not exist near the physical region. The clarity of these works and the convincing arguments put forward should have been sufficient to discourage once and for all any claim for a 3n resonance. As we will see, this is unfortunately far from being the case.

Immediately after GANIL publication of the eventual observation of a tetraneutron, three consecutive papers were published [103,104,105] as a reaction to this surprising result. They agreed in the impossibility of a bound tetraneutron, based on the fact that the modifications of the NN (and eventually NNN) interaction required to obtain a 4 n bound state were so large that the description of all the surrounding nuclear chart ( 4 He, 4 H, 5 He...) would be destroyed, and not by a small quantity but by several tens of MeV. If the conclusion was unanimous, the methods used were however quite different.

In 2003 Timofeyuk [103] used Hyperespherical Harmonics (HH) formalism and a quite simplistic NN interaction (Reichstein and Tang [106]) to compute the lowest diagonalized hyper-radial potential V KK ′ (ρ) in the model space K max = 16, and all of them where found to be monotonously repulsive for the 2n, 4n, 6n and 8n systems. The possibility of a narrow resonance was also excluded in any of these systems. In a previous preprint [107] the same author had proposed alternative transfer reactions to investigate multineutron states and evaluated some cross section using the DWBA approximation.

Later in 2003 Bertulani et al.

[104] considered a totally different approach. They computed the interaction between two 2n dimers artificially bound in a simplistic NN model (Volkov [108]), hoping to find a possible attraction among them, as one could expect between two bosonic-like ( 1 S 0 ) systems. Whatever they tried, they systematically obtained a dimer-dimer repulsion. They concluded about the impossibility to explain any 4n bound or resonant state within this model, "although more complex variational approaches still can be explored", as well as any possible T = 2 4n forces. Despite the simplicity of the model and the crude approximations used in obtaining the dimer-dimer potential, their conclusion turns out to be quite general as it will be discussed later in the framework of EFT (Effective Field Theory) for fermions in the unitary limit. A few months earlier Pieper [105] had applied Green Function Monte-Carlo (GFMC) methods to study the possible existence of a 4 n bound state. He used the NN Argonne V18 interaction supplemented with the Illinois IL2 3NF [109], which had been very successful in reproducing the binding energies of A = 2-10 nuclei. The conclusion was the same than in previous works: any attempt to force a 4 n bound state would have "devastating effects" in the description of the nuclear chart. However, Pieper surprisingly claimed a possible 4 n resonance at about only 2 MeV. These two conclusions seemed somehow contradictory: the 4n Hamiltonian could hardly accommodate at the same time the absolute impossibility of a bound state and a near-threshold resonance.

The possible resonant state was determined by computing the energy of 4 n bound states in a Woods-Saxon trap (4) with fixed a [105]. Using different values of R, the potential strength V 0 was progressively decreased until the 4n became unbound, and the last energy values were extrapolated towards the absence of trap at V 0 = 0 (Fig. 15). All radii led to similar extrapolated values of E ∼ 2 MeV. Pieper's work had a great influence in this topic and it was taken as a definite proof to justify further experimental, and some theoretical, investigations. It is worth noting, however, that his own statement was very cautious: "since the GFMC calculation with no external well shows no indication of stabilizing at that energy, the resonance, if it exists at all, must be very broad".

In An interesting critical review of 4n results can be found in Sec. 4 of that paper. They explained why no indication of 4n FSI was found in the pion DCX reactions, and concluded that any resonant state, if existing at all, must be very broad. Furthermore, they computed with a reaction model how the position of the broad 4n peak can be shifted depending on the reaction used to observe it. They pointed out, with illustrative examples from standard Quantum Mechanics textbooks, that sharp resonance-like structures in the continuum can be formed even if only repulsive potentials are present in the system. In 2005 Lazauskas et al., motivated by the paradoxical results of Pieper, undertook a series of works to locate the nearest 3n [84] and 4n [85] resonances in the complexenergy plane using a realistic nn Reid93 interaction with CSM and ACCCM, adapted to properly deal with the continuum. The technique used for solving the three-and four-body problem was the Faddeev-Yakubovsky (FY) equations in configuration space, which were already tested in similar systems like 4 He, p+ 3 H, n+ 3 He, 2 H+ 2 H and n+ 3 H [111]. The only approximations made were the number of grid points and the partial-wave expansion in both interactions, and the FY amplitudes which were carefully checked. Several J π states were examined, being the nearest ones 3/2 -, 1/2 -, 1/2 + for 3n and 0 + for 4n.

The trajectories in the complex plane were plotted as a function of the strength of an artificial three-and fourbody force of the type (3), starting from a bound state (W 0 = W c 0 ) and moving to the continuum by decreasing the strength until the physical value at W 0 = 0. All of them were located very far from the physical region well before W 0 was removed, lying in the third quadrant of the second sheet, i.e. Im(E) < 0 and Re(E) < 0, often known as "sub-threshold resonances" (Fig. 9). We have displayed in the upper panel of Fig. 16 the J π =3/2 -, 1/2 -, 1/2 + trajectories of 3n obtained with the CSM and ACCCM. The lower panel shows the 0 -trajectory of 4n using AC-CCM, where the convergence of the method as a function of the degree (N, M ) in ( 5) can be observed.

The conclusion was that there could be no observable 3n or 4n resonance, a fortiori bound state, with realistic nuclear forces, in close agreement with the previous 3n work [102] despite the different interaction and methods used. It was pointed out that Pieper's conclusions were certainly due to the extrapolation procedure used.

For the sake of completeness, we note that in 2008 Anagnostatos [112] claimed that 4 n and 6 n could be particlestable. Within the framework of the so-called Isomorphic Shell Model, mixing semi-classical and classical mechanics, he concluded that, even if the results depended on the strength of V nn , at least low-lying resonant states could be reasonably expected, with 6 n being more favored than 4 n. These conclusions remain however very doubtful, for the very peculiar approach used to solve the few-body quantum mechanical problem in a classical way.

Post RIKEN 2016

After RIKEN publication of a new tetraneutron signal [4], the theoretical interest on this field was renewed and even enhanced. We remind that this new result was in fact also compatible with a bound state, while GANIL result was found to be also compatible with a low-lying resonance [50], and they were therefore in agreement.

In the same 2016, Hiyama et al. [86] published a joint effort between the authors of Refs. [84,85] and theorists from RIKEN, working independently and using different approaches for solving the few-nucleon problem. The work from 2005 [85] let no hope for the RIKEN signal to be attributable to a 4n resonance, and it was natural to clarify not only the tension with the experimental results but also with Pieper's conclusion.

The philosophy, as well as the computational methods, were different. From the physical point of view, the strategy was to keep unchanged realistic V nn forces and the well controlled V nnn forces in the T = 1/2 total isospin channel, and let as unique free parameter to reproduce RIKEN results the strength of V nnn in the T = 3/2 channel. The introduction of the T = 3/2 3NF was required to account for the systematic underestimation by the UIX model [109] of the binding energies in the He isotopic chain, although it was found to be smaller than the T = 1/2. The advantage of this approach was to preserve the good description by the modern nuclear Hamiltonians of all the systems not (or little) sensitive to this 3NF component.

From the technical point of view, the 4n solutions in configuration space were obtained in two different ways: by using the variational Gaussian Expansion Method [113,114,115,116,117,118] and by solving FY equation in configuration space [84,85]. The complex energy solutions where computed using both the CS and ACCC methods. The 3NF, for both T channels, were a sum of two hyperradial Gaussians (one attractive and one repulsive): For the T = 3/2 channel, the T = 1/2 parameters were used except for the attractive strength W 1 (3/2), the only parameter adjusted to reproduce RIKEN results. The critical value of W 1 (3/2) at which 4n becomes bound was determined for J π = 0 + , 2 + , 1 + , 0 -, 1 -... (by order of appearance). By decreasing this value, all of them moved into the resonance region, i.e. into the second Riemann sheet of the energy surface with E = E R -iE I . It was found however that the strength required to reproduce RIKEN results for the most favorable case (0 + ) was 36, -30] MeV, which is totally unphysical.

V 3N ijk = 3/2 T =1/2 2 n=1 W n (T )e -(
W 1 (3/2) ∈ [-
To have an idea, we can compare it with the corresponding W 1 (1/2) ≈ 2 MeV, while T = 3/2 should be smaller [109] and subdominant in an EFT expansion [121]. Taking into account that W 1 (1/2) is responsible for about 5 MeV of extra binding in 4 He, the additional energy required for the 4n system to become a resonant state is of the order of 100 MeV. Furthermore, it was found that starting from W 1 (3/2) ≈ -18 MeV the well-known unbound states 4 Li, 4 He * and 4 H become bound.

The trajectory of the 4 n(0 + ) resonant state in the complex-energy plane as a function of W 1 (3/2) is shown in Fig. 17. Already at half the strength required to reproduce RIKEN result, the real part of the resonance is E R ≈ 6 MeV and its width Γ ≈ 10 MeV. But even these values are unphysically large, and they completely miss the description of the well-known n+ 3 H cross-section [86]. The 3n was also studied in detail, and the conclusions were that the corresponding resonant states were located even further from the physical region than the 4n ones.

In 2017 the same authors completed the previous work [123,124] by computing in their numerical simulations the strength function used by RIKEN experimental group, following their suggestion. While the works of 2005 [84,85] had located the pole position of the nearest 3 n and 4 n resonances generated by realistic interactions, the question asked in 2017 was whether or not one could account for a 4n near-threshold resonance by modifying the most unknown part of the nuclear Hamiltonian: the 3NF for T = 3/2. Despite the different motivations and methods, all the conclusions were in line with the pre-GANIL era works: any near-threshold observable resonance in the 3 n and 4 n system (and a fortiori a bound state) must be definitely excluded. Concerning the RIKEN result, it was pointed out that there exist quantum mechanical mechanisms, other than S-matrix poles, that could generate a low-energy enhancement in a scattering cross-section.

All in all, it seemed clear at this stage that Pieper's prediction of a few MeV tetraneutron resonance should be considered as an isolated case in the theoretical landscape, maybe a consequence of the trap artifact or of his cavalier extrapolation to the continuum, but in any event devoid of any physical ground. However, if this was clear for some people, it was not clear for everybody.

Indeed, in the same 2016 Shirokov et al.

[125] found a 4 n resonance with parameters E R = 0.84 MeV and Γ = 1.38 MeV, in perfect agreement with RIKEN results. They used their non-local NN interaction JISP16 [126] and attempted several methods to describe the 4 n state predicted by Pieper [105]. The JISP16 V N N is expanded in terms of HO radial basis R nl as:

V LL ′ SJ (r, r ′ ) = n(L) n=0 n ′ (L) n ′ =0 R nL (r) C LL ′ SJ nn ′ R nL ′ (r ′ ) with n(L) = n ′ (L) = 4.
It is thus a rank-4 separable interaction, whose C nn ′ coefficients and b parameter of the HO were adjusted to successfully reproduce the lowenergy properties of NN and some light nuclei. It can be considered as a high-precision potential.

First they tried to solve the 4n problem with the No-Core Shell Model (NCSM) by scaling the interaction and tracking the lowest state as a function of the scaling factor s. They performed extrapolations to the unbound regime similarly to Pieper [105] and found results in quantitative agreement, i.e. a resonance at E R ∼ 2 MeV without any information on the width. Then they used No-Core Gamow Shell Model (NCGSM) techniques, which directly provide the resonant parameters in the complex energy plane. In this way they obtained a relatively broad resonance at E R ∼ 2.5-3 MeV and Γ ∼ 2.5-6 MeV. They noted however that the resonance parameters decreased when the size of the NCGSM basis increased.

Finally they used extended NCSM techniques to deal with the 4n in the continuum. In particular, in order to calculate the 4 → 4 S-matrix resonant parameters, they used a Single-State Harmonic Oscillator Representation of Scattering Equations (SS-HORSE), an extension of the NCSM Hamiltonian in an HO basis by an infinite kinetic energy matrix. For this extension they used the so-called "democratic decay approximation", which describes the 4n continuum in terms of an hyper-spherical harmonics (HH) basis, limited in practice to its "minimal approximation", i.e. with hyper-spherical momentum K = 2.

The SS-HORSE approach provides the 4 → 4 scattering matrix S(E), and the corresponding phase shifts δ(E), at the positive eigen-energies of the NCSM Hamiltonian, from which they extrapolated the position of the S-matrix poles following some analytic expressions. These results depend critically in the parametrization of S(E) in terms of four adjusted parameters. The values so obtained for the resonance parameters of the S-matrix pole were considered "surprisingly small". Their final result, almost exactly the same than RIKEN one, was obtained by further adding a false pole plus an additional fitting parameter, and was presented as a "prediction".

The findings of Shirokov, in line with Pieper's, were in strong contradiction with all the previous ones obtained with direct FY and Gaussian methods, and brought even more confusion to an already quite confuse situation.

In 2017 Gandolfi et al.

[127] used Pieper's methods (Woods-Saxon trap, QMC techniques and extrapolation to the continuum) and reached a similar result: a 4n nearthreshold resonance at E R = 2.1(2) MeV. Furthermore, they concluded that a 3n resonance was "lower in energy than the tetraneutron one", at E R = 1.1(2) MeV. This was in fact the main question raised in the paper title, and with their positive answer the authors encouraged further experimental measurements. The resonance widths could not be estimated, but this work improved in several aspects the previous QMC calculations: the accuracy of the computed binding energy, the development of a consistent NN and 3NF chiral EFT interaction at nextto-next-to-leading order (N2LO), and the introduction of a quadratic instead of a linear extrapolation to reach the continuum. However, as we already noted and will develop below, this approach suffers from two inherent redhibitory drawbacks: the use of a trap itself and the procedure of extrapolation to the continuum. One interesting conclusion was that, due to the extreme diluteness of the system, the roles played by the 3NF and the details of the NN one were very small. As a check of their trap results, they repeated the calculation by introducing an overall scaling factor s on the nn interaction. Surprisingly, they found that a value of only s = 1.3 was enough to bind the 4n, in sharp contradiction with all the previous works where an s value 2-3 times bigger had been found. Pieper had also concluded that changes more dramatic than 30% were needed. Even so, the extrapolated result at the physical value s = 1 was E = 2.0(1.0) MeV, consistent with their trap results. Also in 2017, Fossez et al.

[128] used two ab initio techniques and several NN interactions (chiral EFT and JISP16) to study the existence of a 4n resonance, which "would deeply impact our understanding of nuclear matter". The first technique was the NCGSM (as in Ref. [125]), which is a NCSM extension into the complex-energy plane with couplings to the continuum by means of the Berggren basis, that incorporates some selected resonances. The second technique was the Density Matrix Renormalization Group (DMRG), an alternative way to solve the nuclear many-body problem in the continuum. Both methods were supplemented by the use of natural orbitals and a new identification technique for broad resonances.

The nn interaction was scaled by a global factor s in order to bind the 4n system, which was then decreased until it reached the physical state s = 1. The results depend strongly on the particular technique (NCGSM or DMRG) and the different approximations made (see Fig. 18). For example, at s = 2 the 4n system is bound in all the approximations, but its binding energy ranges from ≈ 2-6 MeV. Only in the less accurate calculation (NCGSM-2p2h) the physical value could be reached, with E R ≈ 7 MeV and Γ ≈ 4 MeV. The other calculations could not go below s ≈ 1.5. They observed a clear trend in growing width as the 4n continuum was more accurately described. They even noted that "the opening of new decay channels and the presence of continuum states in the configuration mixing above the threshold is expected to make the width explode when s → 1, especially in the DMRG results where all decay channels are open".

In view of this dispersion it is very difficult for the reader to draw a conclusion. In any event, the authors confirmed the existence of an S-matrix pole of the 4n system associated with J π = 0 + . They concluded that while the energy position of a 4n resonance might be compatible with the RIKEN value when s → 1, including more than two particles in the continuum suggested that the width would be larger than Γ ≈ 4 MeV.

Although sharing some methodological aspects, this work is much more careful in the analysis than Shirokov's, and is free from daring extrapolations. The only objections we can think of, besides the limited model space and the dependence on the HO parameter mentioned by the authors, is the use of a global scaling factor s and the anomalous transition from an artificially bound state into the continuum that it generates. The appearance of 2 n and 3 n open channels can substantially modify the properties of the computed resonance when the system goes into the continuum. One may find intriguing that in all the considered schemes the critical scaling factor determining the 4n threshold lies in the range s ∼ 1.5-1.9, sensibly smaller than in most previous works (except Gandolfi's s = 1.3 [127]). It is also unclear whether the Berggren basis used to extend NCSM into the continuum is adapted to describe a broad resonance, or if the results are conditioned by the particularities of the Berggren basis itself.

In 2018 Deltuva studied 3n [129] and 4n [130] states in the framework of AGS equations in momentum space [131]. This formalism, equivalent to FY equations, provides rigorous solutions for the three-and four-body problem and describes very accurately all the corresponding nuclear systems around the breakup thresholds [132,133,134]. An important advantage of the AGS approach with respect to previous calculations is its ability to compute the resonance position and its effect on the scattering amplitudes that could lead to observables in a scattering process, in particular the 3n → 3n and 4n → 4n transition amplitudes.

The trineutron [129] was studied using several realistic nn interactions and with a scaling factor s α only in the non-zero angular momentum partial waves, in order to avoid a bound 2 n in 1 S 0 that would open a threshold in the 3n system and complicate the analysis. For sufficiently large values of s α , a 3 n bound state was found for several J π and at critical values within s ≈ 6-8, depending on the nn interaction model. When decreasing s α , a resonance behavior was observed in all the transition amplitudes, as shown in Fig. 19 for the 3 PF 2 case. Clearly, the resonant behavior for s = 6.2 weakens when decreasing s and completely disappears for the physical interaction strength. Note that even with s = 4 the amplitude is totally flat. The corresponding pole trajectories were also extracted, and for the Reid93 interaction they matched the results from FY calculations [84]. Deltuva concluded that there are no physically observable 3n states consistent with the presently accepted interaction models, and that the pole trajectories were similar for all the realistic nn potentials and in close agreement with Refs. [84,86].

The tetraneutron [130] was studied using the same rigorous methodology. The 4n transition operators T βα in the 0 + state were calculated using a physical nn potential, i.e. s α = 1. As in the 3n case it showed no indication of resonances, but for sufficiently large s α a resonant behavior was clearly seen in all matrix elements. The tetraneutron became bound at s = 5.29, thus lower in energy than the trineutron, bound only at s > 6. The results were found to be almost independent of the particular nn potential, as in previous works. Interestingly, it was found that even in the absence of any 4n resonance, some 4 → 4 transition operators exhibit a low-energy enhancement (Fig. 20, bottom insert). It was conjectured that they could also manifest in other reactions with the 4n subsystem in the final state, like the 4 He( 8 He, 8 Be) one measured at RIKEN.

The conclusions of these two papers were in close agreement with the exact FY methods [84,85,86,97,98], but in sharp contradiction with Refs. [105,125,128], and even with Ref.

[127] that had used the same EFT potential.

In 2019 Deltuva and Lazauskas shed some light on this striking difference [135,136]. The first paper was a critical Comment [135] to the conclusions of Ref. [127]. In particular, they analyzed quantitatively some serious shortcomings related to confining neutrons in a Woods-Saxon trap (4) and to the misleading extrapolation procedure without taking into account the analytical structure of the thresholds. Two problems were pointed out: In the second paper [136] Deltuva and Lazauskas studied the consequences of a bound2 n, obtained by enhancing V nn with a scaling factor, in the evolution towards the continuum of an artificially bound 4 n state. The results were obtained by combining the FY [84,85] and AGS [130] formalisms. Two different cases were considered:

1. The same scaling factor s was used in all partial waves.

They found that a 0 + bound 4 n indeed emerges for a deep enough 2 n, with s ≈ 2.4-2.7 depending on the particular potential. When slightly decreasing s, this artificially bound 4 n, although having negative energy, is not longer bound but decays into two 2 n, in contradiction with Refs. [127,128] that found a bound 4 n for s ≈ 1.5. Furthermore, by continuing to decrease s, they found that the 4n system will never evolve into a resonance but into a virtual state, similar to the NN 1 S 0 case. The binding is then dominated by the S-wave and it is clear that there is no way to build a 4n resonance with a global scaling factor. 2. The scaling factor for the 1 S 0 interaction was fixed to s = 1.3 (with B nn = 0.316 MeV) and the main extra binding was provided by scaling P-and higher angular momentum partial waves. This breaking of S-wave dominance is essential to produce P-wave resonances in the final state. Such a resonance is indeed formed with s P ≈ 3.9, but when slowly removing the scaling factor the resonant character disappeared well before reaching the physical case s = 1.

Their main result, the fact that one scaling factor cannot generate a 4n resonance, is enlightening and very relevant in view of the many studies that use this artifact.

Later in 2019 Li et al.

[138] contradicted the two previous works. They used ab initio NCGSM and concluded that 3 n and 4 n are both observable resonant states. The 4 n resonance parameters were (E, Γ ) 4n = (2.64, 2.38) MeV. They even claimed that the 3 n was lower in energy and narrower, (E, Γ ) 3n = (1.29, 0.91) MeV, and strongly encouraged experimentalists to search for the trineutron at low energy. These calculations, that pretend to be "nearly exact", are however based in the same trap methods of Refs. [105,127], and we have already discussed the kind of redhibitory problems from which these methods suffer.

The multineutron field is still very active, as shown by the two additional works devoted to this topic that have been already published in 2020, and with which we will close this theoretical review.

In the first one, Higgins et al.

[139] studied the 3n and 4n systems within the adiabatic hyper-spherical framework. The corresponding adiabatic potential-energy curve was analyzed and found to be repulsive, in agreement with Timofeyuk [103]. They concluded that there is no sign of a low-energy resonance for any of these systems. However, they observed in both of them some low-energy enhancement of the Wigner-Smith "time delay", which for a single adiabatic potential is defined as Q(E) ≡ d [δ(E)]/dE, that could provide a hint to understand the GANIL and RIKEN near-threshold enhancements.

Later in 2020 Ishikawa [140], using FY and realistic NN potential (AV18), found no evidence of a 3n resonance. Moreover, he examined different methods used to extrapolate the 3n energy from an artificially bound state to the continuum, and found that the enhancement of nn P-wave alone or the use of 3NF do not generate spurious 3n resonances. However, using external trapping potentials leads to positive 3n energy results, which without a further careful analysis may be considered as a resonance while they are in fact a general defect of the trapping method. In this way, Ishikawa explained the contradictory results of Refs. [127,138] concerning 3n states.

Conclusions from theory

At first glance, the striking disagreement among different theoretical works presented in this section is obvious, shocking, and could be even devastating for the credibility of the field. From a general point of view one can distinguish two main families of results, whose conclusions are orthogonal to each other.

In the first group, the existence of near-threshold 3n and 4n resonant states is claimed. In the case of the tetraneutron, those works seem thus in agreement with the two experimental signals reported at GANIL and RIKEN. However, for the second group this is a scientific nonsense since an explicit calculation of these states finds them very far from the physical region, in the third quadrant of the unphysical energy sheet. According to them, the recent experimental results, if confirmed, in no way can be attributed to a 3n or 4n resonant state and must have another origin, still to be clarified.

If these two points of view can hardly be more opposite, there is at least a general consensus in that the origin of such contradictory conclusions must not be found in the different interactions used. Indeed, for these low-energy physical states with a large spatial extension, any interaction providing acceptable low-energy parameters leads to similar qualitative results. Furthermore, the role of threeneutron forces is negligible. The differences among them must rather be found in the methods used to solve the few-nucleon problem and/or in the way they access the few-neutron continuum.

All results obtained using FY or the AGS equations agree with each other about the non-existence of any nearthreshold3 n and 4 n. These approaches represent days the most rigorous framework for solving the fewbody problem, and include in a natural way all the continuum states, without additional approximations or ansatzs. Their agreement is not only qualitative but also quantitative, in particular when computing the critical values of the scaling factors or the trajectories of the resonance position in the complex-energy surface. Moreover, this agreement is extended to the variational Gaussian Expansion methods, when the resonant states are tackled with CSM or ACCCM, and to HH, where they unanimously find a strong repulsion in the adiabatic potential curves disallowing any bound or resonant state.

Furthermore, we would like to mention the EFT developments in multifermion systems close to the unitary limit, as it is the case of multineutron states, where they found a universal dimer-dimer and fermion-dimer strong repulsion [141,142,143], in contrast with any resonant state and in line with the above mentioned series of works.

This impressive bulk of coherent results is in frank contradiction with the two GFMC results, as well as with the NCSM-SS-HORSE and the two NCGSM works, the latter three being based on many-body techniques.

GFMC [105,127] is accurate for computing bound states and has been able to deal with a larger number of particles than FY and AGS, but it is clearly less adapted to scattering problems, in particular the computation of resonances. Moreover, we have shown that their conclusions are strongly weakened due to two important additional issues: the use of traps and the method employed to extrapolate the (supposed) bound state into the continuum.

NCSM-SS-HORSE [125] is based on several many-body techniques that were not checked in similar few-nucleon scattering problems. Moreover, the 4n continuum is treated a very indirect way, with the "democratic decay" hypothesis for extracting the 4n → 4n amplitude (neglecting the FSI between virtual dineutrons [136]) and a non-trivial extrapolation procedure to extract the complex parameters of an eventual resonance. Finally, they needed the further addition of a "false pole" in order to be in perfect agreement with the previously published RIKEN result.

The more comprehensive NCGSM-DMRG work [128] exhibits a manifest dependence in the approximations used to solve the four-body problem, and their global scaling factor could reach the physical case only for the less accurate of them. However, they observe a clear trend of an increasing width of the state, which according to the authors "exploded when s → 1", and concluded with a resonant state that would be much larger than the reported RIKEN limit. We believe that the use of a global scaling factor is a major problem, which due to the dineutron formation could never generate a resonance but a virtual state [136]. The most recent NCGSM work [138] suffers from most of the previous drawbacks plus the use of a trap.

Finally, any method that does not take properly into account the thresholds can lead to uncontrolled conclusions, and we believe that this is a common problem of the references in this second group.

With respect to the experimental results, we cannot find any theoretical explanation for GANIL and RIKEN events, both corresponding to an enhancement of lowenergy 4n in a given reaction cross-section, either breakup or DCX. We have presented in this section convincing arguments against the attribution of this enhancement to a resonance, at least in the standard way a resonance is understood in nuclear, hadronic and particle physics.

However, as it has been emphasized many times and by many authors [85,86,110,124,130,139], there are other quantum mechanical mechanisms that generate resonancelike structures, even for repulsive interactions. Deltuva's results of Fig. 20 show it explicitly in the 4n → 4n nonresonant transition amplitude. A much simpler case is illustrated in Fig. 21, corresponding to the P-wave scattering of two identical particles ( 2 /m = 1) on a repulsive potential well V (r) = +V 0 × Θ(R -r), where Θ is the Heaviside step function. By changing the parameter R, one can find a rich variety of resonance-like structures, with the corresponding E-dependent phase shift crossing π/2 close to the resonance peak 3 . Ignoring the underlying dynamical content, these structures can be attributed to a resonant state, very close to binding. We believe that this kind of behavior is also manifest in the time-delay analysis of Ref. [139].

Whether or not the experimental search of this quantum mechanical enhancements is pertinent enough to jus- tify great investments can be submitted to debate. However, associating these enhancements to physical resonant states would be misleading. 4 Sit finis libris... "Sit finis libris, non finis quaerendi" (this may be the end of the book, but not the of the quest), was an elegant way to close an open question in the scholastic times. Despite all the controversies and debates, this quest is definitely a fascinating story. Think of the distant 1960s, with experimentalists trying to do 'magic' with pions in order to create neutron matter in the laboratory; others going to nuclear reactors, looking for their grail in such a radioactive environment; others manipulating ion beams, hoping to add an remove nucleons against all odds. In parallel, theoreticians starting the first three-body calculations, despite lacking at that time the sophisticated tools and inputs that we have seen were needed. Sixty years later the quest is still fascinating... but still open. If from the theoretical point of view the question of a narrow 3n and 4n resonance seems closed, there is an urgent need for a definitive experimental conclusion concerning the GANIL and RIKEN events.

Nuclear theory, and in particular the description of light nuclei from the underlying forces between nucleons, has seen an enormous progress working hand in hand with experiments. The main theoretical models have been born from the known properties of stable nuclei, and then they have been refined and updated with the increasing knowledge about more and more exotic isotopes. In the field of neutral nuclei, however, the theory advances almost blind. The only two experimental signals [START_REF] Marqués | [END_REF]4] are still weak and do not provide any firm and precise observable that could benchmark the different methods and techniques.

In this respect, it should be a priority to confront the controversies surrounding the tri-and tetraneutron calculations with the help of high-statistics, unambiguous experimental results. New techniques or experiments that would provide still another weak and/or ambiguous signal do not seem the best way to unlock progress in the field. It is our hope that the new experiments already (or soon to be) undertaken at RIKEN will bring the very much needed reliable reference points for the theoretical conclusions to be adapted. If the positive signals were refuted, then the calculations predicting resonances should be reevaluated. If they were confirmed though, then the detailed characteristics of the signals (such as energies, widths, or angular correlations) should be explained.

Once the 3n and 4n systems will be clarified, both from theory and experiment, it will be easier to move on the heavier systems, in particular 6n and 8n. Even if the lightest multineutron states were not observable, a pertinent question remains: when increasing the mass number, can at some point a multineutron manifest as a bound or resonant state? Even for the most reluctant theories, and of course for experiments, this is still an open question. The analogy between neutrons and atomic 3 He is very suggestive, for their fermionic character and for the form of their interaction, with a hard repulsive core and a small attractive well. The dineutron (with a critical enhancement factor s ≈ 1.1) is even closer to binding than the 3 He dimer (s ≈ 1.3), and we know that a 3 He drop becomes bound beyond N ∼ 30 atoms [144]. On the other hand, the nn P-waves are smaller and its centrifugal barrier is less peripheral than in the helium atomic case [85,145], and that could make a difference.

The known as "helium anomaly", i.e. the fact that 8 He is more bound than 6 He, has long been a clue towards the search for light multineutrons, taken as a hint for additional binding due to the increasing number of neutrons on top of 4 He. Less known but maybe more spectacular would be the analog "hydrogen anomaly", provided it is established. The With the recent progress in the exact calculation of the five-body problem, 7 H is now accessible almost ab initio (treating the triton as a particle) for the theory and, as we have seen at the end of Sec. 2, a related experiment that should provide very high statistics and resolution has been already undertaken. This super-heavy isotope of hydrogen could in this way be the key for the next steps in the field, both for the role of the tetraneutron at the t + 4n threshold but even for the hexaneutron at the p + 6n threshold. Going beyond A = 6 represents still a too important obstacle for experiments, and for the moment should be left for the theory, once the lighter multineutrons have helped to provide a solid base for the models.

In any event this domain will remain fascinating since, paraphrasing the authors of Ref.

[128], it "would deeply impact our understanding of nuclear matter".

A Two neutrons in a trap

In order to illustrate the kind of problems related to the extrapolation from bound states to the continuum, we present in this Appendix a detailed solution of the simplest case, in the same spirit as it was considered in Ref.

[135]: two neutrons confined in a trap.

Let us consider two neutrons interacting via V nn and furthermore submitted to a potential trap having, as in Refs. [105,127,138], a Woods-Saxon form:

V T (r i ) = V 0 1 + e r i -R a (7)
The total two-body Hamiltonian is:

H = - 2 m n ∆ + V nn (r 2 -r 1 ) + V T (r 1 ) + V T (r 2 ) (8)
The particle coordinates r i can be expressed in terms of the relative (r = r 2r 1 ) and center-of-mass (2R = r 1 + r 2 ) coordinates as:

r 1 -r 2 = 2r 1 -(r 1 + r 2 ) ⇒ r 1 = - r 2 + R r 2 -r 1 = 2r 2 -(r 1 + r 2 ) ⇒ r 2 = + r 2 + R
Assuming R ≡ 0 to get of the center-of-mass motion, the intrinsic two-body Hamiltonian takes the form:

H = - 2 m n ∆ + V nn (r) + 2V T r 2 (9)
The numerical solution of this problem is trivial and can be obtained with any required accuracy. For illustrative purposes, we have considered a phenomenological CD MT13 nn interaction: At first glance it seems that both series of points can be indeed linearly extrapolated (solid lines) towards a positive value (empty square at V 0 = 0), which could be interpreted, as in Refs. [105,127,138], as a 2n resonance with positive energy. Furthermore, it is worth noting that even in the absence of any interaction, the results of the linear extrapolation (dashed lines) display a similar trend and would lead to a very disturbing resonant state between two non-interacting neutrons! One could argue that in the case of two neutrons in a 1 S 0 state, the extrapolated results of Fig. 22 correspond in fact to the nn virtual state with positive energy, as it is often considered. We remind, however, that the virtual state, corresponding to a pole in the negative imaginary axis k 0 = +iκ, has negative energy k 2 0 = -κ 2 . The pole trajectory in the complex k-plane and E-surface of an artificially bound nn (k i and E i ∈ E I ) into the continuum (k f and E f ∈ E II ) is summarized in Fig. 23. As one can see, both bound and virtual (unbound) states correspond to negative-energy values in the E-surface, though living in different Riemann sheets.

V nn (r) = V r exp(-µ r r) r -V a exp(-µ a r) r (
If the extrapolation analysis sketched in Fig. 22 is done carefully, in particular using much smaller binding energies, one can see that a strong curvature in the E(V 0 ) dependence starts at B ≈ 0.1 MeV until the continuum threshold E = 0 is reached at a critical value V c 0 (R), which increases with R. For instance, V c 0 ≈ 0.1 MeV for R = 5 fm (see the zoom displayed in Fig. 24).

After the threshold, the V 0 -dependence represented in Fig. 24 changes drastically the curvature, and turns downwards toward negative values of E as can be see in Fig. 25, and as expected from the location of a virtual state in the continuum illustrated in Fig. 23. The highly non-trivial trajectory in the (E, V 0 ) plane displayed in Fig. 25, results from the ACCC extrapolation given by equation (5). To be faithful, a set of bound state values (V n 0 , E n ) must be computed with very high accuracy, as accurate should be the critical values of the potential strength at threshold V c 0 , when the system enters the continuum. As an example of the accuracy required, the values used for the case R = 5 fm are listed in Table 1.

The trajectory in the (E, V 0 ) plane is represented in the complex k-plane in the left panel of Fig. 23. It is the straight line formed by the union of γ ∪ {0} ∪ γ ′ . The corresponding path in the E-surface is Γ ∪{0} ∪Γ ′ , shown in the right panel, where E i belongs to E I and E f to E II . The transition from bound state to continuum is at k = 0 and coincides with the transition E I → E II . Note that both initial (k i ) and final (k f ) states have negative energy (E i and E f respectively), despite the fact that E f is in the continuum, but belong to different Riemann sheets.

From the above results it is clear that a simple polynomial extrapolation cannot describe the non-trivial evolution of a bound state into the continuum displayed in Fig. 25. This is in fact the main conclusion of this Appendix. ACCCM can do it accurately by taking into account the analytical properties related to the thresholds. However, the method requires the previous computation of several bound state values with high accuracy and very close to E = 0. This can be easily done in the two-body case we have considered above, although it becomes more and more involved when increasing the number of particles. This level of accuracy is not accessible to all the methods used to solve the few-nucleon problem, but it Fig. 25. Energy of 2n in a trap, from bound state to continuum. The corresponding trajectory in the complex-momentum plane (γ ∪γ ′ ) and the E-surface (Γ ∪Γ ′ ) is displayed in Fig. 23. V0 E 0.20 -0.036938 0.15 -0.009817 0.14 -0.006220 0.13 -0.003367 0.12 -0.001331 0.115 -0.00064676 0.10384(1) 0 Table 1. Energy values of 2n in a trap used in the ACCC extrapolation (5), as a function of the potential strength V0 for R = 5 fm. The last point corresponds to the threshold extrapolated value.

has been successfully applied to A = 3, 4, 5 in Ref. [84,85,86,87].

The case presented above corresponds to an S-wave virtual state, which is different from the realistic 3n and 4n P-wave resonances displayed in Figs. 14,16. They have however in common that the trajectories of the corresponding pole positions end at negative energy, Re(E) < 0, in the third quadrant of the second Riemann sheet E II , as illustrated in Fig. 23. Of course, it is clear that a polynomial extrapolation of the results displayed in Fig. 22 (or Fig. 15) is forced to end at positive energy values, and is therefore unable to reproduce the realistic behavior of the 3n and 4n resonant states.

Furthermore, it is worth noting that even in the case of positive energy P-wave resonances, with Re(E) > 0 and Im(E) < 0 like the one represented in Fig. 17, a polynomial extrapolation disregarding the analytic threshold structure can lead to wrong results. This can be illustrated by considering the same nn interaction (10) but suitably scaled in order to generate a P-wave nn resonant state. For instance, with a scaling factor s = 3.0 one has a resonance at E = 4.40 -3.32 i MeV, that is a state with Γ ≈ 6.6 MeV, still smaller than the realistic case of Fig. 17. To complete our results for nn S-wave, we have computed the energies of this P-wave resonant state in the trap, as a function of V 0 for two different parameters of the trap size R, and performed a quadratic extrapolation towards the 'free' case V 0 = 0. Results are displayed in Fig. 26. They show a significant disagreement between the extrapolated values and the physical resonance energy, indicated by a blue solid circle. We cannot exclude, as noted in Ref.

[135], that for very narrow P-wave resonances the difference could be small, but the polynomial extrapolation is not reliable in a general case. This could be the case in the unphysical example (P-wave resonance in an S-wave two-Gaussian potential) considered in Ref.

[127] in order to justify their extrapolation method. The exact solution of two neutrons in a trap can be straightforwardly extended to the 3 n case. The Jacobi coordinates of the 3n system are expressed in terms of the particle coordinates r i as:

y 1 = 2 √ 3 r 1 - r 2 + r 3 2 = 2 √ 3 r 1 + r 1 2 - r 1 + r 2 + r 3 2
Since in the center of mass r 1 + r 2 + r 3 = 0, one has:

r 1 = y 1 √ 3
and equivalent expressions for the other r i . The three-body intrinsic Hamiltonian takes the form:

H = H 0 + V nn (x 1 ) + V nn (x 2 ) + V nn (x 3 ) + W (y 1 , y 2 , y 3 ) with:

W (y 1 , y 2 , y 3 ) = V T y 1 √ 3 + V T y 2 √ 3 + V T y 3 √ 3
Since the different sets of Jacobi coordinates are related by the linear transformations: y 2 = c 21 x 1 + s 21 y 1 y 3 = c 31 x 1 + s 31 y 1 the problem is equivalent to adding a three-body force: W (x 1 , y 1 ) = W [y 1 , y 2 (x 1 , y 1 ), y 3 (x 1 , y 1 )] which in the case of three identical particles can be easily solved with a single Faddeev equation. This calculation requires however some care due to the presence of open 3 n → 2 n + n channels for some values of the trap parameters (V 0 , R), but will not be developed here.

  [nb/sr/(MeV/c)]

Fig. 4 .

 4 Fig.4. Schematic view of Détraz's experiment[36]. In the collision of protons and tungsten a neutral cluster was supposed to be produced, go through a screen, and induce ( A n, xn) reactions on natural zinc. The new72 Zn isotope was detected through its decay into72 Ga * with a germanium detector.

Fig. 5 .

 5 Fig.5. Energy spectrum of 10 C from the 7 Li( 7 Li, 10 C)4n reaction at θ = 7.4 • . Known contaminant reactions are indicated either explicitly or with an asterisk. The red curve corresponds to five-body phase space, and the range in yellow to the region expected for a bound tetraneutron. Adapted from Ref.[44].

Fig. 7 .

 7 Fig.7. Missing-mass spectrum of the 4 He( 8 He, 8 Be)4n reaction. The solid (red) curve represents the sum of the direct decay of correlated 2n pairs plus the estimated background. The dashed (blue) curve represents only the latter, multiplied by a factor of 10 in order to make it visible. The 4 events at threshold are highlighted in yellow. Adapted from Ref.[4].

Fig. 8

 8 Fig.8. The solid histogram represents the number of experiments reviewed that were searching for multineutrons (36 total, mainly for tri-and tetraneutrons) as a function of the year of publication of the results. In yellow are those searching specifically for the tetraneutron. The stars represent the three positive signals reported, the empty one that was refuted[36] and the two solid ones that have not been contested yet[START_REF] Marqués | [END_REF]4]. The pale Gaussians guide the eye through the recurring pattern.

Fig. 9 .

 9 Fig. 9. Complex momentum plane and energy manifold.

Fig. 10 .

 10 Fig. 10. The principle of the Complex Scaling Method.

Fig. 11 .

 11 Fig. 11. Scaling factors (2) needed to bind 2n in the 1 S0 state for some nn interactions. AV18 and Reid93 give s 1.08.

Fig. 12 .

 12 Fig. 12. Possible decay channels of 3n and 4n in a trap.

Fig. 13 .

 13 Fig.13. Schematic view of the ACCC method. An ensemble of bound state positions (λn, kn) is computed for λ > λc, a critical branching point where k = 0, and the trajectory of the S-matrix resonance pole for λ < λc is determined.

Fig. 14 .

 14 Fig. 14. S-matrix pole trajectory of 3n in J = 3/2 -with a rank-2 separable potential. Adapted from Ref. [102].

Fig. 15 .

 15 Fig. 15. Energies of 4n in a Woods-Saxon trap (4) as a function of the potential strength V0. The results obtained with several values of R are used to extrapolate linearly into the continuum. Adapted from Ref. [105].

  2004, Grigorenko et al. published a work devoted to 4n and 5 H in the framework of the HH method [110].

Fig. 16 .

 16 Fig.16. Selected resonance trajectories for the 3n (upper) and 4n (lower) states as a function of the strength W0 of the artifact. Results were obtained by combining ACCC and CS (solid points) methods. Adapted from Refs.[84,85].

Fig. 17

 17 Fig. 17. 4 n resonance trajectory for the J π = 0 + state. Circles correspond to AV8 ′ potential and triangles to INOY04 [122]. The 3NF strength W1(3/2) was changed from -37 to -16 MeV in steps of 1 MeV for AV8 ′ , and from -36 to -24 MeV in steps of 2 MeV for INOY04. The RIKEN resonance region is indicated by the arrow at the top.

  Fig. 17. 4 n resonance trajectory for the J π = 0 + state. Circles correspond to AV8 ′ potential and triangles to INOY04 [122]. The 3NF strength W1(3/2) was changed from -37 to -16 MeV in steps of 1 MeV for AV8 ′ , and from -36 to -24 MeV in steps of 2 MeV for INOY04. The RIKEN resonance region is indicated by the arrow at the top.

r 2 ij +r 2 jk +r 2 ki )/b 2 n

 22 P ijk (T )(6) where P ijk (T ) is a projection operator for the total threenucleon T state. The parameters (W n , b n ) for the T = 1/2 channel were already fixed to reproduce the binding energies of some selected light nuclei and 4 He(e, e ′ ) 4 He * electromagnetic transition form factors [119], in conjunction with the AV8 ′ version of the Argonne NN potential [120].

Fig. 18 .

 18 Fig.18. Evolution of the energy and width (shaded area) of the 4n system with the scaling of the N3LO interaction from 2 to 1. The circles represent the NCGSM results with two neutrons in the continuum, which is used to generate the NCGSM results based on natural orbitals with two (triangles) and three (squares) neutrons in the continuum. The DMRG results without truncations are represented by stars. The RIKEN experimental energy is indicated by the solid circle, and the yellow area shows the uncertainty. Adapted from Ref.[128].

Fig. 19 .Fig. 20 .

 1920 Fig. 19. 3n transition amplitude for different values of the scaling factor in 3 PF2 partial wave, taken from Ref. [129].

Fig. 21 .

 21 Fig. 21. Resonant P-wave cross-section in a repulsive potential of the form V (r) = +V0 when r < R, for different values of R.

10 )

 10 reproducing the nn scattering length a nn = -18.59 fm and effective range r nn = 2.94 fm values, with the parameters V r = 1438.72 MeV, µ r = 3.11 fm, V a = 509.40 MeV, µ a = 1.55 fm, and 2 /m n = 41.44250 Mev fm 2 [83].The binding energies for two values of the trap size (R = 3 and 5 fm) and fixed diffusiveness a = 0.65 fm are plotted in Fig.22as a function of the strength 0 .

Fig. 22 .Fig. 23 .

 2223 Fig.22. Energy of 2n in a trap(7) as a function of the strength parameter V0 for two values of the trap size R. A linear extrapolation (solid lines) converges towards a positive energy value. The free case Vnn = 0 (dashed lines) displays a similar trend.

Fig. 24 .

 24 Fig.24. Zoom of Fig.22showing the behavior of the energy when approaching the continuum threshold at V c 0 and E = 0.

Fig. 26 .

 26 Fig.26. Energy of 2n P-wave resonant state in a trap(7) as a function of the strength parameter V0 for two values of the trap size R. Quadratic extrapolations (dashed lines) of the computed binding energies in the trap converge at V0 = 0 towards positive energy values, sensibly different from the exact result (blue circle).

  1. For some trap parameters [127], e.g. R = 6 fm and a = 0.65 fm, 2 n starts being bound at V 0 ∼ 0.09 MeV. Thus, when 4 n is supposed to become bound in the trap, 2 n is already bound by several MeV (see Appendix for details) and 4 n should in fact decay into 2 n(s) and maybe even 3 n, provided the theoretical dynamical method was able to account for open scattering channels. In any case, the 4n states of Refs. [105, 127] fulfilling E 2n < E 4n were not really bound states,

	but some discretized continuum states that do not
	evolve into a resonance. The extrapolation of their en-ergies into the continuum is misleading 2 .
	2. Even assuming that the 4n state in the trap was really
	bound, the evolution into a continuum state involves
	branching at each threshold with discontinuity in the
	second derivative of energy with respect to a strength
	parameter. Polynomial extrapolations neglect this dis-
	continuity and are thus conceptually incorrect. The
	non-trivial evolution of the simple 2n case was con-
	sidered in detail, and is reproduced in our Appendix.
	Gandolfi et al. argued [137] that they had not claimed
	the existence of 3n or 4n resonances, nor seen evidence of
	open decay channels in the neutron trap, and that their
	extrapolation had worked in some selected examples.

  Waiting for the next theoretical paper claiming near-threshold 3n and 4n resonances, it seems clear that after all these decades of efforts and discussion a general consensus based on rigorous solutions of the problem is emerging to conclude... what was already clear from the first Glöckle's work in 1978.

  4,5,6 H isotopes are unbound by several MeV and exhibit broad resonances, but 5 H ground state is narrower than the other two, as demonstrated by recent ab initio calculations of these resonances[87]. Moreover,7 H is reported to lie almost at the 4n threshold[146], although it has in common with the tetraneutron that the experimental signals to date are weak, ambiguous and sometimes contradictory.

Note that there is a misprint in Table1of Ref.[97], the real part of the resonance energy appears with opposite sign.

It is indeed very disappointing that this trivial kinematical fact could have been ignored not only by the authors but also passed through the referee procedure of top level journals.

The same potential was considered in[110] to illustrate an S-wave resonance-like behavior.