Homework 4; due Friday, October 18PY 502, Computational Physics (Fall 2013)Department of Physics, Boston UniversityInstructor: Anders SandvikThe deuteron, the nucleus of which consists of a neutron and a proton, constitutes the simplestcomposite nuclear system. In this assignment you are to solve the radial Schr¨odinger equationfor the bound state of the neutron-proton system, using a Yukawa potential supplemented by ahard-core short-range repulsion. You will adjust the width and depth parameters of this potentialto fit the known binding energy and radius of the deuteron.Radial Schr¨odinger equation for the neutron-proton systemThe wave function for a two-particle system with a central potential V (r) can be written in theformΨL,Lz,n(~x) = RL,n(r)YL,Lz(φ, Θ), (1)where YL,Lz(φ, Θ) are the spherical harmonics and the radial function satisfies the equation−¯h22m1rd2dr2r +L(L + 1)¯h22mr2+ V (r)RL,n(r) = EL,nRL,n(r), (2)where m is the reduced mass,m =M1M2M1+ M2. (3)Defining the functionUL,n(r) = rRL,n(r), (4)a simpler radial equation is obtained;−¯h22md2dr2+L(L + 1)¯h22mr2+ V (r)UL,n(r) = EL,nUL,n(r). (5)This is of the same form as the one-dimensional Schr¨odinger equation, apart from the fact thatthere −∞ < x < ∞ but here r ≥ 0, and the presence of the repulsive ”centrifugal barrier” whicheffectively contributes to the potential energy.For the deuteron L = 0,1and the radial function U(r) = U0,0(r) can be written asd2U(r)dr2= αV (r)E− 1, (6)where we have definedα = E2m¯h2. (7)The deuteron has only one bound state, with energy E = −2.226 MeV. The neutron and protonmasses are almost equal; Mp= 1.6726 · 10−27kg and Mp= 1.6749 · 10−27kg. Using these values,1Actually, due to a small non-central nuclear force component, L is strictly not a conserved quantum number forthe deuteron; a small amount of L = 2 wave function is mixed with the L = 0 state. We will here neglect this.1we get α = −5.3667 ·1028m−2. In order not to have to work with the very small numerical valuescorresponding to the short inter-nuclear distances expressed in meters, we change the unit of lengthfrom m to fm (1 fm = 10−15m) in Eq. (6), leading tod2U(r)dr2= βV (r)E− 1U(r), β = −0.053667. (8)Nuclear potentialA description of nuclear systems in terms of point particles governed by static central potentialsis not completely correct, but nevertheless is important as a first approximation. Several types ofmodel potentials are used, among them the Yukawa potentialVY(r) = −V0e−r/ar/a. (9)At very short distances the potential should become strongly repulsive, which is not a feature ofthe Yukawa potential. A hard-core (infinite barrier) repulsion can be included to accomplish this.Then the full potential isV (r) = ∞, for r < r0, V (r) = VY(r) = −V0e−r/ar/a, for r ≥ r0. (10)This potential will be used here. The three parameters—the hard-core radius r0, the range a, andthe depth parameter V0—can be adjusted so that known properties of the deuteron are reproduced.Here we shall consider a simplification, fixing the hard-core radius atr0= 0.1 fm. (11)The results are in fact not very sensitive to the exact value of r0. To fix the remaining twoparameters, we will use the binding energy and the radius of the deuteron. The radius is definedin terms of the expectation value of its squarehr2i =14hΨ|r2|Ψi. (12)Here the factor 1/4 comes from the fact that for two particles of equal mass, the distance betweenthem correspond to the diameter of a circular orbit, not the radius (we can here neglect the smallmass difference between the neutron and the proton). Experimentally, one cannot measure theradius directly; different radia can be defined depending on what physical scattering process ismeasured. All of the estimates are, however, close to r = 2 fm, which we will use here.Programming tasksWrite a program that solves the radial wave function written in the form (8) with the given valueof β. For a bound state, with a potential decaying exponentially to zero at long distances, theasymptotic form of the wave function is given byU(r) ∝ e−√|β|r, (r → ∞). (13)2The second boundary condition is simple, due to the hard-core;Ψ(r0) = 0. (14)In this case, it is best to start the integration from the outside, at some longest distance rmaxfromthe center (where rmaxis an input value to be read in by the program, and you have to figureout by experimentation what a suitable value is), at which the wave function is well approximatedby the form (13). The integration is done inward, and at the last point r0the second boundarycondition (14) should be satisfied. Actually, provided that rmaxis sufficiently large, the initialcondition at this distance plays a very minor role (i.e., the resulting wave function in the regionwhere it is large depends very little on it). Instead of using (13) for the two starting values U(rmax)and U(rmax− ∆r), it is therefore also fine to choose two arbitrary (preferably small,  1) valueswith U (rmax) < U(rmax− ∆r).Normally, when solving the Schr¨odinger equation we are interested in finding the energy eigenvalues.Here we are considering the corresponding inverse problem; we know the binding energy E = −2.226MeV [given in the form of the constant β and in the ratio with the potential in Eq. (8)] and theradius r =phr2i = 2 fm. We want to find the potential that gives rise to a ground state with thisenergy (and no excited bound states).For a given value of the Yukawa range parameter a in Eq. (10), your program should extract theratio V0/E for which a bound state is obtained. The energy of the bound state is negative (relativeto the potential at r = ∞, which here is 0), and hence the ratio −V0/E = V0/|E| has to bepositive. Using values V0/|E| = 0, ∆V, 2∆V, ..., first search for two values between which U (r0)changes sign. Then use bisection to find the V0/|E| for which the boundary condition U(r0) = 0is satisfied. Knowing |E| = 2.226 MeV, you then have the potential depth parameter V0that givesthe correct binding energy for the range parameter a used.2You can then calculate the radiususing the wave function corresponding to these parameters, according to (12). Here you shouldkeep in mind Eqs. (1) and (4) and note that the angular part Y of the wave function is normalizedand does not enter explicitly in an expectation value of an operator not involving the angles.Write the program in such a way that a number