\section*{Funding}
Russian Science Foundation (RSF), project №19-72-10069 (experimental part), Ministry of Science and Higher Education of the Russian Federation, program №FFNN-2022-0016 for the Valiev Institute of Physics and Technology, Russian Academy of Sciences and the Foundation for the Advancement of Theoretical Physics and Mathematics BASIS (project nos. 20-1- 1-34-1, 21-1-3-40-1) (theoretical part).

\section*{Acknowledgments}
This research was performed according to the Development program of the Interdisciplinary Scientific and Educational School of Lomonosov Moscow State University "Photonic and Quantum technologies. Digital medicine". N. N. Skryabin and I. V. Dyakonov acknowledge support from the Russian Roadmap for Quantum computing (Contract No. 868-1.3-15/15-2021 dated October 5, 2021 and Contract No. P2154 dated November 24, 2021)


\section*{Supplemental document}
See Supplement for supporting content. 


\bibliographystyle{unsrt} 

\section{Derivation of the first-order autocorrelation function of  quasi-thermal light}

Let us derive the correlation function of the quasi-thermal field obtained using a disk rotating with an angular velocity $\Omega$ (Fig.~1 in the main text). To describe the field scattered at an angle to the disk surface normal, we follow the paper~\cite{Crosignani1971}. While in that work a plane wave fallen on the disk was bounded by a rectangular window, for our case we consider Gaussian beam, which amplitude is


\begin{equation}
    E(x,y)=E_0 \exp \left(-\frac{(x-L)^2+y^2}{2 w^2}\right).
\end{equation}
Here,  $w$ is a waist of the incident laser beam, which is focusing on the disk surface, $E_0$ corresponds to the field in the center of the beam waist. Also for simplicity we consider the case when the center of the beam waist is at point $\{x=L, y=0, z=0\}$.

Taking into account the intermediate formulas in \cite{Crosignani1971} and equation (13), the dependence of the first-order correlation function on time delay can be written as
\begin{multline}
        g^{(1)}(\tau) \varpropto \exp[\text{i}\omega \tau] \int\limits_{-\infty}^{+\infty}dx''dy''\exp\left[-\text{i}k_y y''(\cos\left(\Omega\tau\right)-1) \right]\exp\left[\text{i}k_y x''\sin\left(\Omega\tau\right)\right]\times\\
    \times E(x'',y''-\Omega\tau L)E(x'', y''),
\end{multline}
where $\omega$ is the frequency of the laser radiation.
After integrating and using series expansion of the trigonometric functions, one has:

\begin{equation}\label{eq:2_3}
    g^{(1)}(\tau) \varpropto \exp\left[-\cfrac{L^2\Omega^2\tau^2 }{4 w^2}\right]\exp \left[\text{i}(\omega+k_y L \Omega)\tau \right] =\exp\left[-\cfrac{V^2\tau^2 }{4 w^2}\right]\exp \left[\text{i}(\omega+\Delta)\tau \right] .
\end{equation}
Here, $V$ is speed of the disk point at the center of the waist, the first factor is the standard Gaussian envelope corresponding to the correlation properties of the thermal state, with the correlation time $w/V$, which corresponds to the time of passage of the disk point through the beam waist, and the second factor is sinusoidal oscillation, where $\omega$ is the frequency of the radiation incident on the disk, and $\Delta = k_y V$ corresponds to the Doppler shift associated with a non-zero projection of the disk velocity on the direction of light propagation.


\section{Intensity correlation function of two interfering quasi-thermal fields}


The Gaussian (thermal) distribution is remarkable by its higher-order moments can be expressed in terms of first-order moments.

Let's put the intensities at inputs $1$ and $j$ as $I_1$ and $I_j$, and intensities at outputs $1$ and $i$ as $I_1^{\text{(out)}}$ and $I_i^{\text{(out)}}$. In these terms the output second-order cross-correlation function equals:

\begin{multline}\label{eq:124}
    \left\langle {I_1^{\text{\text{(out)}}}I_i^{\text{\text{(out)}}}(\tau )} \right\rangle  = \left\langle {E_1^{\text{(out)}*} E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}{E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle  = \\
    = {\left\langle {E_1^{\text{(out)}*} {E_1^{\text{\text{(out)}}}}} \right\rangle } {\left\langle {E_i^{\text{(out)}*} (\tau ){E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle } + \left\langle {E_1^{\text{(out)}*} {E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle \left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}} \right\rangle \\ = \left\langle I_1^{\text{\text{(out)}}} \right\rangle \left\langle I_i^{\text{\text{(out)}}} \right\rangle +\left\langle {E_1^{\text{(out)}*} {E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle \left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}} \right\rangle
\end{multline}

We are interested in the last term. Let's describe it in more detail:
\begin{multline}
    \left\langle {E_1^{\text{(out)}*} {E_i^{\text{(out)}}}(\tau )} \right\rangle  = {\left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{(out)}}}} \right\rangle ^ * } = \\ = \left\langle {\left( {M_{11}E_1^*  + M_{1j}E_j^* } \right)\left( {M_{i1}{E_1}(\tau ) + M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}{E_j}(\tau )} \right)} \right\rangle  = \\
 = M_{11}M_{i1}\underbrace {\left\langle {E_1^* {E_1}(\tau )} \right\rangle }_{\left\langle {I_1} \right\rangle g_1^{(1)}(\tau )} + M_{11}M_{ij}{{\mathop{\rm e}\nolimits} ^{\text{i}\varphi }}\underbrace {\left\langle {E_1^* {E_j}(\tau )} \right\rangle }_0 +\\+ M_{1j}M_{i1}\underbrace {\left\langle {E_j^* {E_1}(\tau )} \right\rangle }_0 +M_{1j}M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}\underbrace {\left\langle {E_j^* {E_j}(\tau )} \right\rangle }_{\left\langle {I_j} \right\rangle g_j^{(1)}(\tau )},
\end{multline}

\begin{multline}\label{eq:2_4}
    { {\left\langle {E_1^{\text{(out)}*} {E_i^{\text{(out)}}}(\tau )} \right\rangle } ^2} = { {M_{11}M_{i1}\left\langle {I_1} \right\rangle g_1^{(1)}(\tau ) + M_{1j}M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}\left\langle {I_j} \right\rangle g_j^{(1)}(\tau )}^2} = \\
 \left\langle {I_1} \right\rangle \left\langle {I_j} \right\rangle M_{11}  M_{1j}  M_{ij} M_i1e^{-\text{i} \varphi } g_1^{(1)}(\tau ) g_j^{(1)}(\tau )^*+
 \left\langle {I_1} \right\rangle \left\langle {I_j} \right\rangle M_{11}  M_{1j}  M_{ij} M_i1 e^{\text{i} \varphi } g_j^{(1)}(\tau ) g_1^{(1)}(\tau )^*+\\
 \left\langle {I_1} \right\rangle^2 M_{11}^2 M_{i1}^2 g_1^{(1)}(\tau ) g_1^{(1)*}(\tau )+\left\langle {I_j} \right\rangle^2 M_{1j} ^2 M_{ij}^2 g_j^{(1)}(\tau ) g_j^{(1)*}(\tau )
 .
\end{multline}
 	
Now we replace the first-order correlation functions here with \eqref{eq:2_3}. Let them differ in Doppler shift for fields at different inputs. This can be achieved by collecting radiation scattered on the same RGGD at different angles. Then assuming
 	
\begin{gather}
    g_1^{(1)}(\tau ) = \exp \left(  - \frac{V^2\tau ^2}{4w^2} \right)\exp \left[ {  \text{i}\left( {{\omega} + {\Delta _1}} \right)\tau } \right],  \\
    g_j^{(1)}(\tau ) = \exp \left(  - \frac{V^2\tau ^2}{4w^2} \right)\exp \left[ {  \text{i}\left( {{\omega} + {\Delta _j}} \right)\tau } \right],\\
\text{and } d = {\Delta _j} - {\Delta _1}
\end{gather}
for the normalized cross-correlation function we have:
\begin{multline}\label{eq:g2t}
    g^{(2)}(\tau) = 1 + \exp{\left(-\frac{\tau ^2 V^2}{2 w^2}\right)} \frac{1}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)} \times \\
    \left[\left\langle I_1 \right\rangle^2 M_{11}^2 M_{i1}^2 +\left\langle I_j \right\rangle^2 M_{1j}^2 M_{ij}^2 +2 \left\langle I_1 \right\rangle \left\langle I_j \right\rangle M_{11} M_{i1} M_{1j} M_{ij} \cos \left(d \tau +\varphi \right)\right] .
\end{multline}
In the special case of the ideal symmetric beam splitter ($M_{11} = M_{12} = M_{21} = M_{22} = 1/\sqrt{2}$ and the only MEP equals $\pi$) \eqref{eq:g2t} is reduced to Eq.~(4) of the main article:
\begin{equation*}\label{eq:g2_ideal_BS}
       g^{(2)}(\tau) = 1+ \exp\left(-\frac{\tau ^2}{2 \sigma^2}\right)\times
       \left[\frac{ \left< I_1 \right>^2+\left< I_2 \right>^2}{\left(\left< I_1 \right>+\left< I_2 \right>\right)^2}-\frac{2 \left< I_1 \right> \left< I_2 \right> }{\left[ \left< I_1 \right>+\left< I_2 \right>\right] ^2}\cos \left(d \tau \right)\right].
\end{equation*}

In general case, \eqref{eq:g2t} can be represented as Eq.~(5) from the main article:
\begin{equation*}\label{eq:g2_approx}
    g^{(2)}\left( \tau  \right) = C + \exp \left( { - {\frac{\tau ^2}{2\sigma^2}}} \right)\left[ {A + B\cos \left( {d \tau  + \varphi_{\text{thermal}} } \right)} \right],
\end{equation*}



where parameters:

\begin{equation}
\begin{aligned}
A &= \frac{\left\langle I_1 \right\rangle^2 M_{11}^2 M_{i1}^2 +\left\langle I_j \right\rangle^2 M_{1j}^2 M_{ij}^2}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)};\\
B &= \frac{2 \left\langle I_1 \right\rangle \left\langle I_j \right\rangle M_{11} M_{1j} M_{ij} M_{i1}}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)};\\
C &=1;\\
\sigma &= \frac{w}{V}.
\end{aligned}
\end{equation}



\section{Implementation of the beam splitter with adjustable matrix element phase}

The arrangement of the polarization elements (Fig.~2 in the main text) in the order PBS, HWP at -22.5$^\circ$, tilted Glan-Thompson prism, and quarter-wave plate QWP at 45$^\circ$ and second PBS allows us to change the MEP of the transfer matrix from $0$ to $\pi$. Swapping the outputs of the circuit, we get the same matrix element phase (MEP) with the opposite sign. Thus, we can change the MEP from $-\pi$ to $\pi$. The  Fig.~\ref{fig:phase_eta} shows the dependence of the MEP and minimum intensity transmission coefficient $T_{min}$ on the prism transmission coefficient of horizontally polarized beam $\eta$. It can be seen that the MEP varies from $0$ to $\pi$, without zeroing out the transmission anywhere.

When calculating this dependence, we used Jones calculus:
\begin{equation}
   \vec{H}=\left(
\begin{array}{c}
 1 \\
 0 \\
\end{array}
\right),
\Vec{V}=\left(
\begin{array}{c}
 0 \\
 1 \\
\end{array}
\right)
\end{equation}

Glan-Thompson prism matrix can be represented as

\begin{equation}
    \hat{P}=\left(
\begin{array}{cc}
 \sqrt{\eta} & 0\\ 
 0 & 1
\end{array}
\right),
\end{equation}
where $\eta$ is transmission coefficient for muted polarization, which can be changed by tilting the prism around the vertical axis. In our case, in the initial position of the prism, when $\eta=0$, the horizontal polarization is not transmitted.

Matrices of wave plates are represented as
\begin{equation}
    \hat{HWP} = \left(
\begin{array}{cc}
 \cos ^2(\alpha)-\sin ^2(\alpha) & 2 \sin (\alpha) \cos (\alpha) \\
 2 \sin (\alpha) \cos (\alpha) & \sin ^2(\alpha)-\cos ^2(\alpha) \\
\end{array}
\right),
\end{equation}


\begin{equation}
    \hat{QWP} = \left(
\begin{array}{cc}
 \cos ^2(\beta)+\text{i} \sin ^2(\beta) & (1-\text{i}) \sin (\beta) \cos (\beta) \\
 (1-\text{i}) \sin (\beta) \cos (\beta) & \sin ^2(\beta)+\text{i} \cos ^2(\beta) \\
\end{array}
\right),
\end{equation}
where $\alpha$ and $\beta$ are angles of fast axes.

The resulting matrix $\hat{U}$ will be

\begin{equation}
    \hat{U} = \hat{HWP} \times \hat{P} \times \hat{QWP}.
\end{equation}
Then we can find an equivalent matrix with zeroed first row and column and determine the MEP of the element $\{2,2\}$ and the minimal transmission coefficient $T_{min}$.






\begin{figure}[htbp]
\centering\includegraphics[width=0.7\linewidth]{Supplement_phase_realization.pdf}
\caption{MEP value $\varphi/\pi$ and minimum intensity transmission coefficient  $T_{min}$ vs Glan prism transmission coefficient $\eta$ for horizontal polarization, which can be tuned by rotation of the prism. }\label{fig:phase_eta}
\end{figure}



\bibliographystyle{unsrt} 
