1
$\begingroup$

Any help to integrate this equation:

$$ \begin{eqnarray} p(k,q)&=&-16 H^2 q^2\int_{q_a}^q \int_{q_a}^q \frac{1}{\sqrt{H^2-k_1{^2}}}\sinh(\sqrt{H^2-k_1{^2}}(q-b))\frac{1}{H b} \Big[m^2 y(b)^2 \nonumber\\ &-&\frac{1}{H b} y^{\prime 2}(b)\Big]+\frac{1}{\sqrt{H^2-k_2{^2}}}\sinh(\sqrt{H^2-k_2{^2}}(q-f)) \frac{1}{H b} \Big[m^2 y^2(f) \nonumber\\ &-&\frac{1}{H f} Y^{\prime 2} (f) \Big] d b ~d f \end{eqnarray} $$

where:

$y(b) = c_1 b^{\frac{1}{2}+\sqrt{1+\frac{m^2}{H^2}}}+ c_2 b^{\frac{1}{2}-\sqrt{1+\frac{m^2}{H^2}}},$

$y(f) = c_1 f^{\frac{1}{2}+\sqrt{1+\frac{m^2}{H^2}}}+ c_2 f^{\frac{1}{2}-\sqrt{1+\frac{m^2}{H^2}}}.$

Here is my naive trial:

I let $c_1=c_2=1$. Also let $k_1=k_2=k$. Let $m=H=100$ GeV. Also take $q_a=-100$.

y[b_] := c1*b^(1/2 + Sqrt[1 + m^2/H^2]) + c2*b^(1/2 - Sqrt[1 + m^2/H^2]);     

    y[f_] := c1*f^(1/2 + Sqrt[1 + m^2/H^2]) + c2*f^(1/2 - Sqrt[1 + m^2/H^2]); 

eq[k_,b_,f_,q_] := 
  -16 H^2 q^2 ( (1/Sqrt[H^2 - k1^2])*Sinh[Sqrt[H^2 - k1^2]*(q - b)]*(1/(H*b))*
       (m^2*y[b]^2 - (1/(H*b))*D[y[b], b]^2) + (1/Sqrt[H^2 - k2^2])*
       Sinh[Sqrt[H^2 - k2^2]*(q - f)]*(1/(H*b))*(m^2*y[f]^2 -  (1/(H*f))*D[y[f], f]^2) ) 

Then the integration should be taken in the limit where $q\to 0$

i1[k_?NumericQ, f_?NumericQ, q_?NumericQ] := NIntegrate[eq[k, b, f, q], {b, -100, q}]

i2[k_?NumericQ, q_?NumericQ] := 
 NIntegrate[i1[k, f, q], {f, -100, q}]

eqq[k_] = Limit[i2[k, q], q -> 0]

Plot[eqq[k], {k, 0, 1}]

But this code doesn't work!

In the equation's reference 2309.11694 $H=H^\star=10^{14}=m=-q_a$ GeV, and the plot of $p(k,q)$ as a function of $k/H$ is given by:

enter image description here

$\endgroup$
8
  • $\begingroup$ In your code functions y1,y2 are not defined. $\endgroup$ Oct 22, 2023 at 14:04
  • $\begingroup$ @AlexTrounev. You mean D[y[b], b] ? Ok, edit it. There are only y[b] and y[f]. $\endgroup$
    – Dr. phy
    Oct 22, 2023 at 14:29
  • $\begingroup$ Upper-case Y in expression is not defined. Maybe a typo. $\endgroup$
    – josh
    Oct 22, 2023 at 15:06
  • $\begingroup$ Additional context from Physics.SE: Calculating the power spectrum of the gravitational waves $\endgroup$ Oct 22, 2023 at 15:54
  • 1
    $\begingroup$ You're calling NIntegrate on an integral with symbolic parameters. You need to specify your parameter values c1 = c2 = 1; m = H = 100; , and set k1 and k2 to k. I am also not sure if Limit will work with inexact input $\endgroup$
    – ydd
    Oct 22, 2023 at 16:28

0