No, that setup is not stable

There do not appear to be any stable orbits for another moon, of any size, outside the orbit of Luna

I ran a simulation in rebound to see what would happen with this moon setup. I used a grid of different behaviors to check how the second moon would react in different scenarios.

Here was my setup in rebound:

import rebound
from math import sqrt

for m_selene in [3/2, 1, 1/2, 1/4, 1/10, 1/100, 1/1000]:
    for a_selene in [3/2, 2, 3, 5, 10]:
        for e_selene in [0, 0.01, 0.02, 0.05, 0.1, 0.25]:
            sim = rebound.Simulation()
            sim.integrator = 'whfast'
            sim.units = ('AU', 'days', 'Msun')


            sim.add(m=1)
            sim.add(m=0.000003004, a=1, e=.016709)
            sim.add(primary=sim.particles[1], m=0.000000037, a=0.00257, e=.0549)
            sim.add(primary=sim.particles[1], m=0.000000037*m_selene, a=0.00257*a_selene, e=e_selene)

            sim.integrate(10*365.2563)
            earth_luna = sqrt((sim.particles[2].x - sim.particles[1].x)**2 + (sim.particles[2].y - sim.particles[1].y)**2)
            earth_selene = sqrt((sim.particles[3].x - sim.particles[1].x)**2 + (sim.particles[3].y - sim.particles[1].y)**2)
            print(m_selene, a_selene, e_selene, "{0:.2f}".format(earth_luna), "{0:.2f}".format(earth_selene))

First I ran the simulation for just the earth-moon system to ensure it was stable, and it was. Then I ran it again with all objects (Sun, Earth, Luna "” the real moon, Selene "” a new outer moon).

I modified the above program to find 'breakpoints' where different behavior was observed. This is what I got:

• Let m_selene be the mass of Selene as a proportion of Luna's mass (so 0.5 means Selene is half the size of the moon).
• Let a_selene be the semi-major axis of Selene as a proportion of Luna's semi-major axis (so 2 means twice as far from the Earth as Luna).
• Let e_selene be the eccentricity of Selene.

Then there are four scenarios.

• Luna is ejected from the system, and Selene becomes the new moon. This occurs only in rare circumstances when Selene is large and close to Luna (m_selene > 1.25, a_selene < 1.75) or when Selene is medium size but with a highly eccentric orbit (m_selene > 0.35, 1.25 < a_selene < 2.25, e_selene > 0.075).

• Luna and Selene are both ejected from the system. This is likely when Luna and Selene are about the same size and close to each other (m_selene > 0.35, a_selene < 2.25). This happens in both highly eccentric and non-eccentric orbits of Selene. It can also happen rarely for smaller masses of Selene (as low as m_selene > .125) but not for larger orbits.

• Luna and Selene are (somehow?!?!) ejected from the solar system. I don't see how this is possible, but it occurs in the case of m_selene = 0.5, a_selene = 1.5, and the eccentricities of selene and the moon are about the same.

• The most common situation is that Selene is ejected from the system and Luna remains. This happens in every (actually all but one) case where m_selene < 0.35 and a_selene > 2.25, and rarely when Selene is larger or closer.

So what you do not notice is any situation where Selene and Luna are both stable in orbit around earth. So this is not an exhaustive numerical solution, nor is it an analytical proof, but from a grid search of a numerical solution, I cannot find any stable orbits for a second smaller moon outside the orbit of Luna.

---

Edit: In response to @Molot's request for orbits inside of Luna, it turns out that most orbits inside Luna's are stable for at least a short time. The grid search only integrates 10 years, so that doesn't prove long term stability, but I ran a simulation for:

• m_selene = 0.5
• a_selene = 0.5
• e_selene = 0.02

and it was stable for 100,000 years. The moon moved slightly outwards over the first 1000 years so that Luna and Selene achieved a resonance, but then the two were stable for as long as the simulation was running.

I can't promise stability for billions of years, and to be truly accurate I'd have to account for Jupiter and Saturn (at least) and the early solar system planet migrations, but it seems plausible that Selene could exist in an orbit inside that of Luna.

Note: If anyone can double-check my numbers and result at the end of this answer, that would be much appreciated.

This looks like a question that can be broken down as per your bullet points, so I think I'll do it like that.

Is earth's gravity strong enough to have a second moon further from us than our current moon?

As JDÅ‚ugosz pointed out, the answer to this depends on whether or not the moon is inside Earth's Hill sphere, the region in which it can hold satellites in stable orbits. The general formula for a body of mass $m$ orbiting a body of mass $M$ with a semi-major axis of $a$ is $$r\approx a\sqrt[3]{\frac{M}{3m}}$$ For Earth, this is^[1] $$r_{\oplus}\approx a_{\oplus}\sqrt[3]{\frac{M_{\odot}}{3m_{\oplus}}}=1.496\times10^{9}\text{ meters}$$ This is roughly one tenth of the distance to the Sun (and four times the Moon's orbital radius). It doesn't take into account influences from the other planets, but that shouldn't be an issue. After all, Venus and Mars are both much further than 0.1 AU away from Earth, even at their closest passes.

Regarding your request for different angular velocities: Two bodies in circular orbits with different radii will always move at different angular and tangential velocities, so you're fine there.

Is there a certain size it would have to be?

So long as the mass of the extra moon is much less than that of the planet, there should be no ill effects on either Earth or the Moon. Any non-zero mass will perturb the orbit of the Earth a bit, so there's no cutoff line. You just have to specify what limit is okay with you.

It's interesting to think of the effects of the moon on the Moon. I have something to say on the subject - something that I've been meaning to write about for a while - but I'll make that a separate section.

Can an orbit like this be stable? (two moons, one half the size of the moon and moving at different speeds, relative to earth)

Well, yes, it can - again, so long as the moon is far enough away from the Moon. The places where the setup are unstable can be found easily - again, this will come later. The short answer, though, is that this should be perfectly safe. Just use Newton's law of gravitation: $$F=G\frac{m_{\text{Moon}}\frac{1}{2}m_{\text{Moon}}}{(r_{\text{Moon}}-r_{\text{moon}})^2}$$ This gives you the force at closest approach - which should be minimal, for large enough values of $r_{\text{moon}}$.

---

Here's the section on stability - the stability of the Moon's orbit - which I alluded to earlier. I'm going to base it strongly off of these lecture notes.

Here, we use the disturbing function, which has the same magnitude as the gravitational potential from the perturbing¹ body - in this case, the extra moon. The general form of it is $$\mathcal{R}=-\frac{Gm_{12}m_3}{R}+\frac{Gm_2m_3}{|\mathbf{R}-\alpha_1\mathbf{r}|}+\frac{Gm_1m_3}{|\mathbf{R}+\alpha_2\mathbf{r}|}\tag{1}$$ where the subscripts $_1$, $_2$, and $_3$ refer to Earth, the Moon, and the moon, respectively, $m_{ij}=m_i+m_j$, $\alpha_i=m_i/m_{12}$.

To analyze this, we have to expand it - basically, write it in another form using summations and functions called spherical harmonics. For those interested, the lecture notes provide a full derivation, but to make it short, the intermediate expansion is $$\mathcal{R}=G\mu_im_3\sum_{l=2}^{\infty}\sum_{m=-l}^l\left(\frac{4\pi}{2l+1}\right)\mathcal{M}_l\left(\frac{r^l}{R^{l+1}}\right)Y_{lm}(\theta,\varphi)Y_{lm}^*(\Theta,\psi)\tag{2}$$ where $$\mathcal{M}_l=\frac{m_1^{l-1}+(-1)^lm_2^{l-1}}{m_{12}^{l-1}}$$ and $Y_{ab}(\beta,\gamma)$ is a spherical harmonic.

Why do we care? Good question. The disturbing function allows us to identify orbital resonances², which can help a system stay stable.

To do this, we have the expand the disturbing function again, this time as a Fourier expansion. We get $$\mathcal{R}=G\mu_im_3\sum_{l=2}^{\infty}\sum_{m=-l}^{l}c_{lm}M_l\left(r^le^{imf_i}\right)\left(\frac{e^{-imf_o}}{R^{l+1}}\right)e^{im(\varpi_i-\varpi_o)}\tag{3}$$ where $$c_{lm}=\frac{(l-m)!}{(l+m)!}[P_l^m(0)]^2$$ where $P_l^m(x)$ is an associated Legendre polynomial and $f$ and $\varpi$ are orbital elements.

This can all be expanded further by way of . . . am I boring you? Well, I'll skip to the good bit. I just wanted to emphasize that this one little function contains a lot of information.

We eventually arrive at a function $\phi_{n'nm}$ that is called the resonance angle, given by $$\phi_{n'nm}=n'M_i-nM_o+m(\varpi_i-\varpi_o)\tag{4}$$ $M$ is the mean anomaly and $\varpi$ is the argument of periapsis.

We can find the rate of change of this with respect to time, $\dot{\phi}$, and then write a function $\mathcal{E}$ as $$\mathcal{E}=\frac{1}{2}\dot{\phi}^2-\omega^2(\cos\phi+1)$$ There are three possible classes of values of $\mathcal{E}$:

1. $\mathcal{E}>0$: The system is circulating.
2. $\mathcal{E}<0$: The system is librating.
3. $\mathcal{E}=0$: The system is at an unstable equilibrium (think of a pendulum pointed straight up).

We can then write down another equation where $$\ddot{\phi}\propto\sin\phi$$ The associated $\mathcal{E}$ tells us if the system is stable.

There's an explicit stability algorithm by the author of the lecture notes in a related paper (see Mardling (2008) using a different expansion. It should be able to give you a yes-or-no answer for resonant stability for given orbits of the Moon and the moon. I'll see if I can apply it here.

The precise algorithm is

(1) Identify which $[n : 1](2)$ resonance the system is near and calculate the distance $\delta_{\sigma n}$ from that resonance: $\delta_{\sigma n} =\sigma-n$, where $n = \lfloor\sigma\rfloor$. (the nearest integer for which $n \leq \sigma$);
(2) Take the associated resonance angle to be zero rather than the definition (2.18) (see discussion below): $\phi_{2n1} = 0$;
(3) Calculate the induced eccentricity from (5.1) and (if $m_1 = m_2$) the maximum octopole eccentricity from (5.3). Determine $e_i = \text{max}[e^{(ind)}_i, e^{(oct)}_i]$ for use in $s^{(22)}_1 (e_i)$;
(4) Calculate $\mathcal{A}_{2n1}$ from (3.6);
(5) Calculate $\mathcal{E}_{2n1}$ and $\mathcal{E}_{2 \ n+1\ 1}$: deem the system unstable if $\mathcal{E}_{2n1} < 0$ and $\mathcal{E}_{2 \ n+1\ 1} < 0$.

I've run some numbers, and gone through a resonance or two. Mardling focused on $[n:1](2)$ resonances, which are important here. I found that for a $[2:1](2)$ resonance, the system is unstable at most eccentricities, but it seems - judging from some of the graphs given - that at higher resonances, the system is perfectly stable.

I ran through a certain combination of values to test for stability; it would be awesome of someone could check them. Here they are, with intermediate values (the subscripts $_i$ and $_o$ refer to the inner and outer satellites, the Moon and the moon).

The given data are (I've chosen the last two): $$m_1=m_{\text{Earth}}=5.9722\times10^{24}\text {kg}=81.285 m_2$$ $$m_2=7.3462\times10^{22}=1m_2$$ $$m_3\equiv\frac{1}{2}m_2$$ $$e_i(0)=0.0549$$ $$e_o(0)=0.01$$ $$a_o=2.75a_i$$ The intermediate values are: $$\sigma=4.5466$$ $$n=\lfloor\sigma\rfloor=4$$ $$\delta\sigma_4=\sigma-4=0.5466$$ $$n=4,n'=1,m=2$$ $$e_i^{(eq)}=0.02206$$ $$A=1.4891$$ $$e_i^{(oct)}=0.09902$$ $$l=l_{min}=2$$ $$s_{224}=0.44833$$ $$F_4^{(22)}(e_o)=2.712\times10^{-6}$$ $$f_4^{(22)}(e_o)=2.6513\times10^{-6}$$ $$\beta_4=5.6513\times10^{-8}$$ $$e_i^{(ind)}=0.0549$$ $$e_i=\text{max}[e_i^{(ind)},e_i^{(oct)}]=e_i^{(oct)}=0.09902$$ $$s_4^{(22)}(e_i)=-0.16443$$ $$\mathcal{M}_2=1$$ $$M_i^{(2)}=0.006114$$ $$M_o^{(2)}=0.0119568$$ $$\sigma_4=1$$ $$c_{(22)}=\frac{3}{8}$$ $$\mathcal{A}_{241}=3.9059\times10^{-8}$$ $$\delta\sigma_{41}=\sigma-n/n'=0.5466$$ $$\bar{\mathcal{E}}_{241}=0.1493857019$$ This is greater than zero, so the system is stable.

---

Update, 28 October 2016

I've been going through and checking some of these numbers again, in the process of trying to write an algorithm around this, and I'm fairly certain that there are some errors, potentially including (and thus starting with) $s_4^{(22)}(e_i)$. This may mean that the system is unstable, as kingledion's results show.

I've been running some simulations of my own, using the European Space Agency's open source ORSA software. I've only done four or five runs over shot periods of time, but they're starting to fall into two categories: Either the secondary moon is slowly ejected, or the two go into separate elliptical orbits, with periodic changes in eccentricity and semi-major axes. For the sake of accuracy, I've only run the simulations over one year, but it seems like there may be stable setups and there may be unstable setups.

Additionally, I haven't included the Sun in this, so I don't know how it could change any of this.

---

^{1 In this case, I'm talking about perturbing the Moon.
2 Specifically, orbit-orbit resonances.}