## Appendix B Derivation of Finite-Size Secret Key Rate  
We assume that the detection methods of [PERSON] and [PERSON] are the same as the asymptotic case and the photon subtraction operation is performed on [PERSON]'s side. In this scenario, \(N\) is the number of entire transmitted signals and \(n\) is the number of signals used for secret key extraction. The rest of the signals \(s\equiv N-n\) are utilized for estimating the parameters of insecure channels that can be implemented by sampling of \(s\) couples of correlated variables \((x_{i},y_{i})_{i=1,2,\cdots,s}\). The finite-size secret key rate can be expressed as  
\[K_{finit}=\frac{n}{N}P_{(m)}^{\mathbb{f}_{1}}[\beta I(A:B)-S_{ \epsilon_{PE}}-\Delta(n)], \tag{12}\]  
where \(\beta\), \(I(A:B)\) and \(P_{(m)}^{\mathbb{f}_{1}}\) are defined the same as above, and \(\epsilon_{PE}\) is the failure probability of parameter estimation and \(\Delta(n)\) is in connection with privacy amplification that is given by  
\[\Delta(n)=(2\dim\mathcal{H}_{B}+3)\sqrt{\frac{\log_{2}(2/\overline{ \epsilon})}{n}}+\frac{2}{n}\log_{2}(1/\epsilon_{PB}). \tag{13}\]  
Here, the parameter \(\epsilon_{PB}\) denotes the failure probability of privacy amplification and \(\overline{\epsilon}\) demonstrates the smoothing parameter. The symbol \(\mathcal{H}_{B}\) represents the Hilbert space corresponding to [PERSON]'s raw key which usually employs binary encoding so we set \(\mathcal{H}_{B}=2\). Because of the unknown channel, we have to estimate some arguments to derive \(S_{\epsilon_{PE}}\) which would minimize the secret key rate. The \(s\) couples of samples can be related as  
\[y=tx+z, \tag{14}\]  
where \(t=\sqrt{T}\in\mathbb{R}\) and \(z\) follow a centered normal distribution with variance \(\sigma^{2}=1+T(\epsilon-3)\). Moreover, we can find the covariance matrix \(\Gamma_{\epsilon_{PE}}\) that minimizes the secret key rate with a probability of at least 1-\(\epsilon_{PE}\) that is given by  
\[\Gamma_{\epsilon_{PE}}=\begin{pmatrix}X^{\prime}\mathbb{I}&t_{ min}Z_{d}^{\prime}\sigma_{z}\\ t_{min}Z_{d}^{\prime}\sigma_{z}&(t_{min}^{2}Y^{\prime}+\sigma_{max}^{2}) \mathbb{I}\end{pmatrix}, \tag{15}\]  
where \(t_{min}\) and \(\sigma_{max}^{2}\) are the minimized \(t\) and maximized \(\sigma^{2}\) of the samples, respectively. For this normal linear model, the maximum-likelihood estimators \(\hat{t}\) and \(\hat{\sigma^{2}}\) can be expressed as  
\[\hat{t}=\frac{\sum_{i=1}^{s}x_{i}y_{i}}{\sum_{i=1}^{s}x_{i}^{2}},\qquad\hat{\sigma}^{2}=\frac{1}{s}\sum_{i=1}^{s}(y_{i}-\hat{t}x_{i})^{2}. \tag{16}\]Furthermore, \(\hat{t}\) and \(\hat{\delta}^{2}\) follow the distribution that  
\[\hat{t}\sim\mathcal{N}(t,\frac{\sigma^{2}}{\sum_{i=1}^{s}x_{i}^{2}}),\qquad\frac {s\hat{\sigma}^{2}}{\sigma^{2}}\sim\chi^{2}(s-1), \tag{10}\]  
where \(t\) and \(\sigma^{2}\) are the authentic values of the parameters. Therefore, \(t_{min}\) and \(\sigma_{max}\) can be calculated as  
\[t_{min} =\sqrt{T}-z_{e_{PE}/2}\sqrt{\frac{1+T(\varepsilon-3)}{sX^{\prime }}}, \tag{11}\] \[\sigma_{max}^{2} =1+T(\varepsilon-3)+z_{e_{PE}/2}\frac{\sqrt{2}[1-T(\varepsilon-3 )]}{\sqrt{s}},\]  
where \(z_{e_{PE}/2}\) satisfies \(1-\text{erf}(z_{e_{PE}/2}/\sqrt{2})=\epsilon_{PE}/2\) and erf is the error function given by  
\[\text{erf}(x)=\frac{2}{\pi}\int_{0}^{x}e^{-t^{2}}dt. \tag{12}\]  
The value of the above-mentioned error probabilities can be assigned as  
\[\overline{\varepsilon}=\epsilon_{PE}=\epsilon_{PB}=10^{-10}. \tag{13}\]  
Finally, the secret key rate of finite-size scenario can be calculated using the derived bound \(t_{min}\) and \(\sigma_{max}\).