## Appendix A Secret Key Rate Analysis of THz-CVQSS in the Finite-Size Scenario  
In the actual CVQSS system, the number of quantum signals transmitted between the dealer and the players is limited at finite-size of data, which has a profound impact on  
Figure 7: The finite-size SKR of inter-satellite THz-CVQSS versus the reconciliation efficiency with different numbers of players. \(\xi_{0}=0.001,\eta=0.6,v_{el}=0.1,V_{M}=7\), \(d=2000\) km, \(f=50\) THz and \(M=10^{10}\).  
the transmission distance and SKR. Thus, in this section, we analyze the finite-size SKR of THz-CVQSS with reverse reconciliation which can be expressed by [40]  
\[R_{f}=\frac{m}{M}\big{(}\beta I_{AB}-S_{\Theta_{PE}}(AE)-\Delta(m)\big{)}, \tag{1}\]  
where, \(m\) denotes the number of quantum signals utilized to share the secure keys. \(M\) denotes the total number of quantum signals exchanged. \(h=M-m\) denotes the remaining number of quantum signals utilized to parameter estimation. \(\beta\) and \(I_{AB}\) are defined in Equation (3). \(S_{\Theta_{PE}}(AE)\) denotes the maximal Holevo information compatible with the statistics except with probability \(\Theta_{PE}\) and can be calculated by the covariance matrix \(\Omega_{\Theta_{PE}}\). \(\Delta(m)\) denotes a security parameter which is related to the privacy amplification and can be expressed by [40]  
\[\Delta(m)=\frac{2}{m}\log_{2}\bigg{(}\frac{1}{\Theta_{PA}}\bigg{)}+7\sqrt{ \frac{\log_{2}\left(\frac{2}{8}\right)}{m}}, \tag{2}\]  
where, \(\Theta_{PA}\) denotes the probability of lack of success in privacy amplification. \(\bar{\Theta}\) denotes the smoothing parameter.  
Next, we will calculate the covariance matrix \(\Omega_{\Theta_{PE}}\). The finite-size SKR \(R_{f}\) is minimized by \(\Omega_{\Theta_{PE}}\) with a probability of at least \(1-\Theta_{PE}\). The covariance matrix \(\Omega_{\Theta_{PE}}\) can be obtained by using \(h\) couples of correlated variables \((x_{i},y_{i})_{i=1...h}\). We consider a normal model for the correlated variables given by  
\[y=qx+p, \tag{3}\]  
where, \(q=\sqrt{T_{1}}\) and \(p\) denotes a centered normal variable with the variance \(\omega^{2}=T_{1}\xi_{tot}+1\), \(\xi_{tot}=\sum_{i=1}^{n}\xi_{i}\). The covariance matrix \(\Omega_{\Theta_{PE}}\) can be calculated by  
\[\Omega_{\Theta_{PE}}=\left[\begin{array}{cc}V\text{I}_{2}&q_{\min}\sqrt{V^{ 2}-1}\sigma_{z}\\ q_{\min}\sqrt{V^{2}-1}\sigma_{z}&\big{(}q_{\min}^{2}(V-1)+\omega_{\max}^{2} \big{)}\text{I}_{2}\end{array}\right], \tag{4}\]  
where, \(q_{\min}\) is the minimal value of \(\hat{q}\) and \(\omega_{\max}^{2}\) is the maximal value of \(\hat{\omega}^{2}\). The estimated value of \(\hat{q}\) and \(\hat{\omega}^{2}\) can be obtained by the maximum likelihood estimation  
\[\hat{q} =\frac{\sum_{i=1}^{h}x_{i}y_{i}}{\sum_{i=1}^{h}x_{i}^{2}}, \tag{5}\] \[\hat{\omega}^{2} =\frac{1}{h}\sum_{i=1}^{h}(y_{i}-\hat{q}x_{i})^{2}.\]  
Through the large number theorem, \(\hat{q}\) and \(\hat{\omega}^{2}\) satisfy the following distributions  
\[\hat{q}\sim N\bigg{(}q,\frac{\omega^{2}}{\sum_{i=1}^{h}x_{i}^{2}}\bigg{)}, \quad\frac{h\hat{\omega}^{2}}{\omega^{2}}\quad\sim\chi^{2}(h-1). \tag{6}\]  
Owing to the estimated values of \(\hat{q}\) and \(\hat{\omega}^{2}\) are true values, i.e., \(E[\hat{q}]=\sqrt{T_{1}}\) and \(E\big{[}\hat{\omega}^{2}\big{]}=T_{1}\xi_{tot}+1\), we then can obtain the expressions  
\[q_{\min} \approx\sqrt{T_{1}}-z_{\delta}\sqrt{\frac{T_{1}\xi_{tot}+1}{hV_{M }}}, \tag{7}\] \[\omega_{\max}^{2} \approx z_{\delta}\frac{\sqrt{2}(T_{1}\xi_{tot}+1)}{\sqrt{h}}+1+ T_{1}\xi_{tot},\]where, \(z_{\delta}\) satisfies \(1-\mathrm{erf}\Big{(}z_{\delta/\sqrt{2}}\Big{)}/2=\Theta_{PE}/2\), here, \(\mathrm{erf}(x)\) denotes the error function given by  
\[\mathrm{erf}(x)=\int_{0}^{x}\frac{2}{\sqrt{\pi}}e^{-u^{2}}\;\mathrm{d}u. \tag{10}\]  
Note that, we set the parameter values as  
\[\bar{\Theta}=\Theta_{PE}=\Theta_{PB}=10^{-10},\quad m=h=M/2. \tag{11}\]