## 4 The Replicated Action for Superfluid Bosons
According to Equation (55), the joint thermal and disorder average is obtained through  
\[\left\langle\,\mathcal{O}\,\right\rangle_{\text{th-dis}}=-\lim_{R\to 0}\frac{ \delta}{\delta f}\frac{1}{R}\mathcal{Z}^{R}[J]=\lim_{R\to 0}\frac{1}{R}\left\langle\, \mathcal{O}[\psi_{\alpha},\psi_{\alpha}^{*}]\,\right\rangle_{\text{rep}}, \tag{68}\]  
where \(\left\langle\bullet\right\rangle_{\text{rep}}\) has to be intended as a thermal average over the replicated partition function in Equation (65) or, equivalently, over the action given by the sum of Equations (66) and (67), i.e.,  
\[\left\langle\,\mathcal{O}\,\right\rangle_{\text{rep}}=\frac{\int\mathcal{D}[ \pi,\varphi]\mathcal{O}\exp\{\hbar^{-1}(S_{g}^{\text{(pure)}}[\pi,\varphi]+S_{ g}^{\text{(dis)}}[\pi,\varphi])\}}{\int\mathcal{D}[\pi,\varphi]\exp\{\hbar^{-1}(S_{g}^{ \text{(pure)}}[\pi,\varphi]+S_{g}^{\text{(dis)}}[\pi,\varphi])\}}. \tag{69}\]  
It is useful to expand the fields in Fourier series such that Equation (66) reads  
\[\langle\mathcal{Z}^{R}\rangle_{\rm dis} =\int\mathcal{D}[\pi,\varphi]\exp\bigg{\{}-\frac{1}{\hbar}\sum_{ \alpha=1}^{R}S_{g}^{\rm(pure)}[\pi_{\alpha},\varphi_{\alpha}]\bigg{\}}\ \times \tag{62}\]0  
where we have defined the vector \(\chi_{\alpha}\) as \(\chi_{\alpha}(\mathbf{q},\omega_{n})^{T}\equiv\left(\varphi_{\alpha}(\mathbf{ q},\omega_{n});\pi_{\alpha}(\mathbf{q},\omega_{n})\right)\) The pure inverse propagator in Equation (70) is given by  
\[\mathcal{M}^{-1}(\mathbf{q},\omega_{n})= \begin{pmatrix}n_{0}\frac{\hbar^{2}q^{ 2}}{m}&-\hbar\omega_{n}\\ +\hbar\omega_{n}&\frac{\hbar^{2}q^{2}}{4 mn_{0}}+\tilde{V}(\mathbf{q}) \end{pmatrix} . \tag{71}\]  
On the other hand, for Equation (63), we have now  
\[S_{g}^{\text{(dis)}}[\pi_{\alpha},\pi_{\beta},\varphi_{\alpha},\varphi_{\beta }]=-\frac{1}{2\beta\hbar L^{d}}\sum_{\mathbf{q}}\sum_{n}\gamma^{2}(\mathbf{q },\omega_{n})\pi_{\alpha}^{*}(\mathbf{q},\omega_{n})\pi_{\beta}(\mathbf{q}, \omega_{n}). \tag{72}\]  
The equation above can be decomposed in its diagonal and off-diagonal contributions, leading us to the final equation for the Euclidean action of Gaussian fluctuations. By renaming \(Q=(\mathbf{q},\omega_{n})\) and using the fact that \(x(-Q)=x^{*}(Q)\), we have  
\[S_{g}^{\text{(rep)}}[\pi,\varphi]=\frac{1}{2\beta\hbar L^{d}}\sum_{Q}\mathsf{ X}^{\dagger}(Q)\mathbb{M}^{-1}(Q)\mathsf{X}(Q). \tag{73}\]  
In the equation above, the sum over the replicas is encoded in the \((R\times 2)\)-dimensional column vector  
\[\mathsf{X}(Q)= \begin{pmatrix}\chi_{1}(Q)\\ \chi_{2}(Q)\\ \vdots\\ \chi_{R}(Q)\end{pmatrix} . \tag{74}\]  
The \(R\times 2\) matrix \(\mathbb{M}^{-1}(Q)\) in Equation (73) is built as  
\[\mathbb{M}^{-1}(Q)= \begin{pmatrix}\tilde{\mathcal{M}}^{-1}(Q)&\mathbb{B}( \gamma)&\ldots&\ldots&\mathbb{B}(\gamma)\\ \mathbb{B}(\gamma)&\tilde{\mathcal{M}}^{-1}(Q)&\ldots&\ldots&\mathbb{B}(\gamma )\\ \vdots&\vdots&\ddots&&\vdots\\ \vdots&\vdots&&\ddots&\mathbb{B}(\gamma)\\ \mathbb{B}(\gamma)&\mathbb{B}(\gamma)&\ldots&\mathbb{B}(\gamma)&\tilde{ \mathcal{M}}^{-1}(Q)\end{pmatrix}  \tag{75}\]with the generating blocks  
\[\bar{\mathcal{M}}^{-1}(Q)= \begin{pmatrix}n_{0}\frac{\hbar^{2}q^{2}}{m}&-\hbar \omega_{n}\\ +\hbar\omega_{n}&\frac{\hbar^{2}q^{2}}{4 mn_{0}}+\bar{V}(\mathbf{q})-\gamma^{2}( \mathbf{q},\omega_{n}),\end{pmatrix}  \tag{76}\]  
\[\mathbb{B}(\gamma)= \begin{pmatrix}0&0\\ 0&-\gamma^{2}(\mathbf{q},\omega_{n})\end{pmatrix} . \tag{77}\]