Growing a set given constraints
I suggest you formulate this as an instance of integer linear programming (ILP).  Let $m$ be the largest number in $S_0$.  For each $i$ such that $1 \le i \le m$, introduce the zero-or-one variable $x_i$, with the intended meaning that $x_i=1$ means that $i \in S$ and $x_i=0$ means that $i \notin S$.

Now your constraints can be readily converted into linear inequalities: e.g., for each $s$ such that $2&lt;s \le m$, we obtain that $x_s=1$ implies $\lor_{s_a,s_b} (x_{s_a}=1 \land x_{s_b}=1)$, where the disjunction is taken over all $s_a,s_b$ such that $1 \le s_a &lt; s$ and $1 \le s_b &lt; s$.  This can be converted into a linear inequality; see Express boolean logic operations in zero-one integer linear programming (ILP).  Also, for each $i \in S_0$, we add the requirement $x_i=1$.  Finally, we minimize $\sum_i x_i$.  Feeding this to an off-the-shelf ILP solver should yield a solution.

This should find the exact optimal solution, as long as the largest number in $S_0$ is not too large.  However, the running time is potentially exponential in $m$, the largest number in $S_0$.  I don't know if there is a polynomial-time solution.

(You could also formulate it as an instance of SAT instead of ILP.  This will require applying the Tseitin transform to convert the implication/disjunction into CNF, and it will require constructing an adder-circuit to add the requirement that the size of $S$ is at most $k$ and then doing binary search over $k$.  Then, you could feed it to an off-the-shelf SAT solver.  I have no idea whether this will work better than the ILP approach.)