Baryonic thermal screening mass at NLO

We determine the resummed 1-loop correction to a baryonic thermal screening mass. The calculation is carried out in the framework of a dimensionally reduced effective theory, where quarks are heavy fields due to their non-zero Matsubara frequencies. The correction due to interactions is computed at O($g^2_{ }$) in the coupling constant. In order to solve a 3-body Schr\"odinger equation, we exploit a two-dimensional generalization of the hyperspherical harmonics method. At electroweak scale temperatures, the NLO correction represents a $\sim 4.6 \%$ increase of the free-theory value $3\pi T$ of the screening mass.


Introduction
The dynamics of QCD at high temperatures plays a rôle in a number of physical processes, ranging from the cosmological evolution of the universe a few microseconds after the Big Bang, to the interpretation of empirical data from heavy ion collision experiments.Due to asymptotic freedom, one may hope that at high enough temperatures, the influence of interactions can be incorporated via a weak-coupling expansion.However, QCD effectively behaves as a three-dimensional Yang-Mills theory in this regime [1], which displays confinement.This implies that the coefficients appearing in the weak-coupling expansion need to be solved non-perturbatively, starting from some order.Sometimes, it is possible to go to a high order before non-perturbative coefficients appear, for instance O(g 6 T 4 ln(1/g)) for the equation of state [2,3].In other cases, such as the Debye screening mass associated with certain gluonic operators, only the level O(g 2 T ln(1/g)) can be reached perturbatively, before non-perturbative effects are met [4,5].
Apart from the appearance of non-perturbative coefficients, another issue with thermal perturbation theory is that the convergence of the weak-coupling expansion appears to be slow, up to very high temperatures.
In the present study, we focus on screening masses in the hadronic sector.Screening masses are the inverses of spatial correlation lengths, which describe how the quark-gluon plasma reacts when a state with given quantum numbers is put into the system.For hadronic observables, screening is very efficient, with a large mass of O(πT ) appearing at leading order.The next-to-leading order (NLO) correction is of O(g 2 T ), without any logarithm.This correction turns out to be perturbative, i.e. infrared (IR) finite, even though individual diagrams do contain IR divergences.
As usual, hadronic operators can be divided into mesonic and baryonic ones.The complete O(g 2 T ) correction to flavour non-singlet mesonic screening masses was determined a long time ago, in the framework of a dimensionally reduced effective theory [6].Recently, lattice calculations of the mesonic masses have reached an unprecedented accuracy over a wide range of temperatures, and allowed for a direct comparison between perturbative and non-perturbative data [7,8].Lattice data were found to be compatible with the O(g 2 T ) correction, even if fast convergence seems to be limited to temperatures which are well above the electroweak scale.
On the other hand, the baryonic sector has been investigated with less precision, both on the analytical and on the numerical side.The analytical result available in the literature is qualitative [9], and all the lattice calculations, both in the quenched approximation [10,11] and in the full theory [12,13], are restricted to very low temperatures, at which a direct comparison with analytical results seems difficult.Moreover, no continuum-limit extrapolation has been performed.Nonetheless, the steady progress in lattice measurements at very high temperatures [8,14,15] motivates a quantitative estimate of the NLO correction to baryonic screening masses, and this is the main goal of this paper.
Our presentation is organized as follows.In sec.2, we define the correlation functions we are interested in, and give the definition of screening masses, which characterize their asymptotic behaviour at large distances.Section 3 introduces the dimensionally reduced effective theory of QCD, as an efficient tool for resummed perturbative computations.In sec.4, the correlation functions from sec. 2 are re-expressed in the language of the effective theory, and they are shown to satisfy a two-dimensional Schrödinger equation.A numerical solution of the Schrödinger equation is presented in sec.5, before we turn to conclusions in sec.6.Further details on the effective theory, on the extraction of a static potential, and on the methods used for the numerical solution of the Schrödinger equation, are reported in four appendices.For QCD, we adopt the notation reported in appendix B of ref. [16] and, unless stated otherwise, refer to it for unexplained notation.

Preliminaries
We are interested in an interpolating operator which carries the nucleon quantum numbers.Maybe the simplest is (we consider N c = 3 colours throughout) where a, b, c are colour indices and C is the charge-conjugation operator, defined in appendix A.2.The contraction with the totally anti-symmetric symbol e abc guarantees gauge invariance.The two-point correlation functions considered are where k 0 = πT is the lowest positive fermionic Matsubara frequency, arising due to antiperiodic boundary conditions in the temporal direction, P ± = (1 ± γ 3 ) /2 is the x 3 -parity projector, the trace is over Dirac space, and r ≡ d 2 r denotes an integration over the transverse spatial directions.The corresponding screening masses, characterizing the longdistance behaviour of eq.(2.2), are defined as In vacuum, due to the spontaneous breaking of chiral symmetry, the positive (N , C + ) and the negative (N * , C − ) parity partners have masses which differ by several hundreds of MeV [17].On the contrary, at high temperatures, thanks to the effective restoration of chiral symmetry, the screening masses associated with the parity partners are expected to become degenerate (for numerical evidence, see refs.[14,18,19]).

Effective action
In QCD at very high temperatures, adopting the imaginary-time formalism, field fluctuations which depend on the temporal coordinate carry a large "mass", and should therefore represent weakly coupled short-distance physics.The strongly coupled long-distance degrees of freedom are given by the zero Matsubara modes of the gauge fields.Given that quarks obey antiperiodic boundary conditions over the temporal direction, they are always "heavy", even if their vacuum mass would vanish.
The gluonic sector of the effective theory contains gauge fields A k , with k = 1, 2, 3, living in three spatial dimensions, whose dynamics is non-perturbative [1].They are coupled to a massive scalar field A 0 , which transforms under the adjoint representation of the gauge group (cf.appendix A for a more detailed discussion).By taking into account these degrees of freedom, the corresponding effective action, which is usually called Electrostatic QCD (EQCD), reads [20,21] where the dots stand for higher-dimensional operators [22].Here i, j = 1, 2, 3 and [D i , D j ] = −ig E F ij with the covariant derivative defined as The matching coefficients m 2 E and g 2 E parametrize the Lagrangian mass squared of the scalar field A 0 and the dimensionful coupling constant of the three-dimensional Yang-Mills theory, respectively.They have been matched to QCD at several orders in perturbation theory, with the leading-order expressions reading [23][24][25] where g is the QCD coupling constant and N f is the number of flavours.
On the other hand, at high temperatures quarks are heavy with masses ∼ πT , due to fermionic Matsubara frequencies.The dynamics of such fields is described by threedimensional non-relativistic QCD (NRQCD) [26][27][28], see appendix A.2 for a derivation.By taking into account the power-counting rules in table 1 of appendix A.2, which imply that normal rather than covariant derivatives can be used in the transverse directions, the effective action for N f = 3 massless fermions in the lowest Matsubara sectors reads where χ and ϕ are two-component spinors defined in the paragraph below eq.(A.9), 1 f is a flavour index, and is a matching coefficient that was computed at 1-loop order in ref. [6].In eq.(3.3), according to the power-counting rules, we neglected higher-dimensional operators, see ref. [26] for some of them.Note that the action in eq.(3.3) displays more symmetries than the full theory, with e.g.χ f and ϕ f being independent of each other at this order.The QCD dynamics at high temperature is now described by S QCD 3 = S EQCD + S NRQCD .

Equations of motion
From the effective action in eq.(3.3) and an infinitesimal transformation of path-integration variables, it is straightforward to see that the three-dimensional fields χ and ϕ for each flavour 2 satisfy for a generic interpolating operator O(y) the equations of motion 1 While χf and χ f are independent variables in the path-integral formulation, in the canonical formulation they are related through χf = χ † f , and similarly for ϕ f . 2 The flavour index is omitted unless it is necessary for the clarity of presentation.
where the derivatives act on the x coordinates.Analogous equations hold for χ and φ.The propagators of χ and ϕ are defined as where in eq.(3.7) the expectation value ⟨•⟩ f indicates the path integral over fermions only.
Then by choosing O = χ, φ at y = 0 in eqs.(3.5) and (3.6), respectively, the fermion propagators satisfy the equations where 1 stands for the identity in spinor and colour indices.Since the fermions had been integrated out, the expectation values in eqs.(3.8) and (3.9) indicate the path integral over the gauge fields.Note that these equations are valid also without integrating over the gauge fields, i.e. for a fixed gauge field background, and that at this order the fermion propagators are diagonal in flavour and spin.

Perturbation theory
The free fermion action is obtained by setting g E = 0 in eq. ( 3.3).The corresponding equations of motion are readily worked out from eqs. (3.8) and (3.9).The free propagators can be written as (the coordinate-space expression is given in eq.(B.1)) where p ≡ d 2 p/(2π) 2 .At next-to-leading order in g E , we can define and analogously for S ϕ (r, x 3 ).By solving eqs.(3.8) and (3.9) and approximating the transverse movement (cf.appendix B), we obtain

S
(1) 4 Baryonic correlators in the effective theory

Contractions
The effective-theory expression for the baryonic interpolating operator from eq. (2.1) is readily obtained by using the definitions in appendix A.2.After taking the Fourier transform in the time direction, we may restrict to the contributions that involve the lowest Matsubara modes propagating in the positive x 3 -direction.This implies that N is represented by two χ fields and one ϕ field.Furthermore, by displacing the fundamental fields in the transverse direction, i.e. by introducing a point-splitting, the Fourier transform in the compact direction of the baryonic interpolating operator leads to where α is a two-component spinor index.To avoid clutter, we have omitted overall factors T3/2 from both operators, originating from eq. (A.7), however they are restored in eq. ( 4.4).
The two-point correlators from eq. ( 2.2) are defined in the effective theory as where 1/T comes from 1/T 0 dx 0 , and We have used the same symbols as in QCD, given that the ambiguity can be resolved from the context.
We remark that, as long as r i ̸ = r j , neither the interpolating operator above eq.( 4.1) nor the correlation function in eq. ( 4.3) is gauge invariant under gauge transformations involving the transverse coordinates.Gauge invariance could be restored by contracting the point-split operator with transverse Wilson lines.However, such transverse Wilson lines play no rôle in the calculation of the screening masses.Indeed, the final result will be gauge independent even without them, as gauge dependence vanishes in the largeseparation limit in the longitudinal direction. 3In order to streamline the presentation, we omit the transverse Wilson lines, and display results only for the Feynman gauge.
By exploiting the antisymmetry of the Levi-Civita symbol, and noting that the propagators are flavour independent, integration over the fermionic fields yields where the Wick contraction is defined as Equation (4.4) implies that our C ± is a sum of two independent correlation functions.This is a consequence of the fact that the action in (3.3) displays "emergent" global symmetries, with the numbers of χ and ϕ-particles separately conserved.Given that the two correlators in (4.4) differ just by a permutation of coordinates, and that in the end all coordinates are set equal (cf.eq.(4.2)), the two correlators yield the same baryonic screening mass.This degeneracy could be broken by higher-dimensional operators in the effective theory [26], leading to a "fine structure" of the screening spectrum, but this is an effect of higher order than our O(g 2 T ).

Free limit
By inserting the free propagators from eqs. (3.10) and (3.11) into eq.(4.5), we see that From here it follows that W (0) satisfies a (2+1)-dimensional Schrödinger equation, Thus, if quarks have small transverse momentum (indeed we will see that parametrically ), the exponential falloff is dominated by 3M = 3πT + O(g 2 ), which then represents the leading contribution to the baryonic screening masses m ± .

Next-to-leading order
To date, the only estimate of O(g 2 ) corrections to a baryonic screening mass in the high temperature regime of QCD is qualitative [9].Here the full O(g 2 ) correction is derived in the same way as for the mesonic case [6], demonstrating in particular its IR finiteness up to this order in the weak-coupling expansion.
The equation of motion for W in the interacting case is readily worked out from eqs. (3.8) and (3.9), and it reads where we have introduced r i ≡ (r i , x 3 ) to simplify the notation.Inserting S χ and S ϕ from eqs. (3.13) and (3.14), respectively, and performing the gluon contractions (cf.appendix C for the gluon propagator and further details on intermediate steps), we get Here, with the notation from eq. (C.4), To extract the screening masses, we take the limit x 3 → ∞, which leads to where r ij ≡ |r i − r j |, and V ± are the static potentials defined in ref. [29], V ± (r) ≡ 4 3 where γ E is the Euler-Mascheroni constant and K 0 is a modified Bessel function.It is appropriate to stress that according to eq. (4.11), the three-body potential receives contributions from two-body interactions only.Finally, by replacing W (0) → ⟨W ⟩ on the right-hand side of eq.(4.9), which is justified at O(g 2 E ) and implements a resummation of potential-like interactions, and taking the already-mentioned limit x 3 → ∞, the equation of motion for the correlator reads where The two contributions in eq. ( 4.4) satisfy the same equation, just with a permutation of coordinates (which in the end are set the same, cf.eq.(4.2)).Therefore, the solutions of both equations yield the same screening mass, which is then also the screening mass extracted from C ± (x 3 ) at O(g 2 E ).We end this section by remarking that the potential U (r 1 , r 2 , r 3 ) from eq. (4.11) is symmetric in the exchange of its first and last coordinate.For the first contribution from eq. (4.4), this corresponds to r 1 ↔ r 3 , which in terms of eq.(4.1) is due to an accidental symmetry in the action from the exchange χ u ↔ χ d .In contrast, for the second contribution from eq. (4.4), this corresponds to r 2 ↔ r 3 , the exchange of two identical χ d particles.

Center-of-mass coordinates
From the discussion in sec.4.3, for large separations in the longitudinal direction, the equation of motion for a generic two-point correlation function related to baryonic interpolating operators (in both parity channels) implies the Schrödinger equation where the potential is given by eq.(4.14).The energy eigenvalue of the ground state yields our screening masses, i.e. m ± = min{E} + O(g 3 T ).
In order to solve the three-body Schrödinger equation, we employ the Jacobi coordinates where R is the position of the center-of-mass in the transverse plane, ξ 1 is the relative separation between two quarks located at r 1 and r 3 , while ξ 2 describes, up to some numerical factor, the relative separation between the quark in r 2 and the center-of-mass of the other pair.In this sense the set of coordinates (ξ 1 , ξ 2 ) describes the relative separations of the underlying two-body problems.With this change of variables the potential only depends on ξ 1 and ξ 2 , given that (5.3) In this way, the Laplace operator can be separated into the center-of-mass and relative motions.The Schrödinger equation for the relative motion can be written as (cf.eq.(D.2)) where it is understood that the static potential and the wave function are expressed in terms of the new coordinates ξ 1 and ξ 2 .

Numerical solution
In order to find a numerical solution to eq. ( 5.4), it is convenient to define the dimensionless transverse coordinates (5.5) Moreover, we express E in terms of a dimensionless eigenvalue Ê, by writing This leads to a Schrödinger equation in terms of dimensionless variables, where V is a rescaled static potential from eqs. (4.11) and (4.12), and ρ is a re-parametrization of the dimensionful quantities of the problem [6], (5.9) Equation (5.7) can be solved numerically by exploiting a two-dimensional generalization of the so-called hyperspherical harmonics method.It is usually employed for threedimensional quantum many-body problems, see ref. [30] for an introduction and appendix D for the two-dimensional generalization.The idea is to reduce the eigenvalue problem in eq.(5.7), which depends on four independent radial variables, to a set of coupled differential equations, which depend on one radial and three angular variables.This leads to where ξ, defined in eq.(D.8), is the radial variable, and L is a non-negative integer number.
The matrix elements V {L,L ′ } are defined in eq.(D.25).They are computed in the basis provided by the eigenvectors of the hyperangular momentum operator, defined in eq.(D.10).Let us stress that those eigenvectors are known analytically, see eq. (D.12), and therefore the matrix elements V {L,L ′ } can be determined numerically with a moderate computational effort.Since we are looking for the ground state of the system, it is natural to restrict ourselves to the case of zero total angular momentum (cf. the discussion above eq.(D.26)).This reduces the calculation of the matrix elements to states associated with even values of L only.
For the numerical evaluation, we do need to restrict to a finite number of eigenstates of the hyperangular momentum, i.e. truncate the basis up to a maximal value L max .The choice of L max is critical and there is no general prescription for this choice, which may be, in general, potential-dependent.
In fig. 1, we show the ground-state eigenvalue as a function of L max .As expected, by increasing the number of states, the result becomes increasingly stable.The relative difference between the lowest eigenvalue extracted with L max = 28 and with L max = 30 is ∼ 10 −6 .The final estimate for the ground-state eigenvalue for N f = 3, is therefore As a further cross-check, we have discretized the Hamiltonian on a mesh grid.The two relative positions from eq. (5.5) live in Cartesian dimensions of extent N a each, where a is a lattice spacing (in units of m −1 E ).The wave functions are then N 4 -component vectors, and the Hamiltonian is an N 4 × N 4 -matrix.For the Laplacian a 5th order discretization is adopted.We extrapolate separately a → 0 with N a fixed (continuum limit), and N a → ∞ (infinite-volume limit), checking that the end result is independent of the ordering of these limits.Including values in the ranges N = 20...64 and am E = 0.16...0.75, we find that the results are consistent with an approximately quadratic dependence on a.The continuum limit yields a result compatible with eq.(5.11), up to all digits shown.Inserting eq. ( 5.11) into eq.( 5.6), leads to our NLO estimate of the baryonic screening masses, min{E} (5.12) The NLO term represents a ∼ 4.6% positive correction to the free-theory result 3πT at the electroweak scale, i.e. g 2 ≈ 1 [24].

Conclusions
We have computed the NLO correction to a baryonic screening mass in thermal QCD.
The computation was carried out in the framework of a dimensionally reduced effective theory, where the lowest Matsubara modes of massless quarks appear as non-relativistic fields, interacting with zero Matsubara mode gluons.The NLO correction in eq. ( 5.12) amounts to a ∼ 4.6% increase of the free-theory value 3πT at g 2 ≈ 1.For comparison, we note that an analogous computation for the mesonic screening masses yields a ∼ 3% positive correction in the static sector [6], while in the non-static sector the effect is ∼ 5% [29].The reason for the difference lies in the form of the static potential in the corresponding Schrödinger equation.The larger correction in the baryonic and non-static mesonic sectors can be traced back to the appearance of the potential V + in eq.(4.14), which is more energetic than V − at short distances.
The computation of our paper was motivated by, and in turn provides motivation for, a non-perturbative determination of baryonic screening masses on the lattice [14,15].By comparing the two results, we can attempt estimate the domain in which tations can function as a quantitative tool for understanding the underlying physics.As discussed in more detail in ref. [15], the upshot is that, assuming that the coefficients of the higher-order corrections are small, the O(g 2 ) term can account for ∼ 90% of the difference between the non-perturbative result and the free-theory value 3πT down to temperatures of a few GeV.This suggests that baryonic observables could be more perturbative than most other quantities that have been studied so far.Of course, to consolidate this picture, it would be valuable to explicitly determine the coefficients of some higher-order corrections.field/operator power counting Power counting for the fermionic fields in eq. ( 3.3) and for their variation in longitudinal (x 3 ) and transverse (r = x ⊥ ) directions.
solving the equations of motion for ψ ↑ 0 and ψ ↓ −1 at O(g 2 E ), we obtain the equations of motion in eqs.(3.5) and (3.6), where we have dubbed χ = ψ ↑ 0 and ϕ = ψ ↓ −1 , respectively.In order to determine NLO corrections to hadronic correlation functions, it is helpful to establish the power counting deriving from the three-dimensional action.The relevant scales of the system are m E , g 2 E and the thermal quark mass k 0 ≈ M , originating from the lowest fermionic Matsubara frequency.The power-counting rules are given in table 1.As a consequence, the action in eq. ( 3.3) is obtained by taking into account spinor fields in the lowest Matsubara sector, including terms up to O(g 0 E ), and labelling, at variance with appendix B of ref. [16], the N f = 3 flavours as f = u, d, s.

B Fermionic propagators in the effective theory
In coordinate space, by making explicit the source and the sink locations, the propagator from eq. (3.10) becomes and similarly but with opposite sign for ϕ.At the next order, from eq. (3.8), the propagator for χ reads implying that χ propagates freely from the source to the position (z, z 3 ), where it emits a gluon, and then again freely propagates up to the sink position.By inserting eq.(B.1) into eq.(B.2) and taking (y, y 3 ) = (0, 0), the result can be rewritten as S (1)  χ (x, x 3 ; 0, 0) = S (0) χ (x, x 3 ; 0, 0) 4 For the notation we recall footnote 1.
From table 1 we see that the quadratic dependence is parametrically of order πT x 3 /[(x 3 −z 3 )z 3 ] ∼ πT g 2 E , implying that fermions probe the variation of gauge fields at distances |z − z 3 x/x 3 | ∼ 1/(gT ) ∼ 1/m E .For the static potential, we need the x 3 → ∞ limit, and then only the position of the gauge field at z 3 ∼ x 3 matters.Then the prefactor is large, and we may evaluate the integral in the saddle point approximation, yielding S (1)  χ (x, x 3 ; 0, 0) ≃ S (0) χ (x, x 3 ; 0, 0) justifying eq.(3.13).

C Evaluation of the static potential
We show here details for the manipulation of eq.(4.8).The free-theory gauge field propagator reads where µ, ν = 0, . . ., 3, a, b = 1, . . ., 8 and, in Feynman gauge, The Debye mass m E is given in eq.(3.2), and λ is an auxiliary mass parameter, used as an IR regulator.The Wick contractions following from eq. (4.8) lead to the colour algebra where T d are Hermitean generators of SU(3), normalized as tr T c T d = δ cd /2.Subsequently, we are left over with the time-integrals (C.4) The propagators are written in a Fourier representation, like in eq.(C.2).In the limit where P denotes a principal value, and pulling the time integral inside the Fourier representation, the Fourier transform becomes 3 ) Carrying out the two-dimensional momentum integral, the potential reads lim where r 12 ≡ |r 1 − r 2 |, γ Euler-Mascheroni constant, and K 0 a modified Bessel function.The combinations that appear in eq.(4.10) are defined as By using the expansion of K 0 (x) for small x, it is straightforward to show that lim This then leads to the explicit expressions for V ± (r) reported in eq.(4.12).

D Two-dimensional hyperspherical harmonics
Let us consider a system of three particles of mass m which are constrained to move in a two-dimensional plane.We assume those particles to interact through mutual two-body potentials which depend only on relative separations.The Schrödinger equation for such a system is where r ij ≡ |r i − r j |.By performing the change of variables in eq.(5.2), the Laplace operator can be written as where R is the center-of-mass coordinate.In such a way, the kinetic energy can be separated into center-of-mass and relative motions.The Schrödinger equation for the relative motion can be written in terms of the new coordinates as By going to polar coordinates, i.e.
where ξ i refers to the absolute value of the 2-vector ξ i and the superscript (j) refers to the components in the two-dimensional plane, the Laplace operator can be expressed as Both for ξ 1 and ξ 2 the angular to the square of the twodimensional quantum-mechanical angular momentum operator, We also note that the scalar products following from eqs. (5.2) and (5.3) read where So far, our coordinates comprise two radial extension ξ i and two angles θ i .The idea underlying the hyperspherical harmonics method is to write the system, by a suitable change of variables, in terms of one radial and three angular variables.This is achieved by going to the so-called hyperspherical coordinates, where ξ is called the hyper-radius and the angular variable ϕ describes the relative length between the two vectors ξ 1 and ξ 2 , i.e. ϕ ∈ [0, π/2].Together with θ 1 and θ 2 , ϕ forms a set of so-called hyper-angles Ω 4 ≡ {θ 1 , θ 2 , ϕ}.With this change of variables, the Laplace operator in eq.(D.5) becomes where we introduced the hyperangular momentum operator In complete analogy with the angular momentum operator, we can define the eigenvalues and eigenstates of the hyperangular momentum operator by assuming that ξ L Y (Ω 4 ) is a harmonic function.This leads to where L, in analogy with the angular momentum quantum number, is called the grand orbital quantum number.Given the form of the operator in eq.(D.10), the corresponding eigenvectors can be written as where Y l i (θ) are the eigenvectors of the L i operators, i.e.
By looking for a solution of the type and by performing the change of variables x = cos 2ϕ, the differential equation in eq.(D.11) can be written as where n is related to the quantum numbers l 1 , l 2 and L by Equation (D.15) has analytical solutions in terms of Jacobi polynomials, where (•) n is the Pochhammer symbol and 2 F 1 is the hypergeometric function.These polynomials are orthogonal in the interval [−1, 1] with respect to the weight function , and their normalization constant is given by The hyperspherical harmonics form a complete set on the surface described by the angles {ϕ, θ 1 , θ 2 } at fixed ξ.As a consequence, the wave functions appearing in eq.(D.where {L} refers to any particular combination of quantum numbers which leads to L = 2n + |l 1 2 |, and the expansion coefficients Φ {L} encode the hyper-radial dependence of the wave function.We can get rid of the first derivative appearing in eq.(D.9) by writing the expansion coefficient as Φ {L} (ξ) = u {L} (ξ) /ξ 3/2 .Therefore, by using eq.(D.11), the Schrödinger equation can be written as are the matrix elements computed in the basis of the hyperspherical harmonics.In this way, the Schrödinger equation has been reduced to a set of coupled onedimensional differential equations which are easier to solve numerically.With this procedure, the bulk of the calculation is given by the computation of the matrix elements by properly truncating the basis up to a certain value L max .The number of states to be included can be reduced by scrutinizing the physics of the problem of interest.Given that our potential V only depends on θ 1 − θ 2 (cf.eq.(D.7)), the total angular momentum in the transverse plane, with its eigenvalue l 1 + l 2 , commutes with the Hamiltonian and represents a conserved quantity.In order to extract the energy of the ground state, it is then sufficient to consider the smallest value of the total angular momentum, i.e. set l 2 = −l 1 in eq.(D. 16).This implies that the grand orbital quantum number is restricted to even values only, L = 2 (n + |l 1 |) . (D.26) Furthermore, given that the potential is invariant under θ 1 → θ 1 + π, we can restrict to wave functions symmetric under this "parity", and choose l 1 even.In general, any fixed L describes a set {L} of various n and l 1 .

Fi gure 1 .
Dependence of the ground-state eigenvalue Ê (open blue circles) on the basis size L max .The calculation is restricted to even values of L, in order to take into account eigenstates with zero total angular momentum.The dashed red line represents the infinite-L max extrapolation.