Problem: A pebble is dropped from a large distance $H$ from the center of a planet of mass $M$ and radius $R$. How much time does it take for the pebble to strike the planet's surface? Assume we can ignore any gravitational influences on the pebble other than the planet's. Also, treat the planet as a perfect sphere with no atmosphere.
Solution: Because $H \gg R$, we cannot approximate the gravitational field as constant like we would for a mass falling near the Earth's surface. Instead, this is more of a celestial mechanics problem where the gravitational force varies with time, which makes it much more interesting. In fact, this problem is a special case of a radial trajectory using Kepler's Equation.
Start with the idea that the pebble's kinetic energy $T$ and its gravitational potential energy $V$ must add to a constant value throughout its fall. The kinetic energy is determined by the pebble's mass $m$ and its instantaneous speed $v$ by \begin{equation} T \; = \; \frac{1}{2} m v ^ 2 \label{eq:kinetic_energy} \end{equation}
We assume that the pebble's mass is very small ($m \ll M$), so that there is no detectable effect on the planet's motion.
Let $x$ be the distance between the pebble and the planet's center at any time $t$. Let $t=0$ indicate the instant in time that the pebble is released from rest. We can arbitrarily choose to define potential energy $V$ as zero when $x = \infty$. This causes $V$ to be negative at every finite distance, which may seem strange, but all that matters for this problem is how $V$ changes as the distance $x$ changes. So long as $V$ decreases at the proper rate as the pebble falls toward the planet, we are fine.
Potential energy is the integral of gravitational force multiplied by the incremental distance $\mathrm{d}r$ over which the pebble is raised or lowered in the gravitational field. The gravitational force is proportional to the planet's mass and the pebble's mass, and is inversely proportional to the square of the distance between their respective centers: \begin{equation} V \; = \; - \int_x^\infty \! \frac{mMG}{r^2} \, \mathrm{d}r \; = \; \int_\infty^x \! \frac{mMG}{r^2} \, \mathrm{d}r \label{eq:potential_energy} \end{equation} Here, $G$ is the universal gravitational constant.
So the pebble's total energy $E = T + V$ is unknown but constant: \begin{equation} E \; = \; \frac{1}{2} m v^2 + mMG \int_\infty^x \frac{\mathrm{d}r}{r^2} \end{equation} Multiply both sides by $\frac{2}{m}$ to remove the pebble's mass from the right side of the equation, and then evaluate the definite integral. \begin{equation} \frac{2 E}{m} \; = \; v^2 - \frac{2MG}{x} \label{eq:e_unknown} \end{equation} We can use initial conditions to replace the unknown constant term on the left side of \eqref{eq:e_unknown} with known values. At the moment the pebble is dropped, its speed is $v=0$ and its distance from the center of the planet is $x=H$. This gives us \begin{equation} \frac{2 E}{m} \; = \; - \frac{2MG}{H} \label{eq:e_solved} \end{equation} Substituting the right side of \eqref{eq:e_solved} back into the left side of \eqref{eq:e_unknown} and solving for $v^2$, we have: \begin{equation} v^2 \; = \; \frac{2MG}{x}-\frac{2MG}{H} \;=\; 2MG \left( \frac{1}{x} - \frac{1}{H} \right) \end{equation} We replace the speed variable $v$ with the equivalent derivative expression $\frac{\mathrm{d}x}{\mathrm{d}t}$ and take the square root of both sides to find \begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} \; = \; \pm \sqrt {2MG \left ( \frac{1}{x} - \frac{1}{H} \right ) } \label{eq:speed} \end{equation} We have reduced the problem down to two variables, $x$ and $t$. All the other values are known constants. Because $x$ decreases as the pebble falls toward the planet, we know $\frac{\mathrm{d}x}{\mathrm{d}t} \lt 0$, so we use the negative value of the square root in \eqref{eq:speed}. Solving for a tiny increment of time $\mathrm{d}t$, we have \begin{equation} \mathrm{d}t \; = \; \frac {-\mathrm{d}x} {\sqrt{2MG \left( \frac{1}{x} - \frac{1}{H} \right)}} \label{eq:delta_time} \end{equation} Answering the original question entails finding the value of $t$ when $x=R$. That value of $t$ tells us how much time elapses between the moment the pebble is dropped and the moment it strikes the surface of the planet. Taking the indefinite integral of both sides of \eqref{eq:delta_time}, we find an expression for $t$: \begin{equation} t \; = \; \int \frac {-\mathrm{d}x} {\sqrt{2MG\left(\frac{1}{x}-\frac{1}{H}\right)}} \label{eq:time} \end{equation} This is a difficult integral. I had to search through integral tables in the CRC Handbook of Chemistry and Physics (1982–1983) to figure out how to solve it. The solution required two different integrals from that table: its entries #238 and #246 are shown here as \eqref{eq:crcint1} and \eqref{eq:crcint2}, respectively: \begin{equation} \int \frac {\mathrm{d}x} {\sqrt{a + bx + cx^2}} \; = \; -\frac{1}{\sqrt{-c}}\sin^{-1}\frac{2cx+b}{\sqrt{b^2-4ac}} \; , \quad (c \lt 0) \label{eq:crcint1} \end{equation} \begin{equation} \int \frac {x \, \mathrm{d}x} {\sqrt{a + bx + cx^2}} \; = \; \frac{\sqrt{a + bx + cx^2}}{c} - \frac{b}{2c} \int \frac{\mathrm{d}x}{\sqrt{a + bx + cx^2}} \label{eq:crcint2} \end{equation} Note that \eqref{eq:crcint1} requires $c$ to be negative; if $c \gt 0$ then an alternate solution for the integral applies. Also note that the rightmost term of \eqref{eq:crcint2} requires use of \eqref{eq:crcint1}. Let's combine these two equations so that the right side of \eqref{eq:crcint2} has no integrals in it. \begin{equation} \int \frac {x \, \mathrm{d}x} {\sqrt{a + bx + cx^2}} \; = \; \frac{\sqrt{a + bx + cx^2}}{c} + \frac{b}{2c} \frac{1}{\sqrt{-c}}\sin^{-1}\frac{2cx+b}{\sqrt{b^2-4ac}} \label{eq:combined} \end{equation}
But how does \eqref{eq:combined} help us evaluate \eqref{eq:time}? We need to do a little algebraic manipulation of \eqref{eq:time} to put it in the right form. \begin{equation} t \; = \; \frac{1}{\sqrt{2MG}} \int \frac{-\mathrm{d}x}{\sqrt{\frac{H-x}{Hx}}} \; = \; - \sqrt{\frac{H}{2MG}} \int \sqrt{\frac{x}{H-x}} \mathrm{d}x \label{eq:manip} \end{equation} We are still not there yet. Because the pebble never reaches the center of the planet, $x \gt 0$ at all times. Therefore it is safe to multiply the integrand by $\frac{\sqrt{x}}{\sqrt{x}}$, giving us \begin{equation} t \; = \; -\sqrt{\frac{H}{2MG}} \int \frac{x \, \mathrm{d}x}{\sqrt{Hx-x^2}} \label{eq:goodt} \end{equation} Now we have something we can work with. The integral in equation \eqref{eq:goodt} matches the left side of \eqref{eq:combined} when we use the following constant substitutions. \begin{align*} a & = 0 \\ b & = H \\ c & = -1 \end{align*} At this point we confirm that $c \lt 0$ as required by \eqref{eq:crcint1}. Replacing the generic constants $a$, $b$, and $c$ from \eqref{eq:combined} with our substitutions gives us \begin{equation} \int \frac{x\,\mathrm{d}x}{\sqrt{Hx - x^2}} \;=\; \frac{\sqrt{Hx-x^2}}{-1} + \frac{H}{-2} \frac{1}{\sqrt{1}} \sin^{-1} \frac{-2x+H}{\sqrt{H^2-0}} \\ \end{equation} \begin{equation} \int \frac{x\,\mathrm{d}x}{\sqrt{Hx - x^2}} \;=\; -\sqrt{Hx - x^2} - \frac{H}{2} \sin^{-1} \frac{H - 2x}{H} \label{eq:intsimp} \end{equation} Looking back at \eqref{eq:goodt}, the expression for $t$ needs the extra negative factor in front of its integral. Multiplying that factor cancels out the negatives on the right side of \eqref{eq:intsimp} to produce \begin{equation} t \; = \; \sqrt{\frac{H}{2MG}} \left[ \frac{H}{2} \sin^{-1} \left( \frac{H-2x}{H} \right) + \sqrt{Hx-x^2} \right] + t_C \label{eq:elimint} \end{equation} We are not done yet because now we have this extra unknown constant $t_C$ from the indefinite integral. (Integral tables omit the arbitrary constant for conciseness, but we can't ignore it here!) We eliminate $t_C$ using initial conditions. When the pebble is dropped, $x=H$ and $t=0$: \begin{equation} 0 \; = \; \sqrt{\frac{H}{2MG}} \left[ \frac{H}{2} \sin^{-1}(-1) + \sqrt{0} \right] + t_C \end{equation} \begin{equation} 0 \; = \; \sqrt{\frac{H}{2MG}} \left[ - \frac{H}{2} \frac{\pi}{2} \right] + t_C \end{equation} \begin{equation} t_C \; = \; \sqrt{\frac{H}{2MG}} \frac{H}{2} \frac{\pi}{2} \label{eq:t0} \end{equation} Substitute the right side of \eqref{eq:t0} back into \eqref{eq:elimint} to obtain: \begin{equation} t \; = \; \sqrt{\frac{H}{2MG}} \left[ \frac{H}{2} \left( \frac{\pi}{2} + \sin^{-1} \left( \frac{H-2x}{H} \right) \right) + \sqrt{Hx - x^2} \right] \label{eq:arcsin} \end{equation} We can simplify \eqref{eq:arcsin} using the identities \begin{equation} \sin^{-1}(-u) \; = \; -\sin^{-1}(u) \end{equation} \begin{equation} \cos^{-1}(u) \; = \; \frac{\pi}{2} - \sin^{-1}(u) \end{equation} Also, I think factoring out $x$ inside the square root makes that part of the equation more intuitive: $x$ is the distance from the center of the planet, and $H-x$ is how far the pebble has dropped. \begin{equation} t \; = \; \sqrt{\frac{H}{2MG}} \left[ \frac{H}{2} \cos^{-1} \left( \frac{2x-H}{H} \right) + \sqrt{x(H - x)} \right] \label{eq:arccos} \end{equation} This tells us the time $t$ when the pebble reaches any position $x$ relative to the planet's center. To find the impact time $t_R$, we let $x=R$, which gives us the answer we were looking for: \begin{equation} t_R \; = \; \sqrt{\frac{H}{2MG}} \left[ \frac{H}{2} \cos^{-1} \left( \frac{2R-H}{H} \right) + \sqrt{R(H - R)} \right] \label{eq:solution} \end{equation}
This concludes the solution of the problem. A different form of this formula appears in this Wikipedia article as the radial free fall of two objects. I find it interesting that the inverse cosine term appears in this formula. It suggests that the time might conceptually include some angular value multiplied by $\frac{H}{2}$, although it is not clear what that might mean.
Numerical Verification: To verify that this formula works correctly, I wrote a C++ program called pebble.cpp that accepts a height above the surface (or $H-R$, from which it calculates $H$) as a command line parameter. The program calculates the impact time using \eqref{eq:solution}, and then it runs a numerical simulation of the pebble falling, approximating the impact time by breaking the fall into thousands of small increments of distance and time. The program prints both time calculations for comparison.
Here are some points to check that both the formula and the program work properly:
Feedback: If you have any questions or comments, please feel free to contact me: