How fast are planets flowing into stars? I got the equations.

[astrostuart]

I am an astronomer with a research question:
I want to evolve this equation for planet distribution:
<br /> \frac{{\operatorname{d}}f(P)}{{\operatorname{d} \operatorname{log}}{P}}=<br /> k_P P^\beta \left(1-e^{-(P/P_0)^\gamma }\right)<br />
as a function of period ``P'',
by using an equation for change of P.

I will actually use
this equation of change of semi-major axis ``a''
(below)
for which I will input assumed values of the
masses and radii of the star and planets, (but transforming between P and a
is not hard. The calculus confuses me, though.)
<br /> \frac{1}{a}\frac{da}{dt}=<br /> -\left(\frac{9}{2}{\left(G/M_*\right)}^{1/2}\frac{R_p^5M_p}{Q_*^{prime}}\right)<br /> a^{-13/2}<br />

<br /> G M P^2= (2 \pi)^2 a^3<br />

I know it's unusual to think of planets "flowing" into stars, but consider this to be in the large limit of what happens to so many planets, that each planet is like a piece of fluid. Perhaps a mathematician can suggest how to say this better.

Thanks

Stuart

 

[meldraft]

What exactly is your question? Is it how to solve the system of equations? Have you tried plugging your system into a numerical solver? It is highly non-linear so I don't see much hope in deriving a closed form solution.

 

[astrostuart]

I am trying to find the rate of planets flowing past radius of the star, R_*.
This would give how many planets flow into the star:

<br /> \frac{\operatorname{d}f(P)}{\operatorname{d} t }<br />
evaluated at the radius of the star, R_*.
(In the above equation, there should be R_*^5,
the radius of the star, not the radius of the planet R_p^5)

so I need to put the tidal migration equation (2nd equation)
into the distribution equation (1st).

Let me try to simplify the problem a little, by first deriving the tidal
migration equation as a function of period "P" instead of the semimajor axis "a", and after that
we can start with just one of the two terms in the occurrence
equation.

I worked to solve for <br /> \frac{\operatorname{d}P}{\operatorname{d} t }<br /> from <br /> \frac{\operatorname{d}a}{\operatorname{d} t }<br />

using
<br /> a=(GM_*)^{1/3}(\frac{ P }{ 2 \pi})^{2/3}

and thus <br /> \frac{\operatorname{d}a}{\operatorname{d}P }=<br /> \frac{2}{3} \frac{ (GM_*)^{1/3} }{ (2 \pi)^{2/3} } \frac{1}{P^{1/3} }<br />
<br /> \frac{\operatorname{d}P}{\operatorname{d}a }=<br /> \frac{3}{2} \frac{ (2 \pi)^{2/3} }{ (GM_*)^{1/3} }P^{1/3} <br />

I get

<br /> \frac{\operatorname{d}P}{\operatorname{d}t }=<br /> \frac{-27}{32 \pi} \frac{ G^{5/3} }{ M_*^{8/3} } <br /> \frac{R_s^5 M_p}{Q_*} P^{-11/3} <br />

I would be grateful if someone could follow along and check this.
Yes, this would be good to put into a symbolic solver.

The next step would be the change in the distribution as a function of time.

 

[astrostuart]

Now, I want P as a function of T using the previous derivation of <br /> \frac{\operatorname{d}P}{\operatorname{d}t }
(I accidentally put an s subscript for the radius of the star above, should be "*")

<br /> \operatorname{d}P=<br /> \frac{-27}{32 \pi} \frac{ G^{5/3} }{ M_*^{8/3} } <br /> %\frac{R_*^5 M_p}{Q_*} P^{-11/3} {\operatorname{d}t }<br />

Can I simply integrate as below? But where do I put the initial distribution,
which I now call P_0, but this is not the same as the badly named
P_0 in the top equation of the first
post, which should really be named P_f for the falloff location.

Integrating, I get
P=- \frac{ 189 \pi^3 }{48} <br /> \frac{ G^{5/3} }{ M_*^{8/3} } <br /> \frac{R_*^5 M_p}{Q_*} P^{-11/3} t <br />

But where would I put the P_0?
I hope I correctly kept track of the terms.

 

[astrostuart]

I think I have an expression on the flow rate, so now my question becomes did I do it right? Is the concept right? Did I keep track of all the terms?

I will put the above dP/dt into the distribution equation using the chain rule:
First, remove the logarithm from the distribution
<br /> \frac{{\operatorname{d} \operatorname{log}}{P}}{{\operatorname{d}}f(P)}=<br /> \frac{1}{P},
so the distribution is (writing P_f instead of P_0),

<br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d} \operatorname{log}}P }=<br /> \frac{1}{P}k_P P^{\beta} \left(1-e^{-(P/P_0)^\gamma }\right)<br />
<br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d}}P }=<br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d} \operatorname{log}}P }<br /> \frac{{\operatorname{d} \operatorname{log}}{P}}{{\operatorname{d}}f(P)}<br />
<br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d}}P }=<br /> k_P P^{\beta-1} \left(1-e^{-(P/P_0)^\gamma }\right)<br />

Using the chain rule <br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d}}t}=<br /> \frac{ {\operatorname{d}}f(P) }{ \operatorname{d} t}<br /> \frac{\operatorname{d} f(P)}{\operatorname{d}f(P)}<br /> gives,

<br /> \frac{ {\operatorname{d}}f(P) }{ {\operatorname{d}}P }=<br /> -\frac{27}{32 \pi^3} \frac{5^{5/3}}{M_s^8/3}\frac{R_s^5 Mp}{Q_s}<br /> k_P P^{\beta-{14/3}} \left(1-e^{-(P/P_0)^\gamma }\right)<br />

I hope that this evaluated for P at the stellar radius gives the rate of planet infall.
That is, P can be found from taking a=R_* in
<br /> G M P^2= 2 \pi a^3<br />
so we would use this "P_*" for P:
<br /> P_*= \left(\frac{2 \pi}{GM}\right)^{1/2} a^{3/2}<br />

If this is right, it is the answer to my original question.
I would like to go further, however, to see if I could get a closed form expression
for how P_f, \beta, and \gamma evolve under the effect of tidal migration equation.
Actually, the one that will change meaningfully is P_f, as \beta changes little, and \beta+ \gamma go to 14/3.
I already have been doing this numerically, and it seems that P_f pretty much increases at an exponentially decreasing rate.
See Taylor (2012) at http://arxiv.org/abs/1211.1984
and Taylor (2013) at http://arxiv.org/abs/1301.4229

I hope the closed form solutions can produce similar results!

 

[meldraft]

I'm a little confused about Pf and P0. If I understood correctly, where I now read P0, it should be Pf? If so, P0 is probably the initial condition you will use to solve your equation. Your final equation is pretty easy to solve (unless I'm missing something), considering it's full of constants and the variables are separated. It degenerates into two integrals, one trivial ( \int constants*P^{constants}dP ), while the other one is of the form \int e^{-P^n}dP. This cannot be solved in terms of elementary functions, but its solution can be defined using the gamma function. You can find its solution here:

http://www.wolframalpha.com/input/?i=int(exp(-x^n),x)

 

[astrostuart]

My mistake, I did not think how P_f for "P falloff" would be confused with "P final", so perhaps P_c for "P change" or "P crossover" since this variable sets where the distribution goes from a steep slope at shorter periods than the "crossover period" to a less steep slope at longer periods.

The form of the distribution function is merely contrived to fit the data, and the only part that matters are for the trivial case without the exponential. The term with the exponential should be taken in the limit of small periods, that is, it will be the same form as the first term, but with a different slope of \beta+\gamma replacing \beta. I would appreciate someone checking that the (now renamed) P_c would get put into the constant.

I appreciate those that have tried to follow this; I can see that I need to do a rewrite using what I've learned here to make it more clear. The fact that it wasn't entirely clear to myself in the beginning is why I came to this forum.