Consider a particle constrained to a ring of circumference $L$. Following this paper, the position eigenstate on a circle can be expressed in terms of the position eigenstates on the real line as
$$ \left|x\right\rangle _{\mathcal{S}^{1}}=\sqrt{\frac{L}{2\pi\hbar}}\sum_{n=-\infty}^{\infty}\left|x+nL\right\rangle \tag{*} $$
where $x \in \left[0, L\right)$. Here $\left| ~ \right\rangle _{\mathcal{S}^{1}}$ represents a state on the circle, while $\left| ~ \right\rangle$ without the subscript denotes the usual state on the real line $\mathbb{R}$. This relation can be derived using the Dirac comb and the momentum-space resolution of the identity.
The density matrix in the position basis on the circle is
$$ _{\mathcal{S}^{1}}\left\langle x\right|e^{-\beta\frac{\hat{p}^{2}}{2m}}\left|y\right\rangle _{\mathcal{S}^{1}}. $$
Using the completeness relation for the (quantized) momenta $p_n = 2\pi\hbar n / L$ it can be proven that this expression reduces to
$$ \frac{1}{L} \sum_{n=-\infty}^{\infty} e^{-\beta \frac{p_n^2}{2m} } e^{i p_n (x - y) / \hbar}, $$
which is the standard expression for the density matrix of a particle on a ring.
An alternative derivation uses the expansion $(*)$ to replace, say, the ket $\left|y\right\rangle _{\mathcal{S}^{1}}$. In that case, the resulting Dirac comb inside the momentum integral picks out the discrete allowed momenta, leading to the same result:
$$ \begin{align*} _{\mathcal{S}^{1}}\left\langle x\right|e^{-\beta\frac{\hat{p}^{2}}{2m}}\left|y\right\rangle _{\mathcal{S}^{1}} & =\sqrt{\frac{L}{2\pi\hbar}}\sum_{q=-\infty}^{\infty}{}_{\mathcal{S}^{1}}\left\langle x\right|e^{-\beta\frac{\hat{p}^{2}}{2m}}\left|y+qL\right\rangle \\ & =\sqrt{\frac{L}{2\pi\hbar}}\sum_{q=-\infty}^{\infty}\int_{-\infty}^{\infty}dp{}_{\mathcal{S}^{1}}\left\langle x\right|e^{-\beta\frac{\hat{p}^{2}}{2m}}\left|p\right\rangle \left\langle p\right|\left.y+qL\right\rangle \\ & =\sqrt{\frac{L}{2\pi\hbar}}\sum_{q=-\infty}^{\infty}\int_{-\infty}^{\infty}dpe^{-\beta\frac{p^{2}}{2m}}{}_{\mathcal{S}^{1}}\left\langle x\right|\left.p\right\rangle \frac{1}{\sqrt{2\pi\hbar}}e^{-ip(y+qL)/\hbar}\\ & =\sqrt{\frac{L}{2\pi\hbar}}\sum_{q=-\infty}^{\infty}\int_{-\infty}^{\infty}dpe^{-\beta\frac{p^{2}}{2m}}\frac{1}{\sqrt{L}}e^{ipx/\hbar}\frac{1}{\sqrt{2\pi\hbar}}e^{-ip(y+qL)/\hbar}\\ & =\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}dpe^{-\beta\frac{p^{2}}{2m}}e^{ip(x-y)/\hbar}\sum_{q=-\infty}^{\infty}e^{-ipqL/\hbar}\\ & =\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}dpe^{-\beta\frac{p^{2}}{2m}}e^{ip(x-y)/\hbar}2\pi\sum_{k=-\infty}^{\infty}\delta\left(\frac{p}{\hbar}L-2\pi k\right)\\ & =\frac{1}{L}\int_{-\infty}^{\infty}dpe^{-\beta\frac{p^{2}}{2m}}e^{ip(x-y)/\hbar}\sum_{k=-\infty}^{\infty}\delta\left(p-p_{k}\right)\\ & =\frac{1}{L}\sum_{k=-\infty}^{\infty}e^{-\beta\frac{p_{k}^{2}}{2m}}e^{ip_{k}(x-y)/\hbar} \end{align*} $$
However, if one tries to apply $(*)$ to the bra $_{\mathcal{S}^{1}}\left\langle x\right|$ as well, the calculation yields a product of two Dirac combs inside the same momentum integral, which seems ill-defined and makes the derivation break down. The aforementioned paper does not address this difficulty. Instead, they suggest calculating the density matrix using $\sum_{q=-\infty}^{\infty}\left\langle x\right|e^{-\beta\frac{\hat{p}^{2}}{2m}}\left|y+qL\right\rangle$, i.e., by fixing the bra to the fundamental cell (neglecting its periodic images) and allowing the ket to sum over different periodic images (using $(*)$). This indeed leads to the correct expression, however this approach seems to be an ad hoc one, because they do not explain why the bra has to be fixed.
Question: What is the fundamental reason why the expression $(*)$ should be applied only to one side (bra or ket, but no both) in this context? And how is it justified that using it only once suffices to obtain the correct density matrix on the ring?
