\section{Supplementary material}
\subsection{Differential equations}
\label{app:diff}

In this section, we derive the equations of motion for the expectation values of the operators. As an example, we consider the expectation value of $\hat{\tau}_k$ in the Schr\"{o}dinger's picture:
\begin{equation}
 \tau_k(t)= \langle \psi(t)| \hat \tau_k |\psi(t)\rangle,
\end{equation}
where $|\psi(t)\rangle$ denotes the wave function at time $t$. The Schr\"odinger's equation gives:
\begin{equation}
i \frac{d |\psi(t)\rangle}{dt} = \hat H_{eff}(t) |\psi(t)\rangle.
\end{equation}
The equation of motion for $\tau_k(t)$ is then:
\begin{equation}
-i \frac{d \tau_k(t)}{dt}= \langle \psi(t)| \left[\hat H_{eff}(t), \hat \tau_k \right] |\psi(t)\rangle.
\end{equation}
The commutator on the right-hand side consists of six possible operators: $\{\hat \nu_{k\uparrow},\hat \nu_{k\downarrow},\hat n_{k\uparrow},\hat n_{k\downarrow},\hat \tau_k,\hat \sigma_k \}$. By calculating the derivative of their expectation values to $|\psi(t)\rangle$, we get a set of closed-equations for the expectation values of the operators:
\begin{equation}\label{diffequationgroup}
 \begin{split}
 -i \frac{d \nu_{k \uparrow} }{dt} = & 2\epsilon_{k\uparrow} \nu_{k\uparrow} + \Delta^* \sigma_{-k} - \Delta^* \sigma_k + \alpha_k \tau_{-k} + \alpha_k\tau_k \\
-i  \frac{d \nu_{k \downarrow} }{dt} = & 2\epsilon_{k\downarrow} \nu_{k\downarrow} + \Delta^* \sigma^*_{k} - \Delta^* \sigma^*_{-k} +  \alpha_k \tau_{-k} + \alpha_k\tau_k \\
-i \frac{d n_{k\uparrow}}{dt} = & - \Delta^* \tau_k^* + \Delta \tau_k -\alpha_k \sigma_k^* -\alpha_k \sigma_k \\
-i \frac{d n_{k\downarrow}}{dt} = & -\Delta^* \tau_{-k}^* + \Delta \tau_{-k} + \alpha_k \sigma^*_k + \alpha_k \sigma_k \\
-i \frac{d\tau_k}{dt} = & (\epsilon_{k\uparrow}+\epsilon_{k\downarrow}) \tau_k + \Delta^* n_{-k\downarrow} - \alpha_k \nu_{k\downarrow} + \Delta^* n_{k\uparrow} \\ & - \alpha_k \nu_{k\uparrow} - \Delta^* \\
-i \frac{d\sigma_k}{dt} = & (\epsilon_{k\uparrow}-\epsilon_{k\downarrow}) \sigma_k - \Delta^* \nu^*_{k\downarrow} - \Delta \nu_{k\uparrow} - \alpha_k n_{k\downarrow} \\ & + \alpha_k n_{k\uparrow}.
 \end{split}
\end{equation}
The pairing gap $\Delta(t)$ in Eq.~(\ref{diffequationgroup}) can be expressed as
\begin{equation}\label{defdeltaauxiliary}
\Delta(t) = U_f\sum_{k} \tau_k^{\ast}(t),
\end{equation}
which can be solved self-consistently from Eq.~(\ref{diffequationgroup}) and Eq.~(\ref{defdeltaauxiliary}).

\subsection{Calculation of tunneling conductances}
\label{app:method}

To calculate the tunneling conductance, we employ the numerical operator method first introduced in Ref.~\cite{pei12} for the real time Kondo problem, which was then extended to time-dependent systems~\cite{pei13a} and topological superconductors~\cite{pei13b}.

The full Hamiltonian describing the reservoir and the superfluid is expressed as
\begin{equation}
 \hat H_{tot}(t) = \hat H_{eff}(t) + \hat H_{res} + \hat H_V.
\end{equation}
We define the synonyms $\hat d_{j\sigma 1}=\hat c^\dag_{j\sigma}$ and $\hat d_{j\sigma 0} = \hat c_{j\sigma}$ and re-express the conductance as
\begin{equation}\label{defdiffcondapp}
 G(t)= \left( -2 g_V\textbf{Im}\sum_\sigma  \frac{d}{d\mu_R} \langle \hat d_{-1,\sigma, 1}(t) \hat d_{0,\sigma, 0}(t) \rangle \right)_{\mu_R=0}.
\end{equation}
We first calculate the operators $\hat d_{-1,\sigma, 1}(t)$ and $\hat d_{0,\sigma, 0}(t)$ in the Heisenberg picture defined as
\begin{equation}\label{defheisenbergoperator}
\begin{split}
 \hat d_{\alpha} (t)= & \lim_{\tau\to 0} e^{i\hat H_{tot}(t_0) \tau} \cdots e^{i\hat H_{tot} \left(t_{N-1} \right) \tau} \hat d_{\alpha} \\ & \times e^{-i\hat H_{tot} \left(t_{N-1} \right) \tau} \cdots e^{-i\hat H_{tot}(t_0) \tau},
\end{split}
 \end{equation}
where $\alpha$ is the abbreviation of $(j,\sigma, s)$ ($s=0,1$), $\hat H_{tot}$ is the time-dependent Hamiltonian, $N$ is the total number of steps, $t_0=0$ is the initial time, and the intermediate time $t_n=n\tau$, where the time step $\tau=t/N$.

The calculation of Eq.~(\ref{defheisenbergoperator}) is divided into $N$ steps with the time evolving from $t_{N-1}$ to $t_0$. In the {\it n}-th step, we undress the pair of operators $e^{i\hat H_{tot} \left(t_{N-n} \right) \tau} $ and $e^{-i\hat H_{tot} \left(t_{N-n} \right) \tau} $ from $\hat d_\alpha$. Since $\hat H_{tot}$ is quadratic, in principle, the intermediate results at each step can be written as
\begin{equation}\label{heisenbergsolution}
 \hat {\tilde d}_{\alpha} (t_n) = \sum_{\alpha'} W_{\alpha,\alpha'}(t_n) \hat d_{\alpha'},
\end{equation}
with $\hat {\tilde d}_{\alpha} (t_n) $ defined as
\begin{equation}\label{defintermediate}
\begin{split}
 \hat {\tilde d}_{\alpha} (t_n) = & e^{i\hat H_{tot}(t_{N-n}) \tau} \cdots e^{i\hat H_{tot} \left(t_{N-1} \right) \tau} \hat d_{\alpha} \\ & \times e^{-i\hat H_{tot} \left(t_{N-1} \right) \tau} \cdots e^{-i\hat H_{tot}(t_{N-n}) \tau},
\end{split}
 \end{equation}
 and $W_{\alpha,\alpha'}(t_n)$ denoting the propagator. At $n=N$, we get
\begin{equation}
\hat d_{\alpha} (t)=\sum_{\alpha'} W_{\alpha,\alpha'}(t_N) \hat d_{\alpha'}.
\end{equation}

To calculate the propagator $W_{\alpha,\alpha'}(t_N)$, we derive an iterative relation. According to Eq.~(\ref{defintermediate}), we have
\begin{equation}
\hat {\tilde d}_{\alpha} (t_n) =e^{i\hat H_{tot}(t_{N-n}) \tau} \hat {\tilde d}_{\alpha} (t_{n-1}) e^{-i\hat H_{tot}(t_{N-n}) \tau} .
\end{equation}
Substituting Eq.~(\ref{heisenbergsolution}) into the equation above, we find
\begin{equation}\label{iterativerelation}
\begin{split}
& \sum_{\alpha'} W_{\alpha,\alpha'}(t_n) \hat d_{\alpha'} \\ &
= \sum_{\alpha'} W_{\alpha,\alpha'}(t_{n-1})  \bigg( \hat d_{\alpha'}+ i\tau [\hat H_{tot}(t_{N-n}), \hat d_{\alpha'}] \\ & + \frac{(i\tau)^2}{2} [\hat H_{tot}(t_{N-n}),[\hat H_{tot}(t_{N-n}),\hat d_{\alpha'}]] + O(\tau^3) \bigg).
\end{split}
\end{equation}
The calculation of these commutators is straightforward with the results
\begin{equation}\label{commutatorexp}
[\hat H_{tot}(t_j), \hat d_{\alpha}]= \sum_{\beta} G_{\alpha,\beta}(t_j) \hat d_{\beta}.
\end{equation}
Substituting Eq.~(\ref{commutatorexp}) into Eq.~(\ref{iterativerelation}) and comparing the left and the right sides, we obtain the iterative relation
\begin{equation}\label{witerativerelation}
\begin{split}
 W_{\alpha,\alpha'}(t_n)= & W_{\alpha,\alpha'}(t_{n-1}) +i\tau \sum_{\beta} W_{\alpha,\beta}(t_{n-1}) G_{\beta,\alpha'}(t_{N-n}) \\ & -\frac{\tau^2}{2} \sum_{\beta,\beta'} W_{\alpha,\beta}(t_{n-1}) G_{\beta,\beta'}(t_{N-n}) G_{\beta',\alpha'}(t_{N-n})  .
\end{split}
\end{equation}
Here we drop terms of the order $O(\tau^3)$. In principle, the error caused by the discretization of time can be made arbitrarily small by letting $\tau\to 0$. In practice, we adaptively choose $\tau$ to make the discretization error negligible.

According to Eq.~(\ref{witerativerelation}), we start from $n=0$ when $W_{\alpha,\alpha'}(t_0)= \delta_{\alpha,\alpha'}$, and iteratively work out $W_{\alpha,\alpha'}(t_n)$ until $t_N=t$. This approach typically works well, due to the critical fact that the matrix $G_{\alpha,\beta}$ has only very few non-zero elements, such that the number of non-zero $W_{\alpha,\beta}$ which needs to be kept at each step increases slowly according to Eq.~(\ref{witerativerelation}). However, when $N$ is very large, which is inevitable since we must choose a small $\tau$ and at the same time a large $t$ to obtain the asymptotic limit of $G(t)$, there will be too many non-zero $W_{\alpha,\beta}$ which cannot be simultaneously kept. To work around this, we apply a truncation scheme and keep only a fixed number of non-zero propagators (the number is denoted by $M$) with the largest magnitudes at each step. This truncation scheme is necessary for obtaining the steady limit of $G(t)$. We decide the value of $M$ adaptively, i.e.,
set an original $M$ and
increase it until the desired precision is obtained. The logic in doing so is that the truncation error goes to zero in the limit $M\to \infty$.

Substituting the expression of $\hat d_{\alpha}(t)$ into Eq.~(\ref{defdiffcondapp}), we have
\begin{equation}
\begin{split}
 G(t)= & -2 g_V\textbf{Im}\sum_{\sigma\beta\beta'}  W_{-1\sigma 1,\beta}(t) W_{0\sigma 0,\beta'}(t) \\ & \times \left( \frac{d}{d\mu_R} \langle \hat d_{\beta} \hat d_{\beta'} \rangle \right)_{\mu_R=0}.
\end{split}
\end{equation}
Next we calculate the derivative of correlation functions. The correlation functions of the initial state are non-zero only if the two field operators are both in the superfluid or both in the reservoir, since the superfluid and the reservoir are decoupled at the initial time. Considering that the initial state of the superfluid located on sites from $0$ to $L-1$ is independent of $\mu_R$, we immediately get
\begin{equation}
 \frac{d}{d\mu_R} \langle \hat d_{j\sigma s} \hat d_{j'\sigma' s'} \rangle =  0,
\end{equation}
as $j,j' \geq 0$. As $j,j'< 0$, we notice that the reservoir is in the ground state of a free Fermi gas at the chemical potential $\mu_R$, and get
\begin{eqnarray}
\left\{ \begin{array}{c} \langle \hat d_{j\sigma 1} \hat d_{j'\sigma' 0} \rangle = \delta_{\sigma,\sigma'} \displaystyle \frac{\sin((j-j')\theta )}{\pi(j-j')} \\ \langle \hat d_{j\sigma 0} \hat d_{j'\sigma' 1} \rangle = \delta_{\sigma,\sigma'}\left( \delta_{j,j'}-\displaystyle \frac{\sin((j-j')\theta )}{\pi(j-j')} \right) \end{array} \right.
\end{eqnarray}
with $\theta=\arccos(-\mu_R/2g_l)$ and $g_l$ the hopping between two neighbor sites of the reservoir. Finally we find
\begin{eqnarray}
\begin{split}
\left(\displaystyle\frac{d}{d\mu_R} \langle \hat d_{j\sigma s} \hat d_{j'\sigma' s'} \rangle\right)_{\mu_R=0} =& \left[ (1-\delta_{s,s'})(2\delta_{s,1}-1) \right] \delta_{\sigma,\sigma'} \\ & \times \frac{\cos \left( (j-j')\displaystyle\frac{\pi}{2} \right)}{2g_l \pi}.
\end{split}
\end{eqnarray}
\end{document}
