states. The latter is performed by solving the Bogoliubov equations for the linear excitations $[u(z,t),v(z,t)]$ around the dark-dark soliton states $\psi(z,t)=\exp(-i \mu t)[\psi(z) +\sum_\omega(u \,e^{-i\omega t}+ v^* 
e^{i\omega t})]$. These equations are:

$$\begin{aligned}
B 
\left(\begin{array}{c}
u_1 \\
v_1 \\
u_2 \\
v_2
\end{array}\right)
= \omega 
\left(\begin{array}{c}
u_1 \\
v_1 \\
u_2 \\
v_2
\end{array}\right),
\label{eq:BdG}\end{aligned}$$

where: $$\begin{aligned}
B= 
\left(\begin{array}{cccc}
h_1 & g\psi_1^2 & g_{12}\psi_ 1 \psi_2^* & g_{12}\psi_1 \psi_2 \\
-g\psi_1^{*2} & -h_1 &  -g_{12}\psi_1^* \psi_2^* &  -g_{12}\psi_1^* \psi_2 \\
g_{12}\psi_1^* \psi_2 &  g_{12}\psi_1 \psi_2 & h_2 &  g\psi_2^2 \\
-g_{12}\psi_1^* \psi_2^* & - g_{12}\psi_1 \psi_2^* &-  g\psi_2^{*2} & -h_2
\end{array}\right).\end{aligned}$$ Here $h_i=-\frac{1}{2} \frac{\partial^2}{\partial z^2} + \frac{1}{2}\Omega^2z^2 + 2g |\psi_i|^2 + g_{12} |\psi_{j\neq i}|^2 - \mu_i$. They can be seen as an eigenvalue problem with a non-trivial solution given by ${\rm det}|B- \omega I|=0$.

Apart from the perturbative case $g_{12}<<g$, we also explore numerically the stability of overlapped solitons in the whole miscible regime $g_{12}<g$, and show that the finite size of the system determines the stability properties.

Untrapped case
--------------

![(Color online) Spatio-temporal evolution of the density of two untrapped dark-dark solitons with $\mu=1$, $g_{12}=0.015$ and $g=1$ obtained from the numerical solution of the GP Eq. . The top 