%Given the following bellman equation
\begin{align}
    V^i(s^i_t,S_t,D_t) =& \max_{x^i_t \geq 0} ~~ u^i(s^i_t,S_t,D_t) - Fx^i_t + \beta V^i(s^i_{t+1}, S_{t+1}, D_{t+1}) \\
    \text{Subject To:}& \notag\\
    D_{t+1} =& (1-\delta) D_t + m X_t + M\cdot \left( \sum^N_{i=1} s^i_t l^i(s^i_t,S_t,D_t) \right) + g(D_t) \label{law_motion:debris}\\
    s^i_{t+1} =& \left[ 1-l^i(s^i_t,S_t,D_t) \right] s^i_t + x^i_t \label{law_motion:private_stock}\\
    S_t =&\sum_{i=1}^N s^i_t \\
    X_t =&\sum_{i=1}^N x^i_t 
\end{align}
Where $V^i(\cdot)$ represents the value function for the constellation $i$ and $\beta$ represents a common
discount factor across operators.

\subsubsection{Characterizing solutions}
These give the optimality condition:
\begin{align}
    \frac{F}{\beta} =&  \parder{V^i}{s^i_{t+1}}{} 
                    + \parder{V^i}{S_{t+1}}{} 
                    + m\parder{V^i}{D_{t+1}}{} 
\end{align}
Iterating both forward and backward one condition gives the system
\begin{align}
    \frac{F}{\beta} =&  \parder{V^i}{s^i_{t}}{} 
                    + \parder{V^i}{S_{t}}{} 
                    + m\parder{V^i}{D_{t}}{}  \label{EQ:vi}\\
    \frac{F}{\beta} =&  \parder{V^i}{s^i_{t+1}}{} 
                    + \parder{V^i}{S_{t+1}}{} 
                    + m\parder{V^i}{D_{t+1}}{}  \label{EQ:vii}\\
    \frac{F}{\beta} =&  \parder{V^i}{s^i_{t+2}}{} 
                    + \parder{V^i}{S_{t+2}}{} 
                    + m\parder{V^i}{D_{t+2}}{} \label{EQ:viii}
\end{align}


The general envelope conditions are:
\begin{align}
    \parder{V^i}{s^i_{t}}{} - \parder{u^i}{s^i_t}{}=& \beta\left[ 
                \parder{V^i}{s^i_{t+1}}{} \parder{s^i_{t+1}}{s^i_t}{}
                + \parder{V^i}{S_{t+1}}{} \parder{S_{t+1}}{s^i_t}{}
                + \parder{V^i}{D_{t+1}}{} \parder{D_{t+1}}{s^i_t}{}
                \right] \\
    \parder{V^i}{S_{t}}{} - \parder{u^i}{S_t}{}  =& \beta\left[ 
                \parder{V^i}{s^i_{t+1}}{} \parder{s^i_{t+1}}{S_t}{}
                + \parder{V^i}{S_{t+1}}{} \parder{S_{t+1}}{S_t}{}
                + \parder{V^i}{D_{t+1}}{} \parder{D_{t+1}}{S_t}{}
                \right] \\
    \parder{V^i}{D_{t}}{} -  \parder{u^i}{D_t}{} =& \beta\left[ 
                \parder{V^i}{s^i_{t+1}}{} \parder{s^i_{t+1}}{D_t}{}
                + \parder{V^i}{S_{t+1}}{} \parder{S_{t+1}}{D_t}{}
                + \parder{V^i}{D_{t+1}}{} \parder{D_{t+1}}{D_t}{}
                \right] 
\end{align}
Note the linearity of the equations. 
This allows us to rewrite the system of envelope conditions as the following matrix expression.
\begin{align}
    \beta
    \left[
        \begin{matrix}
            \parder{s^i_{t+1}}{s^i_t}{} & \parder{S_{t+1}}{s^i_t}{} & \parder{D_{t+1}}{s^i_t}{} \\
            \parder{s^i_{t+1}}{S_t}{} & \parder{S_{t+1}}{S_t}{} & \parder{D_{t+1}}{S_t}{} \\
            \parder{s^i_{t+1}}{D_t}{} & \parder{S_{t+1}}{D_t}{} & \parder{D_{t+1}}{D_t}{} 
        \end{matrix}
    \right] 
    \left[
        \begin{matrix}
            \parder{V^i}{s^i_{t+1}}{} \\ 
            \parder{V^i}{S_{t+1}}{} \\
            \parder{V^i}{D_{t+1}}{}  
        \end{matrix}
    \right] 
    =&
    \left[
        \begin{matrix}
            \parder{V^i}{s^i_{t}}{} - \parder{u^i}{s^i_t}{} \\ 
            \parder{V^i}{S_{t}}{} -  \parder{u^i}{S_t}{}\\
            \parder{V^i}{D_{t}}{} -  \parder{u^i}{D_t}{}
        \end{matrix}
    \right] \\
\beta A  \nabla_{[s^i_{t+1},S_{t+1},D_{t+1}]} V^i =& \nabla_{[s^i_t,S_t,D_t]} V^i - \nabla_{[s^i_t,S_t,D_t]} u^i \\
    \beta A  \nabla V^i_{t+1} =& \nabla V^i_t - \nabla u_t^i ~~\text{for conciseness} 
\end{align}
Solving for $\nabla V^i_{t+1}$, we get
\begin{align}
    \nabla V^i_{t+1} =& (\beta A)^{-1} (\nabla V^i_t - \nabla u_t^i) \label{EQ:iv}
\end{align}
By iterating \eref{EQ:iv} one period, we get:
\begin{align}
    \beta A  \nabla V^i_{t+2} =& \nabla V^i_{t+1} - \nabla u_{t+1}^i  \\
    \nabla V^i_{t+2} =& (\beta A)^{-1} 
    \left( (\beta A)^{-1} (\nabla V^i_t - \nabla u_t^i)- \nabla u_{t+1}^i \right) \label{EQ:v}  
\end{align}

With \cref{EQ:iv,EQ:v} substituted into the system of \cref{EQ:vi,EQ:vii,EQ:viii}, we can
now solve for the optimal, functional form of $\nabla_{[s^i_t,S_t,D_t]} V^i$.
Substituting this back into \cref{EQ:vi} gives the euler equation for the optimal launch function
$x^i_t(s^i_t,S_t,D_t)$.

\subsubsection{Conditions for existence of a solution}
For any given set of functional forms $l^i(\cdot),g(\cdot)$ and coefficients $m,M$,
one must verify that $A$ is invertible for all
values of the state and choice variables $s^i_t,S_t,D_t$, and $x^i_t$.


For the laws of motion \cref{law_motion:debris,law_motion:private_stock}, the matrix $A$ above is:
\begin{align}
    \left[
        \begin{matrix}
            1- l^i(\cdot) - s^i_t \parder{l^i}{s^i_t}{} 
                & 1-l^i(\cdot) - s^i_t \parder{l^i}{s^i_t}{} - \sum_{j=1}^N s^j_t \parder{l^j}{S_t}{} 
                    & M\left[\parder{l^i}{s^i_t}{} + \sum^N_{j=1} \parder{l^i}{S_t}{} \right] \\
            - s^i_t \parder{l^i}{S_t}{} 
                & - \sum_{j=1}^N s^j_t \parder{l^j}{S_t}{} 
                    & M \sum^N_{j=1} \parder{l^i}{S_t}{} \\
            - s^i_t \parder{l^i}{D_t}{} 
                & - \sum_{j=1}^N s^j_t \parder{l^j}{D_t}{} 
                    & (1-\delta) + M \sum^N_{j=1} \parder{l^i}{D_t}{} + \parder{g}{D_t}{} \\
        \end{matrix}
    \right] 
\end{align}