c The program now prints out everything in rfs.out to 14 digits. c This program has commented out the print statement put in to print out the values of c the modes for all values of k and for all values of t at which they c are computed. c This version has g -> g^2, i.e. the interaction has been changed c from g phi^2 psi^2 to g^2 phi^2 psi^2 c c To be careful I have everywhere changed g to gnew c !$$$c This program computes the backreaction of free scalar fields !$$$c on an inflaton field in flat space in the large N limit which !$$$c is equivalent to the one loop backreaction with just one !$$$c free quantized field present. Quantum fluctuations of the !$$$c inflaton field are ignored in this approximation. !$$$c !$$$c This version is for a flat space background. It is easily !$$$c generalized to other expansions. A version to compute the !$$$c backreaction on the spacetime geometry will come later. !$$$c !$$$c The self-coupling constant for the inflaton is lamphi and the !$$$c coupling between it and the quantum field is g^2. If !$$$c g is set to zero then the phi field obeys a classical !$$$c wave equation. Also is zero in this case if the field !$$$c is in the vacuum state which is what the program is assuming.. !$$$c !$$$ The basic structure of the program is as follows: !$$$ 1. Initial values and values for various parameters are read in. !$$$ 2. A subroutine is called to compute some constants relating to !$$$ intervals used in the numerical integrations of the modes. !$$$ Simpson's rule is being used. !$$$ 3. The initial values for the modes of the quantized fields are !$$$ computed using a second order adiabatic approximation. This !$$$ needs to be changed to fourth order when the backreaction on !$$$ the spacetime geometry is included. !$$$ 4. The differential equation solver is then called. It integrates !$$$ the coupled equations for the inflaton and mode functions. It !$$$ has also been set up to integrate the equations for the scale !$$$ factor. Right now this is set up to be trivial so that the space !$$$ is flat space. However I wrote it this way so that it would be !$$$ easier to generalize the program to include the backreaction on !$$$ the geometry. IMPLICIT None Integer nn,nm,n,i,ip,nk1,nk2,nk3,nt,ntot,nbad,it,nok,ntest,junk Integer iprint,iprt,nprt,istate,iback,ifile,iw Integer ii,iii,nmode Real*8 y,dydt,c1,c2,c3,delk1,delk2,delk3,k,ti,t,delt,eps,h1,hmin Real*8 mu,mphi,xi,phi,phip,phipp,a,ap,W,Wp,lamphi,gnew,gphi2 Real*8 pi,lamk,psi2,psi2a,rhophi,rhopsi,rhopsia,rhotot Real*8 kadiabatic,omega,wm1,wterm,testphipp Real*8 beta1,beta2,beta3,Lam,psi2cutoff,psi2exact,psi2lamdif,psi2modedif Real*8 psi2dif,psi21,psi22,psi23,dif1,dif2,dif3 Real*8 epsilon,psi21exact,psi22exact,psi23exact,T1,T2,T3,T4,T5,psi2int PARAMETER(NN=100000,nm=24000) DIMENSION Y(NN),DYDT(NN),c1(nm),c2(nm),c3(nm),k(3*nm) common/param/pi,mu,mphi,xi,lamphi,gnew Common/in1/lamk,c1,c2,c3,k Common/in2/nk1,nk2,nk3 Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int Common/back/iback Common/fil/ifile EXTERNAL DERIVS EXTERNAL BSSTEP !$$$ !$$$rfs.dat is the file for all the values of parameters to be read into the !$$$program. rfs.tst is a file that includes diagnostic output to determine !$$$how accurately the mode integrations are being done. rfs.out is the main !$$$output file. !$$$ OPEN(15,FILE='rfs.dat') OPEN(16,FILE='rfsint.tst') OPEN(17,FILE='rfs.out') c OPEN(18,FILE='rfsmodesum.tst') OPEN(19,FILE='rfstmn.out') c This file is for the modes which will be followed in time. c Open(20,file='modes.out') !$$$reading the variable junk is just a trick to allow for the labeling of the !$$$data file. !$$$ !$$$mphi is the mass of the inflaton, xi is the coupling of the inflaton to !$$$the scalar curvature, lamphi is the self-coupling of the inflaton, gnew !$$$is the coupling between the inflaton and the quantum fields. !$$$ !$$$The program integrates the mode equation in a two step process. First it !$$$goes from k = 0 to a cutoff lamk. In this stage only the terms in the !$$$WKBdiv approximation that are of order k and k^3 are subtracted off. In !$$$the second stage the integral goes from lamk to an upper limit cutoff. !$$$Here the terms in WKBdiv that go like 1/k are also subtracted off. This !$$$is consistent with the way in which we define the analytic approximation. !$$$Here what is happening is the WKBdiv terms are being subtracted off and !$$$then added back on as part of the analytic approximation. So nk1 is the !$$$number of modes in the first integration and nk2 the number in the second !$$$integration. However there is an overlap at k = lamk since that is the !$$$end and beginning of the two integration ranges. delk1 and delk2 are the !$$$intervals for the two integration regions. !$$$ !$$$phi is the value of inflaton field, phip the value of its first derivative. !$$$The quantity "a" is the scale factor. In this version of course it makes !$$$sense to just set a = 1. !$$$ !$$$eps1 is the target accuracy for the differential equation solver for the !$$$quantities a, ap, phi, phip. These can be less accurate usually than the !$$$modes. eps2 is the target accuracy for the modes. h1 is an initial value !$$$for the differential equation solver to try and hmin is the minimum step !$$$size. The program stops if the step size has to go below this value. c iprt is the spacing between times when you print out the modes spectrum. READ(15,*) junk READ(15,*) mphi,xi,lamphi,gnew READ(15,*) junk READ(15,*) DELk1,nk1,delk2,nk2,delk3,nk3,DELT,NT READ(15,*) junk c READ(15,*) phi,phip,a READ(15,*) gphi2,phip,a phi = dsqrt(gphi2/gnew**2) READ(15,*) junk c READ(15,*) EPS1,EPS2,H1,HMIN READ(15,*) EPS,H1,HMIN Read(15,*) junk c istate is 2 for second order and 4 for 4th order. c iback is 0 for no backreaction on the phi field and 1 for backreaction c Note that even when there is no backreaction I will use the same 4th order adiabiabatic c state as if there was so that the comparison can be direct. Read(15,*) nprt,istate,iback c This is the n for one of the modes to be printed out as a function of time Read(15,*) junk Read(15,*) nmode Write(17,*) ' 0 gnew = ',gnew,' (g*phi)^2 = ',gphi2 If(iback.eq.0) then If(istate.eq.2) then Write(17,*) '0 Background field with adiabatic order 2 states' Else Write(17,*) '0 Background field with adiabatic order 4 states' Endif Else If(istate.eq.2) then Write(17,*) ' 0 Backreaction with adiabatic order 2 states' Else Write(17,*) ' 0 Backreaction with adiabatic order 4 states' Endif Endif pi = 2.d0*dasin(1.d0) Print *,'pi = ',pi !$$$This is the lower limit cutoff for the integrals that go like 1/k in !$$$the analytic approximation. lamk = nk1*delk1 Print *,'lower limit cutoff = ',lamk !$$$c We'll set the parameter mu in the integrals of the analytic approximation !$$$c to be equal to 1. mu = 1.d0 call Intconst(delk1,delk2,delk3) !$$$This is just for debugging purposes. Print *,'c1(1) = ',c1(1),'c1(2) = ',c1(2),' c1(nk1) = ',c1(nk1) Print *,'c2(1) = ',c2(1),'c2(2) = ',c2(2),' c2(nk2) = ',c2(nk2) Print *,'c3(1) = ',c3(1),'c3(2) = ',c3(2),' c3(nk3) = ',c3(nk3) Print *,'k(1) = ',k(1),' k(nk1) = ',k(nk1),' k(nk1+1) = ',k(nk1+1) Print *,'k(nk1+nk2) = ',k(nk1+nk2),'k(nk1+nk2+1) = ',k(nk1+nk2+1),' k(nk1+nk2+nk3) = ',k(nk1+nk2+nk3) !$$$ntest is some parameter in the differential equation solver. This solver !$$$comes from Numerical Recipes. ntest = 0 !$$$c There are 4 variables that we start with. Y(1) = a, Y(2) = ap, !$$$c Y(3) = phi, Y(4) = phip. !$$$c In the sums over modes there is an extra value of n that one must !$$$c be careful of because there is an integral from k = 0 to lamk and !$$$c then one from lamk to the upper limit cutoff for k. Now the value !$$$c of k at the lower limit cutoff is the same for both integrals but !$$$c the integration factors c1 and c2 are different. !$$$c !$$$c There are 4 modes for each k; the real part, the imaginary part !$$$c and their derivatives. !$$$c !$$$c Thus the total number of variables for the differential equation !$$$c solver is 4 + 4*nk1 + 4*(nk2 - 1)+4*(nk3-1) again because of the redundancy !$$$c in the value of k. Thus !$$$c NTOT = 4*nk1 + 4*nk2 + 4*nk3 Y(1) = a !$$$Let's work in a flat space background. ap = 0.d0 Y(2) = ap Y(3) = phi Y(4) = phip !$$$c !$$$c Note that Y(5) = Real psi, Y(6) = Real psip, Y(7) = Im psi, !$$$c Y(8) = Im psip. !$$$ We have the same value of k at n = nk1 and n = nk1+1. !$$$c This is because we want different values for c1(nk1) !$$$c and c2(nk1+1) even though the value of k is the same !$$$c for these since they are at the end of the first integration !$$$c interval and the beginning of the second respectively. c We can have a problem with initial values for the modes if g^2 is too small as c it can be for small g^2 phi_0^2. This checks for that and then sets iw = 1 so that c later on the second order WKB approximation can be used for small values of k. c c Now we keep the value of phipp the same for simplicity. It is too complicated to do anything else. c This is not unreasonable so long as the correction to phipp is small and that is checked c farther below. Here are some reasons why c i) you have a UV cutoff anyway. c ii) you have numerical error anyway c c c This is for a fourth order adiabatic state and the value of phipp is determined c self-consistently unless the starting value for Wm1 < 0 for one of the k's used. c c I do it this way to eliminate infrared divergences. c If(istate.eq.4) then phipp = 1.d0/(1.d0+gnew**2/(48.d0*pi**2))*(-phi + gnew**4*phi**3/(16.d0*pi**2) - - gnew**4*phi**3*dlog(gnew**2*phi**2/mu**2)/(16.*pi**2)) c Now test to see if Wm1 is going to be positive. If it is for the smallest k used then there is c no problem. Omega = dsqrt(k(1)**2 + gnew**2*a**2*phi**2) Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.*Omega**5) If(Wm1.lt.0) then iw = 1 Else iw = 0 Endif Print*,'If iw = 0 then 4th order adiabatic was used everywhere which is what we want. iw = ',iw c Now see if at nk1 it is OK. Omega = dsqrt(k(nk1)**2 + gnew**2*a**2*phi**2) Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**5) If(Wm1.le.0.d0) then Print *,' Wm1 < 0 even at k(nk1)' Stop Endif Else phipp = -phi Endif !$$$ The Do loop sets the intial values for the modes. c If istate = 2 it is 2nd order adiabatic and if istate = 4 it is fourth order c In both cases it is an infrared finite state. n = 0 DO i=5,NTOT-3,4 If(n.ne.nk1) then n = n+1 else n = n+2 endif Omega = dsqrt(k(n)**2 + gnew**2*a**2*phi**2) If(istate.eq.2) then W = Omega Wm1 = 1.d0/Omega else c Note to get the normalization correct I will use Wm1 everywhere so W = 1/Wm1 c This is useful because I used Wm1 to fourth adiabtic order to compute phipp c Also as you can see in the commented statement, if I compute W directly I will c sometimes get a negative number because of the minus sign. There is no such c problem with Wm1 at large k but there can be at small k. c W = Omega - gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**3) Wm1 = 1.d0/Omega + gnew**2*(phip**2+phi*phipp)/(4.d0*Omega**5) c Use second order adiabatic up to and including k(nk1). There will be a sharp transition but c there already is one because of the difference in subtraction terms on either side this value of k c And of course phipp will be wrong but hopefully still close to the correct value If (iw.eq.1.and.n.le.nk1) then Wm1 = 1.d0/Omega Endif W = 1.d0/Wm1 Endif c This value of Wp is 3rd order but in flat space with phip=0 as a starting value it is zero c anyway. Wp = gnew**2/(Omega)*(a**2*phi*phip + a*ap*phi**2) Y(i) = dsqrt(Wm1/2.D0) Y(i+1) = -Wp*(Wm1/2.D0)**1.5D0 Y(i+2) = 0.D0 Y(i+3) = -(W/2.D0)**.5D0 c I want to write the initial values of the modes as well c t = 0 so it is easiest just to set it to that here t = 0.d0 c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3) c We need to print this out twice if n = nk1 because the do loop skips n = nk1+1 but it is easier c to read that value in the correlation function routine. c If(n.eq.nk1) then c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3) c Endif end do !$$$c This is to print out the initial values of app and phipp as well !$$$c as the test integral and psi2 and psi2a. The latter are printed !$$$c out through a command in the test subroutine. TI = 0.d0 T = TI !$$$DERIVS is a subroutine that steps the equations forward in time. !$$$It is written by me and is used for the ODE solver. Here it is !$$$being used to get initial values for certain quantities. !$$$TEST is a subroutine that integrates a simple function that has the same !$$$type of k behavior as does the actual mode integral for the quantity !$$$<\psi>^2. It is therefore useful for determining whether the mode integration !$$$is likely to be accurate or not, and by how much it might deviate !$$$in its accuracty. CALL DERIVS(T,Y,DYDT,NTOT) WRITE(16,*) ' Time dif int 1 dif int 2 dif int 3 diff total' CALL TEST(T,Y,ntot) c This prints out the integrand for the mode sum. This is useful c if gnew = 0 since then the integrand should always vanish. ifile = 0 CALL TestModeSum(T,Y,DYDT,NTOT) rhotot = rhophi+rhopsi c This is just for debugging purposes CALL DERIVS(T,Y,DYDT,NTOT) c c We can also compare the initial value of with its correct value. There is an c error from both the cutoff and the density of modes. I'll compute them separately. c Here psi2cutoff is the exact value we should get with with cutoff. psi2exact is the exact c value we should get if the cutoff was infinity. Lam = k(nk1+nk2+nk3) c This is for the first integration interval which has no 1/k subtraction term. epsilon = k(nk1) psi21exact = (1.d0/(4.d0))*(epsilon*dsqrt(epsilon**2+gphi2)-epsilon**2 1 +gphi2*(-dlog((epsilon+dsqrt(epsilon**2+gphi2))/dsqrt(gphi2))) 2 + (phipp/(6.d0*phi))*epsilon**3/(epsilon**2+gphi2)**1.5d0) dif1 = (psi21exact-psi21)/psi21exact c This is for the second integration interval T1 = 0.5d0*(k(nk1+nk2)*dsqrt(k(nk1+nk2)**2+gphi2) 1 - gphi2*dlog((k(nk1+nk2)+dsqrt(k(nk1+nk2)**2+gphi2)))) T2 = -0.5d0*(k(nk1)*dsqrt(k(nk1)**2+gphi2) - gphi2*dlog((k(nk1)+dsqrt(k(nk1)**2+gphi2)))) T3 = -0.5d0*(k(nk1+nk2)**2-k(nk1)**2) T4 = dlog(k(nk1+nk2)/k(nk1)) T5 = (1.d0/(3.d0*gphi2))*(k(nk1+nk2)**3/(k(nk1+nk2)**2+gphi2)**1.5d0 1 - k(nk1)**3/(k(nk1)**2+gphi2)**1.5d0) psi22exact = (0.5d0*T1+0.5d0*T2+0.5d0*T3+0.25d0*gphi2*T4+(gphi2*phipp/(8.d0*phi))*T5) dif2 = (psi22exact-psi22)/psi22exact T1 = 0.5d0*(Lam*dsqrt(Lam**2+gphi2) 1 - gphi2*dlog((Lam+dsqrt(Lam**2+gphi2)))) T2 = -0.5d0*(k(nk1+nk2)*dsqrt(k(nk1+nk2)**2+gphi2) - gphi2*dlog((k(nk1+nk2)+dsqrt(k(nk1+nk2)**2+gphi2)))) T3 = -0.5d0*(Lam**2-k(nk1+nk2)**2) T4 = dlog(Lam/k(nk1+nk2)) T5 = (1.d0/(3.d0*gphi2))*(Lam**3/(Lam**2+gphi2)**1.5d0 1 - k(nk1+nk2)**3/(k(nk1+nk2)**2+gphi2)**1.5d0) psi23exact = (0.5d0*T1+0.5d0*T2+0.5d0*T3+0.25d0*gphi2*T4+(gphi2*phipp/(8.d0*phi))*T5) dif3 = (psi23exact-psi23)/psi23exact Print *,'psi21 = ',psi21,' psi21exact = ',psi21exact,' dif1 = ',dif1 Print *,'psi22 = ',psi22,' psi22exact = ',psi22exact,' dif2 = ',dif2 Print *,'psi23 = ',psi23,' psi23exact = ',psi23exact,' dif3 = ',dif3 psi2cutoff = (1.d0/(8.d0*Pi**2))*(Lam*dsqrt(Lam**2+gphi2)-Lam**2+gphi2*(-dlog((Lam+dsqrt(Lam**2+gphi2))/Lam) 1 +dlog(dsqrt(gphi2)*2.d0/mu)) + (phipp/(6.d0*phi))*Lam**3/(Lam**2+gphi2)**1.5d0 - gphi2 ) Print*,' Original value of psi2cutoff = ',psi2cutoff If(iback.eq.1) then psi2cutoff = psi2cutoff + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0 1 + gnew**2*phip**2/8.d0 + gnew**2*phi*phipp/8.d0) Print*,' Value of psi2cutoff when 4th order WKB approximation is added on = ',psi2cutoff Endif psi2exact = phipp/(48.d0*pi**2*phi) + (gphi2/(16.d0*pi**2))*(-1.d0 + dlog(gphi2/mu**2)) psi2dif = (psi2exact-psi2)/psi2exact psi2lamdif = (psi2exact-psi2cutoff)/psi2exact psi2modedif = (psi2cutoff-psi2)/psi2cutoff Print*,'psi2 = ',psi2,'psi2a = ',psi2a Print*,'psi2exact = ',psi2exact,' psi2cutoff = ',psi2cutoff Print*,' Total relative difference =',psi2dif Print*,' Cutoff = ',Lam Print*, ' The relative error due to the existence of a cutoff = ', psi2lamdif Print*,' The relative error in the mode integral = ',psi2modedif Write(17,*),'0 Cutoff = ',Lam Write(17,*) ' 0 Relative error due to the cutoff = ', psi2lamdif Write(17,*),'0 Relative error for mode integral = ',psi2modedif c We can check the accuracy of the mode sum by comparing the initial value of phipp c above which is computed exactly analytically with DYDT(4) which is the numerical c calculation of the initial value of phipp. c This only makes sense to do for a 4th order adiabatic state when backreaction is being c taken into account If(istate.eq.4.and.iback.eq.1) then testphipp = (phipp-DYDT(4))/phipp Print *, 'Initial phipp analytic = ',phipp,' Initial phipp numeric = ',DYDT(4), - ' percent difference = ', testphipp Print *, ' If iw = 1 the difference needs to be small since phipp is wrong. iw = ',iw Write(17,*) '0 Initial phipp analytic = ',phipp Write(17,*) '0 Initial phipp numeric = ',DYDT(4) Write(17,*) '0 Percent difference = ', testphipp Endif WRITE(17,*) ' 0 Time phi phip phipp psi2 psi2a' WRITE(19,*) ' Time rhophi rhopsi rhopsia rho_total' c beta2 = |beta|^2 = Carmen and Emil's definition of produced particles WRITE(17,27) T,Y(3),Y(4),DYDT(4),psi2,psi2a WRITE(19,27) T,rhophi,rhopsi,rhopsia,rhotot c !$$$ !$$$Now start the do loop that integrates the equations forward in time. !$$$Ti is the starting time interval for the solver and t is the ending !$$$interval. Delt is the spacing in time between intervals. nt is the !$$$number of intervals. ip = 0 iprint = 1 ifile = 1 iprt = NT/nprt DO IT = 1,NT ip = ip+1 print *,ip TI = T T = T+DELT c CALL ODEINT(Y,NTOT,TI,T,EPS1,EPS2,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) CALL ODEINT(Y,NTOT,TI,T,EPS,H1,HMIN,NOK,NBAD,DERIVS,BSSTEP) CALL DERIVS(T,Y,DYDT,NTOT) CALL TEST(T,Y,ntot) c This prints out the integrand for the mode sum. This is useful c if gnew = 0 since then the integrand should always vanish. If(iprint.eq.iprt) then CALL TestModeSum(T,Y,DYDT,NTOT) iprint = 1 else iprint = iprint+1 endif rhotot = rhophi+rhopsi WRITE(17,27) T,Y(3),Y(4),DYDT(4),psi2,psi2a WRITE(19,27) T,rhophi,rhopsi,rhopsia,rhotot c beta1 = (1.d0/(2.D0*k(1)))*(Y(6)**2+Y(8)**2)+(-Y(6)*Y(7)+Y(8)*Y(5)) c - +(k(1)/2.d0)*(Y(1)**2+Y(3)**2) c If(nmode.le.nk1) then c ii = 4*(nmode)+1 c Else c ii = 4*(nmode-1)+1 c Endif c beta2 = (1.d0/(2.D0*k(nmode)))*(Y(ii+1)**2+Y(ii+3)**2)+(-Y(ii+1)*Y(ii+2)+Y(ii+3)*Y(ii)) c - +(k(nmode)/2.d0)*(Y(ii)**2+Y(ii+2)**2) c iii = 4*(nk1+nk2+nk3-1)+1 c beta3 = (1.d0/(2.D0*k(nk1+nk2+nk3)))*(Y(iii+1)**2+Y(iii+3)**2)+(-Y(iii+1)*Y(iii+2)+Y(iii+3)*Y(iii)) c - +(k(nk1+nk2+nk3)/2.d0)*(Y(iii)**2+Y(iii+2)**2) c Now print out the real and imaginary parts of the modes n = 0 DO i=5,NTOT-3,4 If(n.ne.nk1) then n = n+1 else n = n+2 endif c Note that this way I have all the modes at each value of t c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3) c I need to print it out twice if n = nk1. c If(n.eq.nk1) then c Write(20,*) t,k(n),Y(i),Y(i+1),Y(i+2),Y(i+3) c Endif Enddo If(iback.eq.0.and.dabs(psi2).gt.1.d20) then Print *,'phi2 has reached 1.d20' Stop Endif 27 format(10(e21.14,1x)) end do STOP END Subroutine Intconst(delk1,delk2,delk3) !$$$c This subroutine computes the constants needed for the mode integrations. !$$$It uses Simpson's rule. !$$$c The integrand is assumed to be equal to zero for k = 0. !$$$c If the number of intervals is wrong use the trapezoidal rule. Implicit None Integer n, nm, nk1,nk2,nk3, j Real*8 c1,c2,c3,c12,c22,c32,k,kk,delk1,delk2,delk3,lamk Parameter (nm = 24000) Dimension c1(nm),c2(nm),c3(nm),k(3*nm),c12(nm),c22(nm),c32(nm) Common/in1/lamk,c1,c2,c3,k Common/in2/nk1,nk2,nk3 Common/in3/c12,c22,c32 !$$$This subroutine computes the intervals of integration. !$$$There are two intervals because we have a lower limit cutoff and so want !$$$to integrate a certain set up terms up to that and then all other terms !$$$starting at the cutoff. !$$$Note that for the last interval we want to use the trapezoidal rule !$$$if the number of intervals does not work out right. This changes !$$$the value of c(n)for the last two values of n. The same thing happens !$$$at the end of the second interval of integration. !$$$ !$$$ c c I have now split the second interval into two parts to make the integration more accurate. c kk = 0.d0 j = 0 do n=1,nk1 kk = kk+delk1 k(n) = kk if(n.eq.nk1) then If(j.eq.0) then c1(n) = delk1/2.d0 else c1(n) = delk1/3.d0 endif else if(n.eq.nk1-1) then if(j.eq.1) then c1(n) = delk1*(1.d0/3.d0 + 0.5d0) j = 0 else c1(n) = 4.d0*delk1/3.d0 j = 1 endif else if(j.eq.0) then c1(n) = 4.d0*delk1/3.d0 j = 1 else c1(n) = 2.d0*delk1/3.d0 j = 0 endif enddo !$$$ !$$$Here we start at the lower limit cutoff so the sequence is different. !$$$Also k(nk1+1) = k(nk1) because of the way the limits of integration are. j = 0 kk = kk - delk2 do n=1,nk2 kk = kk + delk2 k(n+nk1) = kk if(n.eq.nk2) then If(j.eq.0) then c2(n) = delk2/2.d0 else c2(n) = delk2/3.d0 endif else if(n.eq.nk2-1) then if(j.eq.1) then c2(n) = delk2*(1.d0/3.d0 + 0.5d0) j = 0 else c2(n) = 4.d0*delk2/3.d0 j = 1 endif else if (n.eq.1) then c2(n) = delk2/3.d0 else if(j.eq.0) then c2(n) = 4.d0*delk2/3.d0 j = 1 else c2(n) = 2.d0*delk2/3.d0 j = 0 endif enddo !$$$ !$$$Here we start at the lower limit cutoff so the sequence is different. !$$$Also k(nk1+nk2+1) = k(nk1+nk2) because of the way the limits of integration are. j = 0 kk = kk - delk3 do n=1,nk3 kk = kk + delk3 k(n+nk1+nk2) = kk if(n.eq.nk3) then If(j.eq.0) then c3(n) = delk3/2.d0 else c3(n) = delk3/3.d0 endif else if(n.eq.nk3-1) then if(j.eq.1) then c3(n) = delk3*(1.d0/3.d0 + 0.5d0) j = 0 else c3(n) = 4.d0*delk3/3.d0 j = 1 endif else if (n.eq.1) then c3(n) = delk3/3.d0 else if(j.eq.0) then c3(n) = 4.d0*delk3/3.d0 j = 1 else c3(n) = 2.d0*delk3/3.d0 j = 0 endif enddo c Now do the intervals when I have 1/2 the number of points c I am going to simply assume there is an odd number of points c and that the intervals work out correctly c You need (nk-1) to be divisible by 4 for the second two bands. c For the first band you need nk1 to be divisible by 4 since you c start with n = 0 (i.e. k = 0) but that doesn't contribute. c Throughout delk is multiplied by 2. j = 0 c We start with n = 2 because of the fact that we are skipping every other step and n = 0 gives c no contribution because k = 0 do n=2,nk1,2 if(n.eq.nk1) then c12(n) = 2.d0*delk1/3.d0 else if(j.eq.0) then c12(n) = 4.d0*2.d0*delk1/3.d0 j = 1 else c12(n) = 2.d0*2.d0*delk1/3.d0 j = 0 endif enddo !$$$ !$$$Here we start at the lower limit cutoff so the sequence is different. !$$$Also k(nk1+1) = k(nk1) because of the way the limits of integration are. c Here we start with n = 1 because the beginning of the interval does contribute j = 0 do n=1,nk2,2 if(n.eq.nk2.or.n.eq.1) then c22(n) = 2.d0*delk2/3.d0 else if(j.eq.0) then c22(n) = 4.d0*2.d0*delk2/3.d0 j = 1 else c22(n) = 2.d0*2.d0*delk2/3.d0 j = 0 endif enddo !$$$ !$$$Here we start at the lower limit cutoff so the sequence is different. !$$$Also k(nk1+nk2+1) = k(nk1+nk2) because of the way the limits of integration are. j = 0 do n=1,nk3,2 if(n.eq.nk3.or.n.eq.1) then c32(n) = 2.d0*delk3/3.d0 else if(j.eq.0) then c32(n) = 4.d0*2.d0*delk3/3.d0 j = 1 else c32(n) = 2.d0*2.d0*delk3/3.d0 j = 0 endif enddo Return End !$$$This is the subroutine called by the differential equation solver to !$$$step the equations forward in time. SUBROUTINE DERIVS(T,Y,DYDT,NVAR) IMPLICIT None Integer n,nk1,nk2,nk3,nvar,nm,i,ntest,nt,iback Real*8 t,Y, DYDT, k, om2, c1, c2, c3,psi2, lamphi, gnew, f2, mphi Real*8 mu,capom2,a,ap,app,phi,phip,xi,pi,psi2a,lamk Real*8 f2p,int1,int2,int3,rhopsi,rhopsia,rhophi,psi21,psi22,psi23 Real*8 int11,psi2int Parameter (nm = 24000) DIMENSION Y(NVAR),DYDT(NVAR),c1(nm),c2(nm),c3(nm),k(3*nm) common/param/pi,mu,mphi,xi,lamphi,gnew Common/in1/lamk,c1,c2,c3,k Common/in2/nk1,nk2,nk3 Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int Common/back/iback !$$$ Y(1) is a, Y(2) is ap, Y(3) is phi, y(4) is phip. DYDT(1) is ap, !$$$ DYDT(2) is ap, DYDT(3) is phip, DYDT(4) is phipp. !$$$ Y(N) is Real(f(k)), Y(N+1) is Real(f'(k)), Y(N+2) IS IM(f(k)) !$$$ Y(N+3) IS IMP(f(k)) AND SIMILARLY FOR THE DYDT'S. !$$$ This version of the program is for flat space. That's why !$$$I explicitly set a = 1, ap = app = 0. a = 1.d0 ap = 0.d0 app = 0.d0 phi = Y(3) phip = Y(4) DYDT(1) = ap DYDT(3) = phip rhophi = (phip**2 + mphi**2*phi**2)/2.d0 !$$$ We have the same value of k at n = nk1 and n = nk1+1. !$$$c This is because we want different values for c1(nk1) !$$$c and c2(nk1+1) even though the value of k is the same !$$$c for these since they are at the end of the first integration !$$$c interval and the beginning of the second respectively. !$$$We first update the values of the derivatives of the modes. n = 0 DO i=5,NVAR-3,4 If(n.ne.nk1) then n = n+1 else n = n+2 endif capom2 = k(n)**2+gnew**2*a**2*phi**2 DYDT(i) = Y(i+1) DYDT(i+1) = -capom2*Y(i) DYDT(i+2) = Y(i+3) DYDT(i+3) = -capom2*Y(i+2) end do !$$$psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and !$$$c2(n) = const. delk2. psi2 = 0.D0 rhopsi = 0.d0 DO n = 1,nk1 i = 4*n + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) INT1 = k(n)**2*f2 - k(n)/2.d0 int2 = k(n)**2*(f2p+(k(n)**2+a**2*gnew**2*phi**2)*f2) - - k(n)**3 - k(n)*gnew**2*phi**2*a**2/2. psi2 = psi2 + int1*c1(n) rhopsi = rhopsi + int2*c1(n) Enddo psi21 = psi2 !$$$Now n is for c2(n) and recall k(nk1) = k(nk1+1) psi22 = 0.d0 DO n = 1,nk2 i = 4*(n-1+nk1) + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1)) int2 = k(n+nk1)**2*(f2p+(k(n+nk1)**2+a**2*gnew**2*phi**2)*f2) - - k(n+nk1)**3 - k(n+nk1)*gnew**2*phi**2*a**2/2.d0 + - gnew**4*phi**4*a**4/(8.d0*k(n+nk1)) psi2 = psi2 + int1*c2(n) psi22 = psi22 + int1*c2(n) rhopsi = rhopsi + int2*c2(n) Enddo psi23 = 0.d0 DO n = 1,nk3 i = 4*(n-1+nk1+nk2) + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) c int11 = k(n+nk1+nk2)**2*f2 int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2)) int2 = k(n+nk1+nk2)**2*(f2p+(k(n+nk1+nk2)**2+a**2*gnew**2*phi**2)*f2) - - k(n+nk1+nk2)**3 - k(n+nk1+nk2)*gnew**2*phi**2*a**2/2.d0 + - gnew**4*phi**4*a**4/(8.d0*k(n+nk1+nk2)) psi2 = psi2 + int1*c3(n) psi23 = psi23 + int1*c3(n) rhopsi = rhopsi + int3*c3(n) Enddo !$$$c Now put in the normalization factor. psi2 = psi2/(2*pi**2*a**2) rhopsi = rhopsi/(4*pi**2*a**4) !$$$c Now add on the analytic approximation. !$$$Note that lamk is the lower limit cutoff. We called it epsilon in our paper. !$$$c Note also that when we do backreaction of the geometry the app will !$$$c be a problem that will have to be dealt with. It will have to be !$$$c computed from the backreaction equations before psi2a is computed. psi2a = -app/(48*pi**2*a**3)-gnew**2*phi**2/(8*pi**2)*(1 - dlog(2*lamk/(a*mu))) psi2 = psi2 + psi2a c rhopsia = - gnew**4*phi**4/(32.d0*pi**2) + gnew**2*phi**4/(32.d0*pi**2) *dlog(2*lamk/(a*mu)) rhopsia = - gnew**4*phi**4/(32.d0*pi**2) + gnew**4*phi**4/(32.d0*pi**2) *dlog(2*lamk/(a*mu)) rhopsi = rhopsi + rhopsia DYDT(4) = -2.d0*ap/a*phip - a**2*(mphi**2 + xi * 6.d0*DYDT(2)/a**3 - + lamphi/6.d0*phi**2 + iback*gnew**2*psi2)*phi c Finally put in the terms from the 4th order WKB approximation for modes above the cutoff If(iback.eq.1) then DYDT(4) = (1.d0/(1.d0+gnew**4*phi**2/(32.d0*pi**2*k(nk1+nk2+nk3)**2)))*(DYDT(4) 1 - 3.d0*gnew**6*phi**5/(64.d0*pi**2*k(nk1+nk2+nk3)**2) 2 - gnew**4*phip**2*phi/(32.d0*pi**2*k(nk1+nk2+nk3)**2)) psi2 = psi2 + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0 1 + gnew**2*phip**2/8.d0 + gnew**2*phi*DYDT(4)/8.d0) Endif c c This is enough if we work in a fixed background. Then I just need c specify DYDT(2) which is a''. DYDT(2) = 0.d0 RETURN END !$$$This subroutine computes the modes sums when half the modes are used. SUBROUTINE TEST(t,Y,nvar) IMPLICIT None Integer nm,nk1,nk2,nk3,j,n,i,nvar Real*8 t,c1,c2,c3,k,w,pi,mphi,xi,lamphi,gnew,int,dif1,dif2,dif3,a Real*8 lamk,phi,psi2a,psi2,mu,rhophi,rhopsi,rhopsia Real*8 c12,c22,c32,psi21,psi22,psi23,psi2t1,psi2t2,psi2t3,psi2t Real*8 psi2int,diftot,f2,Y,int1 Parameter(nm = 24000) DIMENSION Y(NVAR) Dimension c1(nm),c2(nm),c3(nm),k(3*nm) Dimension c12(nm),c22(nm),c32(nm) Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int Common/in1/lamk,c1,c2,c3,k Common/in2/nk1,nk2,nk3 Common/in3/c12,c22,c32 common/param/pi,mu,mphi,xi,lamphi,gnew a = 1.d0 phi = Y(3) J = 0 psi2t1 = 0.d0 DO n = 2,nk1,2 i = 4*n + 1 f2 = (Y(i)**2+Y(i+2)**2) INT1 = k(n)**2*f2 - k(n)/2.d0 psi2t1 = psi2t1 + int1*c12(n) Enddo psi2t2 = 0.d0 c We are starting with n = 0 really although that gives zero contribution. c Since we skip every other mode we should start thus with n = 2. c For the other bands we start with n = 1 because the beginning of the band does contribute DO n = 1,nk2,2 i = 4*(n-1+nk1) + 1 f2 = (Y(i)**2+Y(i+2)**2) int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1)) psi2t2 = psi2t2 + int1*c22(n) Enddo psi2t3 = 0.d0 DO n = 1,nk3 i = 4*(n-1+nk1+nk2) + 1 f2 = (Y(i)**2+Y(i+2)**2) int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2)) psi2t3 = psi2t3 + int1*c32(n) Enddo psi2t = (psi2t1+psi2t2+psi2t3)/(2*pi**2*a**2) + psi2a dif1 = (psi21-psi2t1)/psi21 dif2 = (psi22-psi2t2)/psi22 dif3 = (psi23-psi2t3)/psi23 diftot = (psi2-psi2t)/psi2 Write(16,17) t,psi2t1,psi2t2,psi2t3 WRITE(16,17) t,dif1,dif2,dif3,diftot 17 format(8(e12.5,2x)) RETURN END c c This subroutine is used to check the convergence of the mode sum. c c In version 2 there was a program that was checking to see what c happens to the integrand when you should have the modes canceling c with the counterterms. The upshot seemed to be that the cancelation c gets worse at large k because the integrand goes like k^2 or something c like that. c SUBROUTINE TestModeSum(T,Y,DYDT,NVAR) IMPLICIT None Integer n,nk1,nk2,nk3,nvar,nm,i,ntest,nt,ifile,nfile,iback Integer nn,j,nnk2,nnk3 Real*8 t,Y, DYDT, k, om2, c1, c2, c3, psi2, lamphi, gnew, f2, f2p,int1,int2, mphi Real*8 mu,capom2,a,ap,app,phi,phip,phipp,xi,pi,psi2a,lamk Real*8 difk,difk1,psi2dif,psi2old,inttest,inttest1,cc,rhophi,rhopsi,rhopsia Real*8 beta2,psi21,psi22,psi23,psi2int Parameter (nm = 24000) DIMENSION Y(NVAR),DYDT(NVAR),c1(nm),c2(nm),c3(nm),k(3*nm),cc(3) common/param/pi,mu,mphi,xi,lamphi,gnew Common/in1/lamk,c1,c2,c3,k Common/in2/nk1,nk2,nk3 c Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia Common/der/psi2,psi2a,rhophi,rhopsi,rhopsia,psi21,psi22,psi23,psi2int Common/fil/ifile Common/back/iback ifile = ifile+1 nfile = ifile+30 If(ifile.eq.1) then Open(nfile,file='integrand.1') else if(ifile.eq.2) then Open(nfile,file='integrand.2') else if(ifile.eq.3) then Open(nfile,file='integrand.3') else if(ifile.eq.4) then Open(nfile,file='integrand.4') else if(ifile.eq.5) then Open(nfile,file='integrand.5') else if(ifile.eq.6) then Open(nfile,file='integrand.6') else if(ifile.eq.7) then Open(nfile,file='integrand.7') else if(ifile.eq.8) then Open(nfile,file='integrand.8') else if(ifile.eq.9) then Open(nfile,file='integrand.9') else if(ifile.eq.10) then Open(nfile,file='integrand.10') else if(ifile.eq.11) then Open(nfile,file='integrand.11') else if(ifile.eq.12) then Open(nfile,file='integrand.12') else if(ifile.eq.13) then Open(nfile,file='integrand.13') else if(ifile.eq.14) then Open(nfile,file='integrand.14') else if(ifile.eq.15) then Open(nfile,file='integrand.15') else if(ifile.eq.16) then Open(nfile,file='integrand.16') else if(ifile.eq.17) then Open(nfile,file='integrand.17') else if(ifile.eq.18) then Open(nfile,file='integrand.18') else if(ifile.eq.19) then Open(nfile,file='integrand.19') else if(ifile.eq.20) then Open(nfile,file='integrand.20') else if(ifile.eq.21) then Open(nfile,file='integrand.21') Endif !$$$ Y(1) is a, Y(2) is ap, Y(3) is phi, y(4) is phip. DYDT(1) is ap, !$$$ DYDT(2) is ap, DYDT(3) is phip, DYDT(4) is phipp. !$$$ Y(N) is Real(f(k)), Y(N+1) is Real(f'(k)), Y(N+2) IS IM(f(k)) !$$$ Y(N+3) IS IMP(f(k)) AND SIMILARLY FOR THE DYDT'S. phi = Y(3) phip = Y(4) phipp = DYDT(4) a = 1.d0 !$$$ We have the same value of k at n = nk1 and n = nk1+1. !$$$c This is because we want different values for c1(nk1) !$$$c and c2(nk1+1) even though the value of k is the same !$$$c for these since they are at the end of the first integration !$$$c interval and the beginning of the second respectively. c Write(18,*) 't = ',t c Write(18,*) ' n k(n) Real f Imag f int1 int2 psi2 psi2dif' !$$$ psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and !$$$ c2(n) = const. delk2. cc(1) = c2(1) cc(2) = c2(2) cc(3) = cc(1) psi2 = 0.D0 psi2old = psi2 DO n = 1,nk1 i = 4*n + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) INT1 = k(n)**2*f2 - k(n)/2.d0 c This was used to check that the sums c This is the % difference. c inttest = f2 - 1.d0/(2.d0*k(n)) c - + gnew**2*phi**2*a**2/(4.d0*k(n)**3) c difk = inttest/f2 c difk1 = -1.d0 c This is the integrand for the energy density of the psi field int2 = k(n)**2*(f2p+(k(n)**2+a**2*gnew**2*phi**2)*f2) - - k(n)**3 - k(n)*gnew**2*phi**2*a**2/2. psi2 = psi2 + int1*c1(n) c c This is to check the convergence of the modes. c psi2dif = (psi2-psi2old)/psi2 psi2old = psi2 c beta2 = |beta|^2 beta2 = (1.d0/(2.D0*k(n)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(1)*Y(2)+Y(3)*Y(1)) - +(k(n)/2.d0)*(Y(i)**2+Y(i+2)**2) c Write(nfile,37) t,n,k(n),Y(i),Y(i+2),int1,int2,psi2,psi2dif Write(nfile,37) t,n,k(n),Y(i),beta2,int1,int2,psi2,psi2dif Enddo c psi21 = psi2 37 Format(1x,d12.5,1x,i5,1x,7(d12.5,1x)) !$$$Now n is for c2(n) and recall k(nk1) = k(nk1+1) c This is to make the intervals work better. c c We need an odd number of n's because there are 3 n's in each interval. c We integrate in each interval from nn-1 to nn+1. nnk2 = ((nk2-1)/2)*2 c psi22 = 0.d0 DO nn = 2,nnk2,2 Do j=1,3 n = nn-2 + j i = 4*(n-1+nk1) + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) int1 = k(n+nk1)**2*f2 - k(n+nk1)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1)) c inttest = f2 - 1.d0/(2.d0*k(n+nk1)) + c - gnew**2*phi**2*a**2/(4.d0*k(n+nk1)**3) c difk = inttest/f2 c This is to test whether a different subtraction would lead to c a faster convergence of the mode sums. c inttest1 = f2 - 1.d0/(2.d0*dsqrt(k(n+nk1)**2 + c 1 gnew**2*phi**2*a**2)) c difk1 = inttest1/f2 int2 = k(n+nk1)**2*(f2p+(k(n+nk1)**2+a**2*gnew**2*phi**2)*f2) - - k(n+nk1)**3 - k(n+nk1)*gnew**2*phi**2*a**2/2. + - gnew**4*phi**4*a**4/(8.*k(n+nk1)) psi2 = psi2 + int1*cc(j) c psi22 = psi22 + int1*cc(j) End Do psi2dif = (psi2-psi2old)/psi2 psi2old = psi2 beta2 = (1.d0/(2.D0*k(n+nk1)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(1)*Y(2)+Y(3)*Y(1)) - +(k(n+nk1)/2.d0)*(Y(i)**2+Y(i+2)**2) c Write(nfile,37) t,n,k(n+nk1),Y(i),Y(i+2),int1,int2,psi2,psi2dif Write(nfile,37) t,n,k(n+nk1),Y(i),beta2,int1,int2,psi2,psi2dif Enddo !$$$Now n is for c3(n) and recall k(nk1+nk2) = k(nk1+nk2+1) !$$$ psi2 is <\psi^2>. Recall that c1(n) = const. delk1 and !$$$ c2(n) = const. delk2, c3(n) = const*delk3 cc(1) = c3(1) cc(2) = c3(2) cc(3) = c3(1) c This is to make the intervals work better. c c We need an odd number of n's because there are 3 n's in each interval. c We integrate in each interval from nn-1 to nn+1. nnk3 = ((nk3-1)/2)*2 c psi23 = 0.d0 DO nn = 2,nnk3,2 Do j=1,3 n = nn-2 + j i = 4*(n-1+nk1+nk2) + 1 f2 = (Y(i)**2+Y(i+2)**2) f2p = (Y(i+1)**2+Y(i+3)**2) int1 = k(n+nk1+nk2)**2*f2 - k(n+nk1+nk2)/2.d0 + - gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2)) c inttest = f2 - 1.d0/(2.d0*k(n+nk1+nk2)) + c - gnew**2*phi**2*a**2/(4.d0*k(n+nk1+nk2)**3) c difk = inttest/f2 c This is to test whether a different subtraction would lead to c a faster convergence of the mode sums. c inttest1 = f2 - 1.d0/(2.d0*dsqrt(k(n+nk1+nk2)**2 + c 1 gnew**2*phi**2*a**2)) c difk1 = inttest1/f2 int2 = k(n+nk1+nk2)**2*(f2p+(k(n+nk1+nk2)**2+a**2*gnew**2*phi**2)*f2) - - k(n+nk1+nk2)**3 - k(n+nk1+nk2)*gnew**2*phi**2*a**2/2. + - gnew**4*phi**4*a**4/(8.*k(n+nk1+nk2)) psi2 = psi2 + int1*cc(j) c psi23 = psi23 + int1*cc(j) End Do psi2dif = (psi2-psi2old)/psi2 psi2old = psi2 beta2 = (1.d0/(2.D0*k(n+nk1+nk2)))*(Y(i+1)**2+Y(i+3)**2)+(-Y(i+1)*Y(i+2)+Y(i+3)*Y(i)) - +(k(n+nk1+nk2)/2.d0)*(Y(i)**2+Y(i+2)**2) c Write(nfile,37) t,n,k(n+nk1+nk2),Y(i),Y(i+2),int1,int2,psi2,psi2dif Write(nfile,37) t,n,k(n+nk1+nk2),Y(i),beta2,int1,int2,psi2,psi2dif Enddo !$$$c Now put in the normalization factor. psi2 = psi2/(2*pi**2*a**2) !$$$c Now add on the analytic approximation. !$$$Note that lamk is the lower limit cutoff. We called it epsilon in our paper. !$$$c Note also that when we do backreaction of the geometry the app will !$$$c be a problem that will have to be dealt with. It will have to be !$$$c computed from the backreaction equations before psi2a is computed. psi2a = -app/(48*pi**2*a**3)-gnew**2*phi**2/(8*pi**2)*(1 - dlog(2*lamk/(a*mu))) psi2 = psi2 + psi2a c Now add on the 4'th order wkb approximation for modes above the cutoff If(iback.eq.1) then psi2 = psi2 + (1.d0/(4.d0*pi**2*k(nk1+nk2+nk3)**2))*(3.d0*gnew**4*phi**4/16.d0 1 + gnew**2*phip**2/8.d0 + gnew**2*phi*phipp/8.d0) Endif RETURN END