I think $|\Psi\rangle_S$ and $|\Psi\rangle_I$ are used as follows. First $H$ is split into the two parts:
$$H = H_0 + H_1$$
where $H_1$ contains the interesting interacting parts and $H_0$ are to be eliminated wherever possible.
Then we have
$$ \boxed{\color{blue}{i\frac{d\Psi}{dt}} = H\Psi = (H_0 + H_1)\Psi}$$
Now consider the quantity $|\Psi\rangle_I = e^{iH_0t}|\Psi\rangle_S=e^{iH_0t}\Psi$:
$$ i\frac{d}{dt}\left(\color{red}{e^{iH_0t}\Psi}\right) = i \frac{d}{dt}(e^{iH_0t}\Psi) = \color{orange}{-H_0e^{iH_0t}\Psi} + \color{orange}{e^{iH_0t}H_0\Psi} + e^{iH_0t}H_1\Psi = H_1 \left(\color{red}{e^{iH_0t} \Psi }\right)$$
Solving this first order differential equation for $e^{iH_0t}\Psi$,
$$ i\frac{d}{dt}\color{red}{e^{iH_0t}\Psi} = H_1\color{red}{e^{iH_0t}\Psi}$$
$$ e^{iH_0t}\Psi = -i\int_{t=0}^{t1} H_1e^{iH_0t}\Psi dt + c$$
By physical considerations,
$$\boxed{c = 1}$$
$W$ on the right hand side can be integrated similarly up to some $t_2$ less than $t_1$, and again ad infinitum:
$$ = 1 - i\int_{t=0}^{t_1} H_1 \left(\color{purple}{1 - i\int_{t_2=0}^{t_2} H_1e^{iH_0t}\Psi dt_2}\right) dt_1 $$
$$= 1 - i\int_{t_1=0}^{t_1} H_1(t_1)dt_1 - H_1(t_1) i\int_{t_1=0}^{t1}\int_{t_2=0}^{t_2} H_1(t_2) e^{iH_0t}\Psi dt_2 dt_1$$
$$= 1 - i\int_{t_1=0}^{t_1} H_1(t_1)dt_1 - H_1(t_1) i\int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2} H_1(t_2) \left(\color{purple}{1 - i\int_{t_3=0}^{t_2} H_1e^{iH_0t}\Psi dt_3}\right) dt_2 dt_1$$
$$= 1 + (- i)\int_{t_1=0}^{t_1} H_1(t_1)dt_1 + (-i)^2\int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2}H_1(t_1) H_1(t_2) dt_2 dt_1 + (-i)^3\int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2}\int_{t_3=0}^{t_3} H_1(t_1) H_1(t_2) H_1(t_3) e^{iH_0t}\Psi dt_3 dt_2 dt_1$$
In general,
$$ e^{iH_0t}\Psi = \sum_{n=0}^\infty (-i)^n \int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2}\int_{t_3=0}^{t_3}...\int_{t_n=0}^{t_n}H_1(t_1)H_2(t_2)H_3(t_3)...H_n(t_n)$$
Term $(-i)^2$:
$$ T_2 = \int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2}H_1(t_1) H_1(t_2) dt_2 dt_1 $$
$$ = \int_{t_1=0}^{t}\int_{t_2=0}^{t}\theta(t_1-t_2) H_1(t_1) H_1(t_2) dt_2 dt_1 $$
where $ \theta(t_1-t_2) = 1$ only if $t_1-t_2 > 0$.
Term $(-i)^3$:
$$ T_3 = \int_{t_1=0}^{t_1}\int_{t_2=0}^{t_2}\int_{t_3=0}^{t_3} H_1(t_1) H_1(t_2) H_1(t_3) dt_3 dt_2 dt_1 $$
$$ \int_{t_1=0}^{t}\int_{t_2=0}^{t}\int_{t_3=0}^{t} \theta(t_1-t_2)\theta(t_2-t_3) H_1(t_1) H_1(t_2) H_1(t_3) dt_3 dt_2 dt_1 $$
since $\theta(t_1 - t_2)\theta(t_2 - t_3) = 1$ only if $t_1 > t_2 $ and $ t_2 > t_3$.
In general,
$$ \small{e^{iH_0t}\Psi = \sum_{n=0}^\infty (-i)^n \int_{t_1=0}^{t}\int_{t_2=0}^{t}...\int_{t_n=0}^{t}\theta(t_1-t_2)\theta(t_2-t_3)...\theta(t_{n-1}-t_n)H_1(t_1)H_2(t_2)...H_n(t_n)}$$
is numerically solvable to high summation orders.

$$\boxed{\Psi_{x<0} = ae^{ikx} + be^{-ikx}}$$
$$\boxed{\Psi_{x>a} = fe^{ikx} + ge^{-ikx}}$$
since
$$ -\frac{d^2\Psi}{dx^2} + (0)\Psi$$
$$ = -\frac{d^2}{dx^2}(ae^{ikx} + be^{-ikx})$$
$$ = -\frac{d}{dx}(ikae^{ikx} + (-ik)be^{-ikx})$$
$$ = -(ik)^2ae^{ikx} - (-ik)^2be^{-ikx}$$
$$ = k^2ae^{ikx} + k^2be^{-ikx}$$
$$ = k^2\Psi = E\Psi$$
for any mixtures of $a,b,f,g$.
$$\boxed{\Psi_{0<x<a} = ce^{\kappa x} + de^{-\kappa x}}$$
since
$$ -\frac{d^2\Psi}{dx^2} + (V_0)\Psi$$
$$ = -\frac{d^2}{dx^2}(ce^{\kappa x} + de^{-\kappa x}) + V_0\Psi$$
$$ = -\frac{d}{dx}(\kappa ce^{\kappa x} + (-\kappa)de^{-\kappa x}) + V_0\Psi$$
$$ = -(\kappa)^2ce^{\kappa x} - (-\kappa)^2de^{-\kappa x} + V_0\Psi$$
$$ = (-\kappa ^2+V_0)\Psi = E\Psi$$
since the lack of $i$ reduces the $E =-\kappa^2+V_0$ requirement.
Ensuring continuity at boundary $x=0$:
$$\Psi_{x<0} = \Psi_{0<x<a}$$
$$ae^{ik(0)} + be^{-ik(0)} = ce^{\kappa(0)} + de^{-\kappa(0)}$$
$$\boxed{a+b = c+d}$$
and
$$ \frac{d}{dx}\Psi_{x<0} = \frac{d}{dx}\Psi_{0<x<a}$$
$$ \frac{d}{dx}(ae^{ikx} + be^{-ikx}) = \frac{d}{dx}(ce^{\kappa x} + de^{-\kappa x})$$
$$ ikae^{ik(0)} -ikbe^{-ik(0)} = \kappa ce^{\kappa(0)} - \kappa de^{-\kappa(0)}$$
$$ \boxed{ik(a - b) = \kappa (c - d)}$$
By inspection,
$$ \boxed{c = \frac{a+b}{2} + i\frac{k(a-b)}{2\kappa}\\d = \frac{a+b}{2} - i\frac{k(a-b)}{2\kappa}}$$
Ensuring continuity at boundary $x=a$:
$$ \Psi_{0<x<a} = \Psi_{x>a}$$
$$ \boxed{ce^{\kappa a} + d e^{-\kappa a} = fe^{ika} + ge^{-ika}} $$
and
$$ \frac{d}{dx}\Psi_{0<x<a} = \frac{d}{dx}\Psi_{x>a}$$
$$ \frac{d}{dx}(ce^{\kappa x} + de^{-\kappa x}) = \frac{d}{dx}(fe^{ikx} + ge^{-ikx})$$
$$ \kappa ce^{\kappa x} -\kappa de^{-\kappa x} = ikfe^{ikx} - ikge^{-ikx}$$
$$ \boxed{-i\frac{\kappa}{k} ce^{\kappa a} + i\frac{\kappa}{k} de^{-\kappa a} = fe^{ika} -ge^{-ika}} $$
By inspection,
$$ \boxed{f = \frac{e^{-ika}}{2}\left[ce^{\kappa a}\left(1 -i\frac{\kappa}{k}\right) + de^{-\kappa a}\left(1 + i\frac{\kappa}{k}\right)\right] \\g = \frac{e^{ika}}{2}\left[ce^{\kappa a}\left(1 + i\frac{\kappa}{k}\right)
+ d e^{-\kappa a}\left(1 - i\frac{\kappa}{k}\right)\right]}$$

Thus, there exist $a,b$ coefficients of $\Psi_{x<0}$ for which there remains non-vanishing $f,g$ coefficients of $\Psi_{x>a}$.
Two-particle wavefunction over 1D barrier:
$$\boxed{\Psi_2(x,y) = e^{px} e^{qy} - e^{qx} e^{py}}$$
is a solution to the constraint equation with barrier potential $V_0$
$$-\frac{1}{2}\left(\frac{d^2}{dx^2}\Psi + \frac{d^2}{dy^2}\Psi \right) + V_0 \Psi = E\Psi$$
since
$$ -\left(\frac{d^2}{dx^2}(e^{px} e^{qy} - e^{qx} e^{py}) + \frac{d^2}{dy^2}(e^{px} e^{qy} - e^{qx} e^{py}) \right) + V_0 (e^{px} e^{qy} - e^{qx} e^{py})$$
$$ =-\left(\frac{d}{dx}(pe^{px} e^{qy} - qe^{qx} e^{py}) + \frac{d}{dy}(qe^{px} e^{qy} - p e^{qx} e^{py}) \right) + V_0 (e^{px} e^{qy} - e^{qx} e^{py})$$
$$ =-\left(p^2e^{px} e^{qy} - q^2e^{qx} e^{py} + q^2e^{px} e^{qy} - p^2 e^{qx} e^{py} \right) + V_0 (e^{px} e^{qy} - e^{-qx} e^{-py})$$
$$=-(p^2+q^2)(e^{px} e^{qy} - e^{qx} e^{py}) + V_0 (e^{px} e^{qy} - e^{qx} e^{py}) $$
$$ (-p^2-q^2+V_0)\Psi = E\Psi$$
allows for tunnelling solution $E = -p^2-q^2+V_0$.
Similarly,
$$\boxed{\Psi_1(x,y) = e^{ipx} e^{iqy} - e^{iqx} e^{ipy} }$$
is a solution to the constraint equation outside the barrier
$$-\left(\frac{d^2}{d{x}^2}\Psi + \frac{d^2}{d{y}^2}\Psi \right) + (0)\Psi = E\Psi$$
since
$$ -\left(\frac{d^2}{d{x}^2}(e^{ipx} e^{iqy} - e^{iqx} e^{ipy}) + \frac{d^2}{d{y}^2}(e^{ipx} e^{iqy} - e^{iqx} e^{ipy}) \right)$$
$$ (p^2 + q^2)\Psi = E\Psi$$
Thus, there exist suitable $p,q$ that will satisfy the $E$ requirement throughout for tunneling through a barrier of any large $V_0$.

The wave function of entangled multi-particles cannot be represented as a product of the wave functions of each particle.

Eg.
Using 4-component ground state wavefunction of hydrogen:
$$ \Psi_1 = \begin{bmatrix} e^{-r} \\ 0 \\ i\cos\theta e^{-r} \\ i\sin\theta e^{i\phi}e^{-r} \end{bmatrix}, \Psi_2 = \begin{bmatrix} 0 \\ e^{-r} \\ i\sin\theta e^{i\phi}e^{-r} \\ -i\cos\theta e^{-r} \end{bmatrix}$$
Define
$$\Psi_{slater}(a, b) = \Psi_1(a)\Psi_2(b) - \Psi_2(a)\Psi_1(b)$$
and any two coordinates
$$ q_1 = [r_1,\theta_1,\phi_1]$$
$$ q_2 = [r_2,\theta_2,\phi_2]$$
Then
$$\Psi_{slater} (q_1, q_2) = - \Psi_{slater} (q_2, q_1)$$
and
$$\Psi_{slater} (q_1, q_1) = 0$$
$$\Psi_{slater} (q_2, q_2) = 0$$
Single slit screen probability:

Double slit controlled initial phase:

Double slit uncontrolled (decoherent) initial phase 30 trials:

Measurement decoherence at slit 1 only:

Decoherence in 1 slit is enough to wash out the entire interference.