Calculation of Dynamical Response Functions Using a Bound-State Method

We investigate a method to extract response functions (dynamical polarisabilities) directly from a bound-state approach applied to calculations of perturbation-induced reactions. The use of a square-integrable basis leads to a response in the form of a sum of δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} functions. We integrate this over energy and fit a smooth function to the resulting stepwise-continuous one. Its derivative gives the final approximation to the physical response function. We show that the method reproduces analytical results where known, and analyse the details for a variety of models. We apply it to some simple models, using the stochastic variational method as the numerical method. Albeit we find that this approach, and other numerical techniques, have some difficulties with the threshold behavior in coupled-channel problems with multiple thresholds, its stochastic nature allows us to extract robust results even for such cases.


I. INTRODUCTION
Various methods have been developed to describe total and differential cross sections via square-integrable wave functions [1][2][3][4][5].Another common approach in nuclear physics is to use the response functions (dynamical polarisabilities) in the calculation of cross sections [6].As far as we are aware, the first method that combined these two ideas was the Lorentz-integral transform (LIT) [2,[7][8][9][10][11][12][13][14][15][16][17].This uses a transformation of the response functions that leads to the solution of a set of inhomogeneous Schrödinger equations.A suitably weighted sum of overlaps between the solutions is then subjected to an inverse transform to obtain the response or dynamical polarisabilities.
This procedure is reasonably robust for problems with a single open channel, such as the photo-dissociation of the deuteron [14], but seems more challenging in other cases (but see [9,13,18]).Most standard nuclearstructure methods can be employed to generate a basis and, as we shall see below, good answers can be obtained with very different approaches in cases with a single open channel.For more than one open channel, however, all methods seem to struggle to some extent, something that we have found to be more obvious in those using the Stochastic Variational Method (SVM) [19,20].
Bound-state methods work with square-integrable basis functions, and any response calculated is naturally discrete.In the manner we apply this, this means that any response function becomes a sum of δ distributions with varying strengths, located at the eigenenergies of the channel Hamiltonian.In the LIT, this sum of δs gets folded with a finite-width Lorentzian [7], giving a continuous function with a small but non-zero width.An inverse transform is then applied to this to get the physical response function.Typically, this inversion is done by fitting to the LIT of a continuous function constructed as a sum of carefully chosen basis functions [9].These basis functions should reflect properties expected of the physical, continuum response, such as the correct threshold behaviour.The choice of width also requires care: if too small, the discrete δs start to be resolved; but if too large, physical features can get washed out.
Despite these subtleties, indications from the literature are that the method can work well for structured bases, such as hyperspherical harmonics (HH) or no-core-shellmodel (NCSM) states [2,[7][8][9][10][11][12][13][14][15][16][17].This type of basis leads to a regularly-spaced spectrum of the kinetic energy operator.In contrast, the stochastically chosen basis of the SVM leads to a rather randomly-spaced spectrum, which can generate additional unphysical structures in the transformed response and so can be less well-suited for inversion of the LIT.
In this paper, we develop a powerful alternative which actually takes advantage of the randomness of the SVM.Instead of folding the response with a smearing function, this uses the integrated response, that is, the integral up to some energy ω of a response function obtained with a square-integrable basis.This turns the sum of δ distributions into a step-wise continuous function of ω.We then fit a continuous function to this, and differentiate to get an approximation to the physical response function.One could argue that this has just replaced one difficult problem (robust inversion of the LIT) to another (fitting of a function that is robust enough to yield a reliable derivative), but we will find that the latter is often an easier one.In this paper, we outline the method and present evidence that it is trustworthy, by comparing it with both an analytically solveable model and results from the LIT.
The paper is organised as follows: in Sec.II we succinctly set out the underlying definitions; most of these can be found elsewhere in detail.We then, in Sec.III show that for a simple Pöschl-Teller potential, all the required quantities can be calculated in analytic form for the full continuum calculation, thus providing an important benchmark for basis-based methods.In the next section, Sec.III B, we solve the problem using first a simple harmonic oscillator basis, and then a suitable Gaußian basis.The latter will be shown to be extremely efficient.These are still relatively trivial problems, and we next compare to the photo-disintegration of the deuteron as discussed by Bampa et al. [14].We see that both using the LIT transforms, and using a simple fit to an integrated response function, give largely identical results, apart from small differences in areas of small response.
Next, we show that our approach well reproduces the response in a 3-particle continuum, Sec.V, in anticipation of the difficulties we will encounter when we turn to a combination of a single particle continuum relative to a two-body bound state and a three-particle continuum in Sec.VI.We compare a hyperspherical-basis calculation and an SVM result for a two-channel model which has a two-body bound state in one of the channels.
Finally, we draw conclusions and give an outlook in VII.

II. THE CALCULATION OF THE CROSS SECTION
Inclusive cross sections due to an external probe (also called "perturbation-induced") are typically of the form [6] where , q with q = |q| are the energy and momentum transfer of the external probe to the system, g is a coupling constant, Ω = θ, φ are the scattering angles, and the f i are kinematic factors.
The key aspect of the calculation is thus the dynamical response functions [10,12].For a single channel with no additional bound states, these read The complete set of scattering eigenstates ψ k are normalised as As we shall see below, this can be calculated either directly, or by a Lorentz integral transform where we take without loss of generality σ I > 0 but do not restrict σ R .Using completeness, Eq. ( 3), we can also express the Lorentz transform using the solution of an inhomogeneous Schrödinger equation, In a finite basis of square-integrable functions, we can write the LIT decomposition in terms of the eigenfunctions φ i and eigenvalues E i of the matrix representation of the Hamiltonian as discussed in [12], i.e., with Thus there is an explicit solution, which does not rely on the transform, as a sum of δ functions Clearly this form is not useful in itself, since we know that the solution obtained by using the scattering states is continuous.However, completeness shows that we can impose the the sum rule which is true for any complete basis, L 2 or not.This is useful to constrain the choice of basis functions.
Ref. [12] suggests that increasing the number of basis functions leads to a denser spectrum, and thus a better LIT transformation; this is actually quite a subtle issue, as shown below.In the standard works on the LIT, one usually uses a nonzero σ I to smooth out small fluctuations when inverting the LIT.This removes detail from the response function, and we expect that there is a price to pay in the inversion-at a minimum the displacement of reaction strength below the reaction threshold.That can be a small price if the response is smooth and phasespace factors suppress the cross section at low energies.Nevertheless, we find it difficult to control the smoothness of numerical results, especially when using Gaußian basis sets.In that case, it is much better to take some form of the limit σ I → 0, and here we will discuss one practical approach to approximate this limit.To this end, we shall mainly investigate the integrated response function rather than F .In the next section, we illustrate this approach for a simple one-dimensional potential.
Naively, we see there is second state at exactly zero energy for λ = 1, which is one of the cases we will consider.This is not a problem since the potential is transparent and has no resonances in the continuum.However, it may be linked to some of the properties of the analytical solution below.We now consider the single bound state case for λ = 1, which has a bound state energy −1/2 with normalised eigenfunction We shall look at the "photo-excitation" cross section for this model potential in the dipole approximation, O(q) = Eqx, and for simplicity we shall choose units where the product of electric field and charge mathcalEq = 1.
In order to evaluate Eq. ( 2) directly, we note that the operator x creates a state in the positive energy continuum.This means that we need the odd-parity positiveenergy spectrum for energy E = k 2 /2, which is known to be [23] We see that these wave functions approach sin(kx ± δ k /2) for x → ±∞, with the phase shift δ k = 2 arctan(1/k).

Momentum-space basis
From the result above, we calculate the generalised Fourier-decomposition of xψ 0 relative to the energy It is quite illustrative to calculate the inverse generalised Fourier transform of ψ 1 (k) in detail, This integral can be tackled by contour integration.Closing the contour in the lower half plane [24], we have contributions from the poles at k = −(2n + 1)i, n ∈ N.
Taking the residues, we recover the original function as expected:

Consistent normalisation
We have not explicitly shown a consistent normalisation of the momentum space states; however since the integral equals the corresponding integral in k-space we see that the normalisation is consistent.

Response function
The benchmark for all future calculations is the explicit analytic result for the response function, which from Eq. (2) becomes where ω = k 2 /2, which is the energy relative to threshold.This has an inverse square-root singularity at threshold, which is due to the constant density of states in one dimension.The integrated response distribution is then  19) (red curve), and the integrated response T , Eq. ( 20) (blue curve).The gray dashed line is the saturation value π 2 /12.see Fig. 1 for a representation of the integral and integrand.The most important message from that graph is the saturation of the integrated response function, which can be easily be proven to equal the norm squared of x |ψ 0 , To understand the analytical form of the LIT transformation, we express the response function as a k space integral, where There are two classes of contributions to the integral, namely one from the pole of the Lorentzian and the other from the poles of the sech 2 , namely the poles of the wave function which lie in the lower half plane at k = −(2n+1)i due to the specific contour chosen, see example above.The latter contribution is given by the sum where ψ (n) is the Polygamma function.The shift of σ by 1/2 should be interpreted as the effect of the bound-state pole at E = −1/2.The poles of the Green's functions give a contribution to the integral in (22) of which is singular in the limit σ I → 0.
Both L 1 and L 2 contain a large contribution from the pole at σ = −1/2; these cancel substantially for small σ, and completely when σ I → 0. Thus the response function can then be recovered from L 2 ; indeed lim Moving into the complex σ plane does therefore lead to a slightly confusing situation: The part of the response function that contributes to the result on the real axis is suppressed as we increase the imaginary part of σ, and a totally separate expression dominates for large σ I .We see that the dominant contributions are caused by the fact that the wave function has poles in the complex k plane.The importance of the contribution of these poles of the wavefunction is in all likelihood specific to the Pöschl-Teller potential, but the fact that the effect of the Green's function dominates on the real axis is probably generic and would explain why using an L 2 basis at small but nonzero σ I is normally effective in practical calculations.

B. L 2 -basis calculations
In practical implementations of the methods discussed before, we normally use a square-integrable basis to find either the LIT transform, or the calculation of the integrated response function as advocated in this work.
We first look at the LIT transform where we solve the equation ( 5), where we now define σ relative to the threshold at zero energy, using a basis of L 2 functions.We shall use either a truncated complete basis consisting of a set of odd harmonic oscillator states, where a geometric progression with ratio α of widths starting at s is used to ensure a numerically acceptable condition number of the overlap (norm) matrix.
For the orthonormal basis (27) we get which gives the equation If we denote the eigenvectors of H in this basis as e n , and the as eigenvalues λ n , we get a LIT transform of the form (6), with γ i = e i • d.
For the Gaußian basis we define a Hamiltonian and overlap matrix by integration from the left with the same basis, The integrated response thus always approaches a limit for large energies, where the last inequality is saturated for a complete decomposition of the state, e.g., n → ∞.
Let us first analyse the LIT transformation, see Fig. 2 for some representative results using a harmonic oscillator basis.Since the Lorentzian ((x − σ R ) 2 + σ 2 I ) −1 is not normalised, the LIT transform diverges in the limit σ I = 0. Thus, for ease of comparison, it is better to look at the LIT transform of the normalised Lorentzian, L(σ) = L(σ)σ I /π.We see that for a coarse spectral spacing (small b) we get a very oscillatory LIT transformation, from which it would be impossible to reconstruct the continuum result for the response function.For smaller spectral spacings we get a convergent result.However, the smaller the value of σ I , the narrower the spacing we need in the spectrum to avoid spurious peaks and oscillations.That raises the question how we can extract the inverse reliably from just one of these transforms.We thus need to answer two questions: 1. How can one ensure a smooth inversion?We already know that the exact inversion of the data is given by a sum of delta distributions at the eigenenergies of the Hamiltonian (30), so for us to obtain a smooth result we can only perform an approximate inversion from a transform with substantial imaginary σ I .
2. How do we disentangle the singular and regular parts of the LIT?We have seen above that we have two contributions, but only one survives in the limit σ I → 0. At finite σ I both contribute, with domination from the "wrong" component at positive energies.
The second problem may be solved by doing the inversion at multiple values of σ I , and interpolating to σ I = 0.This is a non-trivial task in realistic calculations, since oscillatory noise always enters the results.As explained before, some of this may be special for the Pöschl-Teller potential-but if we have problems for one case, there may well be others.
The first one is more troublesome.In some sense what we are trying to do is to invert the LIT by sub-sampling, and solving the result by a continuous approximation, often by choosing a small set of basis functions that impose a shape on the resulting response function.Both of these are very indirect approaches, and it is almost unclear why we need the LIT as an intermediate if the inversion is such a difficult problem to resolve.
So, as an alternative approach, we consider the direct calculation of the response functions, given therefore simply as finite sum of delta distributions.It looks like we loose more than we gain.However, we now analyse the integral of this sum as a function of the "energy" ω = σ R which is continuous, and can reliably be approximated by a smooth function.This integrated response distribution T is continuous and approximated by a smooth differentiable function, which can be used to reconstruct the response function F itself.
We now analyse the results for the two basis sets.We have performed numerical calculations for harmonic oscillator bases with b = 1/2, 1, 2, 4 and 8, see Fig. 3. Since for the numerical calculations the integrated response T makes a finite jump (d • e n ) 2 at eigenvalue λ n , there are multiple ways to represent this.We chose to plot at λ n the midway point between these two values.This seems to be overall a very true representation of the analytical result.As we can see from Fig. 3, we can closely mirror the analytic results.In c), we also show bands obtained from the outer and inner steps; we see a large width for a low spectral density, and a much narrower result for a high one: the width goes down rapidly as b increases.Also, the midpoints track the analytical results remarkably well.
Since we want to yuse the SVM approach, which relies on Gaußian basis functions, later on, we are particularly interested in the efficacy of using Gaußian states as in Fig. 4. We see that in all cases the midpoints of the jumps track the analytical results extremely closely.If we bracket the solutions with the upper and lower values of the jump, we see in the lower panel b) that the width is rather large.That is clearly an enormous over-estimate of the theory error in the results, which are extremely close to the analytical results.We find that 20 Gaußian FIG. 4. The convergence of the cumulative response for a small Gaußian basis (28).In a) we show the results for α = 2 and s = 0.01 (data set i, 12 basis functions), s = 0.001 (data set ii, 16 basis functions) and s = 0.0001 (data set iii, 20 basis functions).As in Fig. 3, the points plotted are the midpoints of the jumps.In b) the lower and upper lines of each colour connect the values just before and just after the jumps, respectively.Colours correspond to those used in a).
basis states can achieve as much accuracy as ≈ 200 orthonormal ones.
So far, we were concerned with a slightly artificial onedimensional world.To show that the method set out above can be effective in a three-dimensional world and a much larger density of states above threshold, we consider the dipole excitation from an s-wave state, chosen as originating from a Pöschl Teller potential, where we choose λ = 1.8 for a single bound state at dimensionless energy −0.32 and wave function We now perform the calculation of the response function as set out above, with a perturbing dipole operator z.We diagonalise the p-wave Hamiltonian in a set of Gaußian basis functions z exp −αr 2 , calculate the integrated response function, and finally fit with a Padé approximant of the form where we take n = 4 in the results shown in Fig. 5.This provides an almost perfect fit.
We use two rather different basis sets but notice little difference in the results: the response functions, as shown in b) coincide within the line width.

IV. DEUTERON PHOTO-DISINTEGRATION
As a further test of our method in a more physical context, we analysed the onebody part of the total E1 deuteron-photodisintegration cross section using the Argonne V18 potential [25].This exactly matches the treatment using the LIT in Ref. [14], allowing us to make direct comparisons with that work.FIG. 6.A comparison of the strength function multiplied with energy-transfer squared, Fig. 8 in Ref. [14].Black line: J = 0, blue line J = 1 and red line J = 2.The dashed lines are the results from Ref. [14], solid lines are our results.
We use a numerical integration by mapping the interval 0 ≤ r < ∞ onto x ∈ [0, 1] via r = βx/(1 − x) and employ Gauß-Laguerre interpolation for a high-order approximation to the derivative, using the fact that we can continue the function as odd beyond x = 0 and x = 1.
As before, we calculate the cumulative response function F by integrating over the δ functions at the eigenenergies of the intermediate channel Hamiltonian.The resulting strength functions F J for intermediate-state total angular momentum J are shown in Fig. 6.We see that our results very closely track the LIT ones from Ref. [14].Small differences should come as no surprise, especially since the representation used in the figure enhances the small tails of these functions, where we expect the methods to show most the largest divergence.Bampa et al. devote substantial work to the LIT inversion, and still see clear dependence on the basis function used.In our method, the fit is relatively robust, and again not as sensitive to the tail since it is very flat in the integrated function.Nevertheless, it is encouraging to see that the differences appear rather small.
A more direct calculation of the basis of the cross section is given by the polarisation functions P J .For a deuteron, the E1 operator allows for the channels J = 0, 1, 2. The imaginary part is directly related to the F J 's, and the real part can be calculated from a dispersion relation.This is a much more robust test of differences between the strength functions.Once again, we see in Fig. 7 that the results of the LIT and of our method are very similar.The main conclusion from the results on the left hand-side, which are direct calculation, is that there are small differences only-essentially, from top to bottom to middle we zoom in by a factor of 10 each time, and even at the largest scale the deviations are small.Some of that seems to be due to a very robust calculation on our end, which explains the differences around 40 MeV.The real surprise is how our simple fit seems to more easily capture the large-energy behaviour which required substantial work in the LIT calculation.That can also be seen in the right column.Since each of these curves is the result of a dispersion integral, it shows how close our results are to the LIT ones over a much larger domain.This gives us some confidence moving forward, even though we have one additional problem to discuss.
We can use this output to calculate the differential photo-dissociation cross section, but little detailed data is available.Instead, we show in Fig. 8 the total photodisintegration cross section.The results from our calculation agrees well with the experimental data extracted from EXFOR [26].As we can see, the agreement is excellent, which is a confirmation both of the stability of the calculation and of the fact that the differences to the LIT variant are small.There is some indication that our results decay slightly too fast at the highest energies.That might not be surprising since our computation only contains the one-body E1 operator.While the importance of magnetic multipoles is well-known to decrease with energy, higher electric multipoles and pion-exchange currents play a more prominent at higher energies.

V. THREE-PARTICLE DECAYS
Clearly the method set out above is very effective for simple structureless two-body problems.One of the questions we have to answer is how it survives in the many-body context, especially where we may have multiple thresholds that complicate the calculations and their interpretation-indeed, this is one of the problems that led to the current approach.The issue already arises for A = 3, where both two-body and three-body breakup channels with different thresholds exist, e.g. 3 He→ ppn or → pd.One issue is the degeneracy of the three-body continuum: Since it is described by two Jacobi coordinates, we expect a substantial degeneracy as a function of the excitation energy, which is apt to lead to some complications.
So let us again study a simplified, exactly solvable, problem which, as we shall show later, is an illuminating simplification of the problem to be studied in detail in future articles.We consider the dipole response function for a simple three-body ground state consisting of two identical bosons interacting with a third distinguishable particle, all three of the same mass.The dipole operator either acts on the third particle, or equivalently on the two bosons: translation invariance shows these two are identical up to a multiplicative constant.This is obviously the bosonic equivalent to the 3 He system mentioned above.
To find a problem that can both be tackled analytically FIG. 7. A comparison of the polarisabilities that enter the calculation of the cross section, Fig. 9 in Ref. [14].Red dashed lines show results from Ref. [14], solid blue lines are results using the present method.We have chosen the scale on the axes to be identical to that in the reference cited, and thus in the figures in the left column the x-axis runs from 0 to 95 MeV, and on the right from 10 to 100.FIG. 8.A comparison of the deuteron cross section as extracted from our method (using the imaginary part of the responses in Fig. 7) versus the data from the EXFOR database [26] (green dots).
and numerically we start from the wave function (see (A2) for the definition of coordinates) and act with the dipole operator We assume that the Hamiltonian describing states in the L = 1 channel is just the free Hamiltonian, which takes the form This can be tackled with results discussed in Appendix A. We can evaluate the effect of the E1 operator using only the "12-3" Jacobi coordinates.We once again apply Eq. ( 5).Looking at the right hand side, we see the effect of the E1 operator on |ψ 0 gives an L = 0 solution in ξ, and a L = 1 state in η, The total response is the norm of this state, In order to find the response function, we need to find the relevant (normalised) continuum solutions, which are The energy for the direct product of these two states is The normalisation can be checked from: Thus This satisfies the consistency check If we now substitute q 2 ξ = ω − 3 4 q 2 η , and integrate over q η , we get, using Indeed, this still integrates to the same value as before.
The calculation for the integrated response can easily be done numerically, and we conclude that the function F rises as ω 3 for small ω, and thus with a quartic power for the integrated response, and decays for large ω as exp(−2ω/3α)ω 3 /2.
As can be seen in Fig. 9, there is a small but definite scatter in the fully numerical (SVM with a suitably random choice of basis functions) calculations, which agrees very well with the analytical result.Thus we have reached stage one: we are confident that we can deal with a three-body continuum correctly.Now we look at coupled channels.

VI. TWO-BODY BOUND STATE AND 3 BODY CONTINUUM
We now increase complexity and add to the three-body system a set of scalar (and thus radial) Pöschl-Teller twoparticle potentials; its details are discussed in detail in the Appendix A. In the spirit of our motivation to study the 3 He/ 3 H systems, two of the particles are assumed to be identical bosons (they could also be fermions with opposite spins), and the third is distinguishable.We assume the potential between identical particles is too weak for a bound state (using λ 1 = 1), but that the ones between any of the two identical and the third particle is strong enough for a bound state (like in the 3 S 1 2N system, using λ 2 = λ 3 = 2.2).Thus where k is "the missing index in ij", e.g., k = 3 for ij = 12.We use p = 2 in the following.
We assume the third particle is charged (or the two identical particles are and the third one is neutral, which is equivalent), and look at excitation due to a dipole operator, relative to the CoM.This gives rise to an L = 1 state in the η 3 coordinate, leaving the other coordinate unaffected.
Further details can be found in the Appendix.

A. Numerical calculation
As we can see in Fig. 10, we get a good description of the integrated response for energies ω < 0, but the results spread out directly above the three-particle threshold at ω = 0.The scatter of the individual calculations is much larger than in the absence of final state interactions, as shown in Fig. 9.So what is the origin of this?First of all, it is not the effect of channel coupling.While it does play a role, the integrated response distribution shows in the worst case scenario a singularity in the second derivative w.r.t.energy, and thus is barely visible in a plot.The lowest points correspond to cases where some response is lacking directly below zero or directly above zero, i.e., in the pure two-body channel, (typically because there is no state near zero in that specific calculation).
It is not immediate clear from the figure that each individual calculation has 4-to-5 states below zero energy.Actually, a lot of the problematic behaviour is due to states that carry no or very little dipole response.That is not a surprise: the response at low energy must be driven by the two-body channel (since the phase space of the three-body channel opens very slowly), but there are a large number of different configurations of (almost) zero response states.
That leads to the question whether this feature is peculiar to the SVM method, or whether more structured computational schemes can resolve this issue.To that end, we also tackle this model using a hyperspherical harmonics approach in the formalism set out in Ref. [27] using a simple harmonic oscillator basis; see Appendix B for details.For the angular-momentum-zero ground state, we need equal angular momentum on both Jacobi coordinates.For the J = 1 state caused by the dipole excitation, which only acts on one of the two Jacobi coordinates, we use angular momenta 0 ⊗ 1 and 2 ⊗ 1, 2 ⊗ 3, 4 ⊗ 3, etc.We use all L's that match the choices made for J = 0, with the odd angular momentum one larger than the largest one in the J = 0 state.The values of L used link to the hyperspherical quantum number K = 2κ+L 1 +L 2 .Thus we see that, with κ = 0, the cutoff in K is linked to the maximal values of L; for example, J = 0 implies L max = K max /2.One encounters some difficulties with the hyperspherical approach: while a small ω is needed for closely spaced states and a description of the 2+1 channel, it makes it also highly non-trivial to describe the ground state in sufficient detail; for the calculation reported here ( ω = 0.04), we use 90 basis states for each value of K, and use K up to 16, and thus L up to 8 for J = 0 and 9 for J = 1).The bound state energy is found as E = −0.648,which is slightly larger than the value using the SVM (E = −0.652).Since the method provides a strict upper bound, the SVM result therefore is objectively better.Obtaining any sensible values for the two-body bound state as a continuum threshold is even more difficult.Its best estimate from the SVM is E = −0.172,with a number of additional states below zero energy, whereas for the HH method we only find two, with the lowest at E = −0.108-andeven that requires substantial work.Again, the SVM result provides a better bound to the true value.As we can see in Fig. 11, the strength calculated in the HH calculation is a few percent larger than that in the SVM ones.We have been able to trace this to the slightly poorer choice of 3-body bound state wave function, as discussed above.Overall we find that there is a great similarity between the two sets of calculations.It is more difficult to find the threshold in the HH calculation-but there exist dedicated methods to improve this [3].That also links to the lack of spectral density below ω = 0. Interestingly enough, the HH calculation shows the same horizontal jumps as the SVM calculations, but more regularly spaced.Specifically, the HH flat region seen near ω = 0 agrees with the average behaviour of the SVM calculations.We thus conclude that for this type of calculation, an SVM method with a suitable basis choice methodology can be highly successful and efficient and makes it probably easier to extract a smooth result-but at the price of the need to do a number of calculations due to the statistical nature of the process.

VII. CONCLUSIONS
In this paper, we have explored a novel approach to extracting response functions from a bound-state method that uses a basis of square-integrable functions.This provides an alternative to the LIT, and may have advantages over that method in extracting responses of systems of three or more particles.
Our approach uses the SVM to solve the inhomogenous Schrödinger equation that gives the response of the wave function to an external field.The stochastically-chosen basis of this method has advantages over the commonlyused hyperspherical harmonics for describing channels with cluster structures, such as appear in break-up of systems of three of more particles, and can describe both thresholds and reponse reasonably well.
Any method using a square-integrable basis generates a response function in the form of a set of δ distributions.The LIT folds this with a Lorentzian to obtain a continuous function of energy.With care, this can be inverted to yield a continuum response function.In contrast, we consider integrated the response function.This is a stepwise continuous function of energy, with steps at random positions depending on the particular SVM basis.By running a sufficient number of independent SVM calculations, we are able to make a robust fit of a smooth function to the ensemble of integrated response functions.The derivative of this function with respect to energy then provides our approximation to the physical response function.
We have successfully tested our method using a simple model of two particles interacting via a Pöschl-Teller potential, where we find excellent agreement with analytic results.We have also applied it to a more realistic two-body system, namely deuteron break-up using the Argonne V18 potential.The agreement with the data is very good, and similar to that of other approaches that consider only the one-body E1 operator.We have also used a dispersive method to calculate the real parts of the dynamical polarisabilities corresponding to the deuteron response functions, finding good agreement with older calculations using the LIT.
Our ultimate goal is to apply the method to photodisintegration of and Compton scattering from light nuclei including 3 H, 3 He and 4 He.As a preparation for this, we have applied it to a simple three-particle model, with an interaction that generates a two-particle bound state.Like the realistic three-body systems, this has twoand three-body continuum channels.Below the threebody threshold, SVM bases give a good description of the strength in the two-body break-up channel, in contrast to HH ones.Above this threshold, the spread of the different bases is wider and so a larger ensemble of them is required for a good fit of a continuous function to the integrated response.The region immediately above this threshold has proved to be a challenging one for both SVM and HH bases, and there are indications that the response there is not well described.We plan to examine this further in future.
One aspect of physics that we have not yet explored in this approach is the treatment of resonant structure in the continuum.This will be needed for applications to 4 He and heavier nuclei.We hope to address this issue in future work, before applying the method to response functions of light nuclei.We can now evaluate the matrix element as in Fig. 12. Below, the square brackets denote a square 9J symbol, also called a recoupling symbol that arises in the recoupling of 4 angular momenta, which results in the expression in Eq. (B18).The overlap and norm are simpler, see Fig. 13, and we find, expressing cos θ 2 as 4π/3Y 1 0 (Ω 2 ) and (where L 1 = 1 in the figure)

FIG. 3 .
FIG. 3. The cumulative response for various values of the the harmonic oscillator length parameter.In a) and b) the points plotted are the midpoints of each jump (see main text for details), the only difference between these two figures is the scale.The number of basis functions used is b = 1/2: 60, b = 1: 50, b = 2: 80, b = 4: 120 and b = 8: 240 In each case, the solid blue curve is the analytical result.In c) the lower and upper lines connect the values just before and just after the jumps, respectively.Colours correspond to a) and b) (we do not show b = 1/2).

FIG. 5 .
FIG. 5. (a) calculation of the integrated response function (dots shown at midpoint of jump) for two different sets of Gaußian basis functions (blue and red dots); solid lines are the fitted smooth response distributions (b) response function obtained by differentiation of the fit from (a)-colours used correspond to the data sets.The two response distributions are visually indistinguishable.

FIG. 9 .FIG. 10 .
FIG.9.Dots: calculation of the integrated response function for a pure three-body continuum for 20 sets of independently chosen Gaußian basis functions; the solid yellow line is the analytic result Eq. (48).
of the Columbian College of Arts and Sciences; and by an Enhanced Faculty Travel Award of the Columbian Col-lege of Arts and Sciences.His research was conducted in part in GW's Campus in the Closet.

FIG. 12 .FIG. 13 .
FIG.12.Reduced matrix element for a general state.A red star denotes a time-reversed state.See[28] for an explanation of the notation used.