\subsection{Supplementary Note 1: Remarks on the convergence towards the completely mixing evolution}

As discussed in the main text, for $q=1$ we may formally study convergence
towards a fully mixing evolution via the
diamond distance $\|\mathcal E_{\rm ensemble}-\mathcal E_{\rm Haar}\|_\diamond$
between the quantum channel that mathematically describes
the average over $m$ discrete unitaries, $\mathcal E_{\rm ensemble}[\rho] =\sum_{j=1}^m U_j \rho U_j^\dagger$,
and the quantum channel that describes continuous averages over Haar-random unitaries
$\mathcal E_{\rm Haar}[\rho] =\int dU \, U \rho U^\dagger$. Such distance can be written as
(see {\it e.g.}~\cite{Watrous2018})
\begin{equation}
	\|\mathcal E_{\rm ensemble} -\mathcal E_{\rm Haar}\|_\diamond  =
	\max_\Psi \| (\openone\otimes\mathcal E_{\rm ensemble}-\openone\otimes\mathcal E_{\rm Haar})(\Psi)\|_1
	\label{diamond},
\end{equation}
where the maximization is over quantum states $\Psi=\ket\Psi\!\bra\Psi$. Although such maximization can be efficiently solved
using semidefinite programming for small dimensional Hilbert spaces, it is generally complex to estimate the
diamond distance from measurement data. Therefore, we study simpler quantities that are related to the diamond distance via bounds.

We first recall one of the possible definitions of the trace norm (see {\it e.g.}~\cite{Bengtsson2017}
Theorem 13.2)
\begin{equation}
	\|\rho-\sigma\|_1 = \max_{E} \sum_i |\Tr[E_i (\rho-\sigma)]|,
	\label{helstrom}
\end{equation}
where the maximum is taken over all possible POVMs such that $\sum_i E_i = \openone$.
With the above definition, we may introduce the following inequality
\begin{equation}
	\|\mathcal E_{\rm ensemble} -\mathcal E_{\rm Haar}\|_\diamond
	\le N \| \chi_{{\rm ensemble}}-\chi_{{\rm Haar}}\|_1~,
	\label{diamondineq}
\end{equation}
where $\chi_{t}=\openone\otimes\mathcal E_t[\ket\Phi\!\bra\Phi]$ is the Choi-Jamio{\l}kowski state
associated to the channel $\mathcal E_t$ and $\ket\Phi=\sum_{j=1}^N\ket{jj}/N$.
Indeed, using \eqref{helstrom} in \eqref{diamond} we may write for two generic channels,
$\mathcal E_0$ and $\mathcal E_1$, the following inequality
(also discussed in \cite{Watrous2018} Exercise 3.6)
\begin{align*}
	\|\mathcal E_0 -\mathcal E_1\|_\diamond  &=
	\max_\Psi \| (\openone\otimes\mathcal E_0-\openone\otimes\mathcal E_1)(\Psi)\|_1
	\cr &=
	\max_{E,\Psi} \sum_i |\Tr \left[E_i(\openone\otimes\mathcal E_0-\openone\otimes\mathcal E_1)(\Psi)\right]|
	\cr &=
	\max_{E,M} \sum_i |\Tr \left[E_i(M\otimes \openone)(\chi_{0}-\chi_{1})M^\dagger\otimes\openone\right]|
	\cr &=
	\max_{M} \|(M\otimes \openone)(\chi_{0}-\chi_{1})M^\dagger\otimes\openone\|_1
	\cr &\le \max_M \|M\|_\infty^2 \|\chi_{0}-\chi_{1}\|_1
	\cr &\le \max_M \|M\|_2^2 \|\chi_{0}-\chi_{1}\|_1
	\cr &= N \|\chi_{0}-\chi_{1}\|_1 ,
\end{align*}
where, without loss of generality, we have set $\ket{\Psi} ={M\otimes\openone}\ket\Phi$, so that
$M$ satisfies $\|M\|_2^2=\Tr M^\dagger M=N$ due to normalization. Indeed, from $\ket{\Psi} =
\sum_{ij} \Psi_{ij} \ket{ij}$ we may set $M_{ij}=\sqrt N \Psi_{ij}$. In the fifth line
we use the property $\|ABC\|_1 \le \|A\|_\infty \|B\|_1\|C\|_\infty$ (see \cite{Watrous2018}
section 2.3.1), and in sixth line $\|A\|_\infty \le \|A\|_2$.

The trace distance in Eq.~\eqref{diamondineq} could be measured experimentally as follows.
The state $\chi_{{\rm ensemble}}$ may be created by first generating path-entangled
photonic states, and then sending half of the entangled state inside the linear-optical network
that implements a quantum walk (or any other unitary). On the other hand, using simple
properties of Haar integrals we find $\mathcal \chi_{{\rm Haar}}=\openone/N^2$.
The resulting  experiment, although technically possible,
is made more challenging by the use of entangled photon sources and by the demand of complete
tomography, for reconstructing $\chi_{{\rm ensemble}}$.
In order to introduce a simpler quantity for experiments we expand $\chi_{\rm ensemble}$
into the basis $\ket i$, where $\ket i$ means that a photon is injected (or detected) in
the $i$th waveguide. Using the definition of the Choi-Jamio{\l}kowski state, we may write
\begin{align}
	\chi_{\rm ensemble} &= \sum_{ij,kl} \frac1N \bra k \mathcal E_{\rm ensemble}[\ket i\!\bra j]\ket l ~ \ket{ik}\!\bra{jl} =
	\nonumber
	\\ &= \sum_{ij,kl} \chi_{ik,jl} \ket{ik}\!\bra{jl}~,
\end{align}
where we have defined  $ \chi_{ik,jl} = \bra k
\mathcal E_{\rm ensemble}[\ket i\!\bra j]\ket l/N$. We may now split
\begin{equation}
	\chi_{\rm ensemble} = \chi_{\rm d} + \chi_{\rm od},
\end{equation}
where $\chi_{\rm d} = \sum_{ik} \chi_{ik,ik} \ket{ik}\!\bra{ik}$ is the diagonal part or $\chi_{\rm ensemble}$
and $\chi_{\rm od}$ its off-diagonal part.
Therefore, using the bound \eqref{diamondineq} and the triangle inequality we may write
\begin{equation}
	\|\mathcal E_{\rm ensemble} -\mathcal E_{\rm Haar}\|_\diamond
	\le N (D_{\rm d} + D_{\rm od})~,
\end{equation}
where $D_{\rm d}$ measures the distance of the diagonal elements from the uniform distribution,
\begin{align}
	D_{\rm d} &= \|\chi_d-\openone/N^2\|_1 = \sum_{ik} \left|\chi_{ik,ik}-\frac1{N^2}\right| =
					\\&= \sum_{ik} \frac1N \left|\sum_{l=1}^m \bra k U_l \ket i\!\bra i U_l^\dagger \ket k -\frac1N\right| =
					\\&= \frac1N \sum_{i,k=1}^N \left|\sum_{l=1}^m |\bra k U_l \ket i|^2  -\frac1N\right|,
\end{align}
and $D_{\rm od}=\|\chi_{\rm od}\|_1 \leq 2$ measures the strengths of the off-diagonal elements.

In order to define a simpler quantity for experimental tests, we consider focus on
\begin{equation}
	\mathcal N =
	\max_\Psi \| (\openone\otimes\mathcal E_{\rm ensemble}-\openone\otimes\mathcal E_{\rm Haar})(\Psi)\|_2,
	\label{diamond2}
\end{equation}
which is analogous to Eq.~\eqref{diamond}, but with the $\ell_2$-norm in place of the $\ell_1$-norm.
The maximization in the above norm is over all possible input states but, due to convexity of matrix norms,
we may restrict such maximization to pure states. Indeed, for any convex function $f(\rho)$ of a quantum state $\rho$,
we may write $f(\sum_i p_i \ket{i}\!\bra{i})\leq \sum_i p_i f(\ket{i}\! \bra{i})\leq  \max_i f(\ket{i}\!\bra{i})$,
so the maximum is always obtained with pure states. Without loss of generality, we can also choose
$\ket{\Psi} = \sqrt{\rho}\otimes\openone\ket{\Phi}$, namely as a purification of the state $\rho$.
Then it is simple to show that $\mathcal N^2 = \Tr[(\rho\otimes\openone)\chi(\rho\otimes\openone)\chi]$,
where $\chi=\chi_{\rm ensemble}-\chi_{\rm Haar}$, which
is convex over $\rho$. As such the maximum is obtained for pure $\rho$ or, in other terms, when
$\Psi$ in Eq.~\eqref{diamond2} is separable. Since $\mathcal E_{\rm Haar}(\Psi)=\openone/N$ we can write
\begin{equation}
	\mathcal N^2 = \max_{\ket{\psi}} \Tr[\mathcal E^2_{\rm ensemble}(\ket\psi\!\bra\psi)]-\frac1N
	= \max_{\ket\psi}\left[\mathcal N_{\rm d}^2 + \mathcal N^2_{\rm od}\right],
\end{equation}
where
\begin{align}
	\mathcal	N_{\rm d}^2 &= \sum_k \left[p_{\ket\psi}(k)-\frac1N\right]^2~,\\
	\mathcal	N_{\rm od}^2 &= \sum_{k,l:k\neq l}  |\bra k\mathcal E_{\rm ensemble}(\ket\psi\!\bra\psi)\ket l|^2,
\end{align}
are respectively the diagonal and off-diagonal contribution and
$p_{\ket\psi}(k) = |\bra k\mathcal E_{\rm ensemble}(\ket\psi\!\bra\psi)\ket k|$.

\subsection{Supplementary Note 2: The numerical analysis of the off-diagonal elements}

We consider the $\Delta \beta$ approach that produces the effective Hamiltonian $H_{\rm eff}$ given in Eq.(3) of the main text. The evolution of such stochastic quantum walks implement the operation $U\rho U^{\dagger}$. Experimentally, the diagaonal elements of $U\rho U^{\dagger}$ can be measured by the probability distribution, while the off-diagonal elements are difficult to measure, but still can be numerically calculated using the $\Delta \beta$ model.

We consider the parameter setting of $\beta_0$ and $C_{ij}$ to be consistent with those in the experiment. The square waveguide lattice has a nearest-neighbor coupling coefficient of 0.225 $\rm mm^{-1}$, and a next-nearest-neighbor coupling coefficient of around 0.09 $\rm mm^{-1}$, as has been characterized for the same waveguide spacing configuration in \cite{Tang2018}. For waveguides further apart, the coupling coefficients are small enough to be not considered. The absolute value of $\beta_0$ does not influence the numerical result. In this work, we use a $\beta_0$ value of 11660 $\rm mm^{-1}$.

Here we calculate the norm for the off-diagonal elements of $\mathbb{E}_i[U_{i}\rho U_{i}^{\dagger}]$ for different $\Delta \beta$ amplitudes, the scale of dephasing, and different numbers of samples. As shown in Fig. S1, increasing the number of samples and the $\Delta \beta$ amplitude could clearly improve the convergence to zero. This suggests that the off-diagonal part of the stochastic quantum walk in the photonic lattice model proposed in this work can reach the uniform distribution to meet the Haar randomness measure.

\subsection{Supplementary Note 3: The $\Delta \beta$ photonic approach}
Before introducing the method of $\Delta \beta$ generation, we would firstly explain how to generally measure the coupling coefficient $C$ using the directional coupler approach~\cite{Szameit2007}. In the directional coupler, light intensity in the two waveguides are denoted as $I(A)$ and $I(B)$, respectively. The intensity would vary sinusoidally with the evolution length $z$, that is, $I(A)=|{\rm cos}(Cz)|^2$, $I(B)=|i {\rm sin}(Cz)|^2$. Therefore,
\begin{equation}
	\frac{I(A)-I(B)}{I(A)+I(B)}=\frac{{\rm cos}^2(Cz)-{\rm sin}^2(Cz)}{{\rm cos}^2(Cz)+{\rm sin}^2(Cz)}={\rm cos}(2Cz).
\end{equation}
Experimentally, we prepare a group of directional couplers of different evolution lengths, and measure the light intensity at each of the two modes (see Fig.S2a). Note that we just use two straight waveguide to form the directional couplers and it does not have a curved part. This avoids unnecessary loss in the bending part and avoids the unwanted coupling in the bending region where two waveguides gradually go apart. We can fit the relationship between $\frac{I(A)-I(B)}{I(A)+I(B)}$ and $z$ with a suitable sinusoidal curve (see Fig.S2b), and the coupling coefficient $C$ can be obtained from the frequency of this sinusoidal curve.

Then $\Delta \beta$ can be measured using this directional coupler approach. One of the waveguide is written using a base speed $V_0=5~\rm mm/s$, and the other waveguide using a different detuned speed $V$ ($V-V_0=\Delta V$) that will lead to a detuned propagation constant $\Delta \beta$ on this waveguide. In the detuned directional coupler, the effective coupling coefficient $C_{\rm eff}$ can be obtained using the same coupling mode method as that for the normal directional coupler~\cite{Szameit2007}, but $C_{\rm eff}$ contains the detuning effect from $\Delta \beta$ through this equation~\cite{Lebugle2015}:
\begin{equation}
C_{\rm eff}=\sqrt{(\Delta \beta/2)^2+C_0^2},
\end{equation}
where $C_0$ is the coupling coefficient for a normal directional coupler. Therefore, $\Delta \beta$ can be calculated when $C_{\rm eff}$ and $C_0$ are both characterized.

We prepared 13 groups of samples with a base speed of $5~\rm mm/s$ and a detuned speed ranging from 5, 7.5, 10, 12.5, all the way to 35$~\rm mm/s$. This corresponds to a $\Delta V$ from 0 to $30~\rm mm/s$. Each group has 12 samples of the same speed setting but different evolution lengths such that we can fit the sinusoidal curve to get $C_{\rm eff}$ or $C_0$ (when $\Delta V=0$).
With these 156 samples, we then obtain $\Delta \beta$ (unit: $\rm mm^{-1}$) against $\Delta V$ (unit: $\rm mm/s$), as ploted in Fig.S2c. There can be a rough fit linearly: $\Delta \beta=0.02\times \Delta V$. The linear relationship is consistent with a previous report\cite{Caruso2016}. Knowing this, we can randomly generate $\Delta \beta$ of $0.01-0.4~\rm mm^{-1}$ by varying $\Delta V$ between $0.5-20~\rm mm/s$.

It is worth mentioning that in this experiment, we set a target $\Delta \beta$ value by strictly following the linear relationship with $\Delta V$. This might cause some systematic error as the real generated $\Delta \beta$ might go above  $0.02\times \Delta V$, e.g., a $\Delta V$ of $20~\rm mm/s$ might have generated a $\Delta \beta$ of $0.5~\rm mm^{-1}$ instead of $0.4~\rm mm^{-1}$. This may explain the slightly better convergence in experiment than in numerical analysis in Fig. 3b of the main text. Nonetheless, there is a good match of the convergence trend overall.

It is also worth noting that although we calibrate $\Delta \beta$, an absolute value of $\beta$ detuning, an increase in writing speed in fact corresponds to a decrease in $\beta$. In this work, plus or minus a $\beta$ detuning has almost no influence in the numerical analysis, but the sign may matter for many other work that require $\beta$ detunings such as the manipulation of site energies, so we believe it's necessary to clarify this for wide readers.

Additionally, our experimental limit for generating a $\Delta \beta$ can reach around $1.0~\rm mm^{-1}$ that requires a speed detuning of around $50~\rm mm/s$. This allows for rich manipulation of $\Delta \beta$ and can make a convergence at a feasible length on current photoic chips.

\subsection{Supplementary Note 4: Waveguide preparation}
For the 17 random settings, each has 8 different evolution lengths. As it is not convenient to measure the evolution patterns in the middle of the waveguide, we have to make 8 samples with a length of 1cm, 2cm, ..., 8cm, respectively. We have ensured the 8 samples follow the same random setting, for instance, $\Delta \beta$ profiles for the 4-cm-long sample are exactly the same with those in the first 4cm of the 8-cm-long sample.

Each sample, regardless of the difference in evolution length or random setting, has the same configuration, that is, a 5 $\times$5 square lattice with the nearest waveguide spacing of 15 $\rm \mu m$ in the vertical direction, and of 13.5 $\rm \mu m$ in the horizontal direction. This ensures the nearest-neighbor coupling coefficient to be consistent in different directions\cite{Tang2018}.

All the waveguides are fabricated using the femtosecond laser direct writing technique\cite{Crespi2013, Chaboyer2015, Tang2018}. We direct a 513-nm femtosecond laser (up converted from a pump laser of 10W, 1026nm, 290fs pulse duration, 1 MHz repetition rate) into a spatial light modulator (SLM) to shape the laser pulse in the temporal and spatial domain. The writing speed is precisely controlled for each segment to ensure the introduction of the target $\Delta \beta$ value. We then focus the pulse onto a pure borosilicate substrate with a 50X objective lens (numerical aperture:~0.55). Power and SLM compensation were processed to ensure the waveguide uniformity.

\subsection{Supplementary Note 5: Measure probability distribution from graphs.} All the experimentally measured evolution patterns with 17 different random settings are provided in Fig.S3-S10. When collecting the data from experiments, we obtained the corresponding ASCII file, which is essentially a matrix of pixels. We created a `mask' containing the pixel coordinate of the circle centre and the radius in pixels for each waveguide, and summed the light intensity for all the pixels within each circle using Matlab. The normalized proportion of light intensity for each circle represents the probability at the corresponding waveguide.

\subsection{Supplementary Note 6: The standard error of the norm}
According to the definition of the norm and our experimental setting, we know how we calculte our norm $N$:
\begin{equation}
	N =\sqrt{\sum_{j=1}^{25}(\frac{\sum_{i=1}^{17}p_{ij}}{17}-0.04)^2},
\end{equation}
where $p_{ij}$ is the probability of $j$th site in $i$th experiment. Replacing it with the mean probability $P_{j}=\frac{1}{17}\sum_{i=1}^{17}p_{ij}$, we have:
\begin{equation}
	N =\sqrt{\sum_{j=1}^{25}(P_{j}-0.04)^2}.
\end{equation}
Using the error propagation formula, we have($\Delta$ means standard deviation):
\begin{equation}
	\Delta N =\sum_{j=1}^{25}(\frac{\partial N}{\partial P_{j}}\Delta P_{j})^2,
\end{equation}
in which $\Delta P_{j}$ is easy to get. Meanwhile:
\begin{equation}
	\frac{\partial N}{\partial P_{j}}=\frac{P_{j}-0.04}{\sqrt{\sum_{j=1}^{25}(P_{j}-0.04)^2}}=\frac{P_{j}-0.04}{N},
\end{equation}
so totally we have:
\begin{equation}
	\Delta N =\sum_{j=1}^{25}(\frac{P_{j}-0.04}{N}\Delta P_{j})^2.
\end{equation}
Therefore, we can calculate $\Delta N$ knowing $N$, $P_{j}$ and $\Delta P_{j}$. Then the standard error of the mean is $\frac{\Delta N}{\sqrt{17}}$.

As for the theoretical result of 17 random samples, we can also similarly work out a standard error, but we consider there is a more influential factor, that is, we may not exactly fabricate a sample with a $\Delta \beta$ amplitude of 0.4${\rm mm}^{-1}$. Therefore, we set a possible range of theoretical results considering a 15$\%$ fluctuation of the $\Delta \beta$ amplitude. We show a shading area with its upper bound and lower bound corresponding to a $\Delta \beta$ amplitude of 0.34${\rm mm}^{-1}$ and 0.46${\rm mm}^{-1}$, respectively. When calculating the shading area, the random setting configuration is kept the same with that for $\Delta \beta$ amplitude of 0.4${\rm mm}^{-1}$, but just proportionally shrink it to 0.34${\rm mm}^{-1}$ or enlarge it to 0.46${\rm mm}^{-1}$, respectively.

\subsection{Supplementary Note 7: The norm of $M$ using $\Delta \beta$ photonic approach for a long evolution length}
As shown in Fig.S11, for a longer evolution length, with a $\Delta \beta$ amplitude of 0.4 $mm^{-1}$, the norm of $M$ for stochastic quantum walk becomes closer to zero. We further show in Fig.S12 that increasing the evolution length can effectively reduce the norm to an extremely small value and the norm would get stable after a certain evolution length. We see in Fig.S12 among lattices of different $\Delta \beta$ amplitudes, the one with a higher $\Delta \beta$ amplitude has faster convergence to the Haar measure, but given the same number of sample random settings, lattices show a consistent threshold of the norm regardless of the $\Delta \beta$ amplitude values. Comparing Fig.S4a and 4b that use 100 and 500 random settings, respectively, we see that increasing the number of samples can further lower down the norm value. This is consistent with what we get in Fig. S1 for off-diagonal elements.

There is not an explicit issue of complexity here in our hardware implementation. We simply build a long-enough chip with random phases and, according to our theory\cite{Banchi2017}, the evolution will converge to a Haar random evolution. Now, when is a chip "long-enough"? This is a hard question: in \cite{Banchi2017} we show that this depends on the spectral gap of a Liouvillian operator, and estimating gaps of many-body operators is hard, see\cite{Cubitt2015}. Nonetheless, there are some noise models whose Liouvillean is exactly solvable, like the one considered in \cite{Banchi2017}. The model considered in this experiment is a noisier version (noise on all waveguides vs noise on just a single waveguide), so the required length would be lower than our theory, which can be considered as an upper bound.

In our photonic setting, we see that generally the lattice with a large $\Delta \beta$ amplitude converges faster than those with a small $\Delta \beta$ amplitude, so we can take advantages of our capacity of large-scale photonic lattice and mature $\Delta \beta$ manipulation techniques to make an efficient implementation of Haar randomness. For instance, for a 5×5 lattice with a $\Delta \beta$ amplitude of 0.6$\rm mm^{-1}$, it can well converge to Haar measure at an evolution length within 10cm, which is experimentally practical.

We further show that the segment length also influences the required evolution length for convergence in Fig. S13. Segment length suggests how frequently we introduce the random $\Delta \beta$ setting, as the value of $\Delta \beta$ will be constant throughout each segment. For instance, 10cm can be 50 segments multiplies a segment length of 2mm, or 20 segments multiplies a segment length of 5mm. We regard the evolution comes to convergence when norm firstly drops below the convergence value, which is around 0.0003 in this case of a fixed number of samples at 100. In Fig. S13c, each dot of convergence length is obtained by averaging 10 groups of setting, each with a fixed number of samples at 100. The error bars in Fig. S13c are the standard deviation of the convergence length for that 10 groups.

In Fig. S13 we see that while varying the segment length would not alter the norm value at convergence, it does change the convergence length, which can reach a minimal value when choosing a proper segment length at around 5-10mm. In the main experiment we use a $\Delta \beta$ of 0.4$\rm mm^{-1}$ and a segment length $\Delta z$ of 2mm. If still using a $\Delta \beta$ of 0.4$\rm mm^{-1}$, but with a $\Delta z$ of 5mm instead 2mm, it can converge at a faster speed, within a feasible chip length of 9.5cm. Changing a segment length is not of any experimental challenge, especially when enlarging it, so even a $\Delta \beta$ of 0.4$\rm mm^{-1}$ can converge at a reasonably short chip length with properly selected segment length.

\subsection{Supplementary Note 8: Convergence for two or multiple photons}
For two travelling particles injected in modes $i$ and $i'\neq i$ we get at the end of the chip the following state
\begin{equation}
	\ket{\psi(U)} = \sum_{kk'} U_{ki} U_{k'i'} a_{k}^\dagger a_{k'}^\dagger\ket 0.
\end{equation}
Setting $\rho(U)=\ket{\psi(U)}\!\bra{\psi(U)}$ and $\rho=\rho(I)$, where $I$ is the identity matrix,
the $\ell_2$ norm becomes
\begin{align}
	\mathcal N_2 =& \left\|\mathbb E_{s}[\rho(U_s)] - \mathbb E_{U\sim\rm Haar}[\rho(U)]\right\|_2^2  \nonumber
							\\=& \mathbb E_{ss'} \Tr[\rho(U_s)\rho(U_{s'})] + \mathbb E_{U,V\sim\rm Haar} \Tr[\rho(U)\rho(V)] -\\&-
	2\mathbb E_{s}\mathbb E_{U\sim\rm Haar}\Tr[\rho(U_s)\rho(U)]\nonumber
							\\=& \mathbb E_{ss'} \Tr[\rho(U_{s'}^\dagger U_s)\rho)] + \mathbb E_{U,V\sim\rm Haar} \Tr[\rho(V^\dagger U)\rho] -\\&-
	2\mathbb E_{s}\mathbb E_{U\sim\rm Haar}\Tr[\rho(U^\dagger U_s)\rho]  \nonumber
							\\=& \mathbb E_{ss'} \Tr[\rho(U_{s'}^\dagger U_s)\rho)] - \mathbb E_{U\sim\rm Haar} \Tr[\rho(U)\rho]
							\label{N2},
\end{align}
where in the last line we used the property that $UV$ and $VU$ are Haar distributed if
$U$ is Haar distributed and $V$ is any unitary matrix.
To evaluate the norm we note that
\begin{align}
	\label{asd1}
	\Tr[\rho(U)\rho] &= |\bra{\psi(I)}\psi(U)\rangle|^2  \\&= \left|\sum_{kk} U_{ki}U_{k'i'} (\delta_{ki}\delta_{k'i'}+\delta_{ki'}\delta_{k'i})\right|^2
	\nonumber
\\
									 &=	 \left| U_{ii}U_{i'i'}+ U_{ii'}U_{i'i} \right|^2.
	\nonumber
\end{align}
By explicit computation using results from \cite{Richard2019}
we find $\mathbb E_{\rm Haar}[\rho(U)\rho] =\frac{2}{N(N+1)}$. Inserting the above results in Eq.~\eqref{N2} we get
\begin{align}
	\mathcal N_2 =& \frac1{m^2}\sum_{s,s'}
	\left| (U_{s'}^\dagger U_s)_{ii}(U_{s'}^\dagger U_s)_{i'i'}+ (U_{s'}^\dagger
	U_s)_{ii'}(U_{s'}^\dagger U_s)_{i'i} \right|^2 - \frac2{N(N+1)}.
\end{align}
As with the single photon experiment, we now focus on measurable quantities, with a similar meaning of the above
formal norm.
By a similar calculation to \eqref{asd1} we find that the probability that the two photons
come out from modes $j$ and $j'$, when they were injected in $i$ and $i'\neq i$ is
\begin{equation}
	p_{j,j'}(U,i,i') = \frac{\left| U_{ji}U_{j'i'}+ U_{ji'}U_{j'i} \right|^2}{1+\delta_{j,j'}},
\end{equation}
and integrating $U$
\begin{align}
	p_{j,j'}^{\rm Haar} = \mathbb U_{\rm Haar} [p_{j,j'}(U,i,i')] = \frac2{N(N+1)}.
\label{p}
\end{align}
The figure of merit for experiments is therefore
\begin{equation}
	\| \mathbb E_{s} [p(U_s,i,i')] - p^{\rm Haar}\|_2^2 =
	\sum_{j,j'=1}^N \left(\frac 1m \sum_{s=1}^m p_{j,j'}(U_s,i,i') - p_{j,j'}^{\rm Haar} \right)^2,
\end{equation}
where $p_{j,j'}(U_s,i,i') $ can be measured in experiments with sample $s$ when two photons are injected
in modes $i$ and $i'$ and detected in modes $j$ and $j'$.
Another measurable quantity is the intensity on mode $j$
\begin{equation}
	I_{j} (U,i,i') = 2p_{j,j}(U,i,i')+\sum_{j' \neq j} p_{j,j'}(U,i,i').
\end{equation}
Using Eq.~\eqref{p} we find a uniform average intensity on an modes
\begin{equation}
	I_{j}^{Haar}=\frac{2}{N}.
\end{equation}

We now numerically investigate the convergence when injecting two indistinguishable photons. We denote $M_2=\mathbb E_{s}[I_{j}]-\frac{2}{N}$, which is consistent with $\mathbb E_{s}[p(U_s,i,i')]-\frac2{N(N+1)}$. According to the Haar measure, the norm of $M_2$ would converge to zero after a certain evolution length. As shown in Fig.S16, there is indeed a convergence to an ultra small value, which is very similar to the single photon scenario presented in Fig.S12.

With a closer comparison with Fig.S12, we can clearly see that the convergence follows a same dependence on $\Delta \beta$ amplitude and number of groups as with the scenario in $q=1$. That is, increasing the $\Delta \beta$ amplitude can speed up the convergence without altering the norm lower bound, while increasing the number of samples can further lower down the norm. Another similarity is that, both $q=2$ and $q=1$ converge at the similar evolution length, e.g. around 30 cm using a $\Delta \beta$ amplitude of 0.3${\rm mm}^{-1}$. A small difference lies in that the lower bound for $q=2$ is slightly higher, e.g., 0.0008 for $q=1$ and 0.001 for $q=2$ with 500 samples. In fact, the lower bound is a less relevant quantity, as it can slightly compensated by increasing the number of samples, and meanwhile, if we want to run a boson sampling experiment we will always use a single sample.

We suggest, from an analytical point of view, for the spectral norm we predict no dependence on $q$, because the spectral norm depends only on the largest singular value, not on the dimension of the operators. The numerical analysis of the experimental data uses on the other hand the Frobenius norm. The latter depends on all singular values and hence may scale with the dimension of the Hilbert space, which is larger for more photons. We suggest this scaling factor should be linear in the dimension, so not too big for the current experiment. Just as what we have numerically shown, with 500 samples, $q=1$ reaches 0.0008 and $q=2$ reaches 0.001. In order to make $q=2$ reach 0.0008, we may use a few more samples. As can be inferred by comparing Fig.16a and b, adding the sample number would slightly enhance the required evolution length, only by a very marginal value. Therefore, we suggest this experimental approach for the convergence to Haar measure is scalable for those beyond single photon scenarios.

More generally we may write the intensity as
\begin{equation}
	I_{j}(U) = \langle \psi (U)|a_{j}^{\dagger} a_{j}|\psi (U)\rangle = \sum_{kk'} U_{k,j}^{*}U_{k',i}\langle \psi |a_{k}^{\dagger} a_{k'}|\psi \rangle,
\end{equation}
so there is:
\begin{equation}
	I_{j}^{Haar} (U) =\sum_{k} \frac{1}{N}\langle \psi (U) |a_{k}^{\dagger} a_{k}|\psi (U) \rangle = \frac{q}{N},
\end{equation}
where $q$ is the number of photons.
As for the uniform probabilities we may use Theorem 4.3 in \cite{Richard2019} to write
\begin{equation}
\begin{split}
	p_{j_{1},j_{2},…,j_{q}}^{Haar} =\frac{1}{C_{N+q-1}^{q}} \langle j_{1},…,j_{q}|P_{s}|j_{1},…,j_{q}\rangle \\ = \frac{1}{C_{N+q-1}^{q}} \equiv \prod_{k=1}^{q} \frac{k}{N+q-k},
\end{split}
\end{equation}
where $P_{s}$ is the projection onto the symmetric subspace, and in the second line we used the fact that multiphoton Fock state are normalized and symmetric. For $q=1$ we get the usual $1/N$, while for $q=2$ we recover Eq.~\eqref{p}.


\section{Pulsars in globular clusters}\label{sec:pulsars}

Globular clusters contain a very rich population of pulsars with 232 known in 36 cluster\footnote{An up-to-date count can be found in the webpage \url{https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html}} and this population  is quite peculiar when compared with the one in the Galactic disk. The abundance of pulsars per unit of stellar mass in globular clusters is $\sim 100$ times higher than that of the disk and the vast majority of the pulsars are  millisecond pulsars (MSP) with rotation periods shorter than a few tens of ms.

The difference between the population of pulsars in globular clusters and in the Galactic disk is caused by the starkly different environments. Globular clusters have no sign of recent star formation and have very high stellar densities in their cores that can reach up to $\sim 10^{4-5}$ $\rm M_{\odot} \, pc^{-3}$ compared to $\sim 1-10$ $\rm M_{\odot} \, pc^{-3}$ in the Solar neighbourhood. This high stellar density favours dynamical encounters between stars that can lead to exchange encounters where neutron stars become part of new binary systems. The exchanges of companions increase the chance of the neutron star going through an accretion phase boosting the number of low mass X-ray binaries (LMXB) that evole into MSPs. Once a neutron star has been recycled into an MSP, it is likely to remain active for a time even longer than the Hubble time due to the low magnetic field.

Due to the peculiar environment within globular clusters, the possibility of going through multiple stages of recycling and the lack of recent star formation, the majority of pulsars have very fast rotational periods usually smaller than 20 ms. The fastest pulsar known, PSR J1748$-$2446ad\cite{Hessels2006}, with a rotational period of 1.395 ms resides in the globular cluster Terzan 5. Not all of the pulsars in globular clusters are, however, very fast. There is a small population of pulsars where the accretion process was interrupted before completion \cite{Verbunt2014}, possibly by external encounters, with period between 20-100 ms called ‘mildly recycled'.
Furthermore, there is, in some specific clusters, a small number of slow pulsars with periods over 100 ms with the record of 1004.04 ms being held by B1718$-$19A in NGC 6342 \cite{Lyne1993}. The presence of these apparently ‘young' pulsars is in contradiction with the known formation scenarios and the lack of new star formation. Apart from dynamic disruption of LMXB \cite{Verbunt2014} other possible explanations have been put forward like accretion-induced collapse, direct collisions with a main sequence star \cite{Lyne1996}, or electron capture supernova of an OMgNe white dwarf \cite{Boyles2011}.

The MSPs populations of globular clusters is also quite different from the one in the Galactic disk.
The percentage of isolated to binary MSPs can vary a lot between clusters with some like NGC 6624 and NGC 6517 having mostly isolated pulsars and others like M62 and M5 having mostly binary pulsars. This is linked to the dynamical parameters of the clusters\cite{Verbunt2014}.

The multiple dynamical interactions create a population of binary MSPs that are in some cases very different from the Galactic disk population. Eclipsing systems with non-degenerate companions like ‘black widows' and ‘redbacks' are common \cite{Lyne1993,Camilo2000,Ransom2005}. Binary systems can acquire high eccentricities and massive companions through successive interactions or through exchange encounters like M15C\cite{Anderson1990}, NGC 1851A \cite{Freire2004}, NGC 6544B \cite{Lynch2012}, and NGC6652A\cite{Decesar2015} . A very peculiar triple system with a white dwarf and a Jupiter-mass companion has also been observed in a globular cluster \cite{Thorsett1999,Sigurdsson2003}.

\subsection{Science with pulsars in globular clusters}

The scientific achievements made possible by this unique population of pulsars together with the possibility of observing a large number of pulsars in the same telescope beam have made globular clusters a high priority target in the past decades.

The internal composition of neutron stars is determined by the equation-of-state that can be probed by measuring their masses. This can be achieved through pulsar timing by measuring the so-called post-Keplerian parameters, which quantify relativistic orbital effects \cite{Antoniadis2013}. One of these parameters is the relativistic precession of the periastron that is relatively easily accessible in highly eccentric systems, and can be used to place constraints on the total mass of the systeem assuming general relativity. Thanks to this method, some neutron stars with the potential of breaking the current mass limit have been identified\cite{Freire2008,Freire2008b}. The importance of determining neutron star masses is not limited to the value of the maximum mass, also the minimum mass and the distribution of masses can be used to test various formation and recycling scenarios\cite{Ridolfi2019}.

Another avenue of research that can be pursued in globular clusters regards the maximum rotation speed. Currently the record of 1.395 ms is held by J1748-2446ad in Terzan 5\cite{Hessels2006} but even faster pulsars might be present in clusters and set more constraints on the equation of state.

Studies of eclipsing pulsars like ‘black widows' and ‘redbacks' are essential to study the possible evolution scenarios of MSPs and to constrain the properties of the companions\cite{vanKerkwijk2000,Bogdanov2005}.

If we change focus from specific pulsars to the entire population of pulsars in a single cluster, the science perspective broadens to include the properties of the globular clusters themselves. The gravitational effects of the globular clusters on the pulsars are accessible through the derivatives of the rotational period\cite{Phinney1993}.
With a significant number of pulsars, this can give us access to the structural parameters of the clusters and explore the possibility of the presence of an intermediate mass black hole in the center\cite{Prager2017,Perera2017,Abbate2018}.

Using the dispersion measures (DM) and estimating the 3-d position of the pulsars inside the cluster, it becomes possible to look for diffuse gas. This method has led to the first detection of ionized gas in a globular cluster\cite{Freire2001,Abbate2018}. If we add the information coming from the rotation measure, it becomes possible to probe also the magnetic fields. This led to the possible discovery of a magnetized outflow from the Galaxy in the direction of the globular cluster 47 Tucanae\cite{Abbate2020}.

The numerous scientific possibilities suggest that globular clusters pulsars should still be considered extremely profitable targets for fututre research. Simulations show that a significant population of potentially interesting pulsars is still undiscovered in the clusters\cite{Bagchi2011,Turk2013}. These could include MSP-black hole systems, which would be invaluable for tests of gravitational theories\cite{Liu2014}. However, despite the efforts, the number of new discoveries in the last ten years has been very low. The ones that were discovered were found in archival data using new search techniques \cite{Pan2016,Cadelano2018,Andersen2018}. This suggests that we had reached the sensitivity limit of the current observing facilities. To uncover the numerous and potentially interesting pulsars new and more sensitive facilities are needed. One such facility is FAST\cite{Nan2011} in China that became operative in 2016 and has already discovered 32 pulsars in globular clusters\cite{Pan2020,Wang2020,Pan2021}. Another new facility that recently came online and has the potential to revolutionize the field is MeerKAT\cite{Booth2012} in South Africa.

\section{Advantages of the MeerKAT Radio Telescope}\label{sec:meerkat}

The MeerKAT radio telescope is an interferometer composed of 64 antennas each with a diameter of 13.5 m. The maximum baseline within antennas is 8 km. Thanks to its location in South Africa at a latitude of $-30^{\circ}$, it can observe the entire Southern sky. When observing with all antennas it can reach telescope gains four times higher than Parkes, the only other radio telescope capable of observing pulsars in cm-wavelengths below a declination of $-45^{\circ}$. The gain is also 1.4 times higher than the Green Bank Telescope giving MeerKAT a significant advantage when looking at the central regions of the Galaxy that are observable at high elevations for long hours.

MeerKAT is currently able to observe in  UHF-band (544-1088 MHz), L-band (856-1712 MHz) showing remarkably RFI-clear bands. In near future also the S-band (1750-3500 MHz) receivers will be available for observations. A description of the capabilities and potentialities of MeerKAT as a pulsar observatory can be found in Bailes et al. 2020\cite{Bailes2020}.

Of the MeerKAT Large Science Projects (LSPs) that are currently using MeerKAT for observations, two include pulsars as the main target: MeerTIME\cite{Bailes2016} and TRAPUM\cite{Stappers2016}.

\subsection{MeerTIME}

MeerTIME\footnote{http://www.meertime.org} is an LSP based on timing already known pulsars with the high sensitivity of MeerKAT that started collecting data in early 2019. It has four main scientific goals: timing relativistic binary pulsars in order to measure the masses of the systems and test theories of gravity\cite{Kramer2021}; the Thousand Pulsars Array that monitors a wide variety of pulsars to study their geometries, magnetospheres and interstellar effects \cite{Johnston2020}; a pulsar timing array that monitors high precision pulsars to contribute to the search for nanohertz gravitational waves [Spiewak et al. submitted]; time pulsars in globular clusters\cite{Ridolfi2021}.

The project plans to observe globular cluster pulsars in order to study the eclipses, measure the masses of the pulsars, and study the properties of the globular clusters themselves.

MeerTIME can observe up to four tied array beams using the Pulsar Timing User Supplied Equipment (PTUSE) machines as the main data acquisition system. These beams can either be in “timing” mode where a pulsar is folded in real time using an ephemeris or in “search” mode where the data is saved in a “filterbank” style file and can be folded later or searched for new pulsars or single pulses from previously known ones. For both modes the observations are done with very high time resolution (9 $\mu$s), full polarimetric information, and real-time coherent de-dispersion.

Observing globular clusters has the advantage that a large number of pulsars fall within the same beam of the telescope and we are therefore capable of observing them simultaneously if the data is recorded in “search” mode. However, the beam of MeerKAT using the full array with a baseline of 8 km would be too small to observe all of the pulsars in some globular clusters so in most cases only the central 44 antennas with a maximum baseline of 1 km are used.

\subsection{TRAPUM}

The TRAnsients and PUlsars with MeerKAT (TRAPUM\footnote{http://www.trapum.org}) LSP focuses on searching for new pulsars. It searches for pulsars in unidentified Fermi sources, supernova remnants, pulsar wind nebulae, in the Galactic plane, in nearby galaxies and in globular clusters. The globular cluster working group is focused on searching both in clusters with known pulsars (which allows de-dispersion at the known DM of the cluster) and in clusters with no known pulsars. The project started taking data in April 2020.

TRAPUM uses Filterbanking Beam-former User Supplied Equipment (FBFUSE) and the Accelerated Pulsar Search User Supplied Equipment (APSUSE), a 60 node computing cluster, to synthesize up to $\sim 400$ beams and record incoherently de-dispersed search-mode data for each one. This allows us to cover a large portion of the sky (up to $\sim 1 - 4$ arcmin in radius depending on the configuration) while keeping the full sensitivity of MeerKAT. In this way it is possible to cover the entire cluster with one pointing and search for pulsars in the center and in the outskirts of the cluster simultaneously. However, due to the very high data rate, it is not possible to retain full temporal and spectral resolution nor full polarimetric information.

Because of the complementarity of science goals, the globular cluster working groups of MeerTIME and TRAPUM have decided to collaborate sharing observing time and resources. This collaboration has already resulted in two publications \cite{Abbate2020b,Ridolfi2021}.

\section{New discoveries}\label{sec:discoveries}

The joint collaboration has discovered 36 new pulsars in globular clusters\footnote{ the full up-to-date list of discoveries can be found at http://www.trapum.org/discoveries.html}. The discoveries are summarized in Table \ref{table1} and \ref{table2} and can be split chronologically into two groups: the ones made with the PTUSE machine before TRAPUM became operative and the ones made possible by the APSUSE machine of TRAPUM.

\begin{table}
\tbl{Properties of the pulsars discovered by MeerTIME/TRAPUM using the PTUSE machine}
{\begin{tabular}{@{}cccc@{}}
\toprule
Pulsar name & Period & DM & notes \\
&  (ms) & (pc cm$^{-3}$) \\
\colrule
Ter 5 an  & 4.8023 & 237.7 &  	Binary \\
M62G	  & 4.6081 & 113.6 &  	Binary \\
NGC 6522D & 5.5369 & 192.7 &  	Isolated \\
NGC 6624G & 6.0912 & 86.20 &  Binary \\
47 Tuc ac & 2.7456 & 24.46 &  Ecl. Binary \\
47 Tuc ad & 3.7460 & 24.41 &  Ecl. Binary \\
NGC 6752F & 8.4854 & 33.20 &  Isolated \\
NGC 6624H & 5.1301 & 86.85 &  Isolated \\
NGC 6440G & 5.2157 & 219.7 &  	Isolated \\
NGC 6440H & 2.8486 & 222.6 &  	Binary \\
M28M	  & 9.5689 & 119.23 &	Binary \\
\botrule
\end{tabular}
}
\label{table1}
\end{table}

\begin{table}
\tbl{Properties of the pulsars discovered by MeerTIME/TRAPUM using the APSUSE machine}
{\begin{tabular}{@{}cccc@{}}
\toprule
Pulsar name & Period & DM  & notes \\
&  (ms) & (pc cm$^{-3}$) \\
\colrule
NGC 6752G &	4.7902	 & 33.27	 & Isolated \\
M28N	  & 3.3429	 & 119.33 & 	Binary \\
NGC 6342B &	2.5683	 & 71.44	 & Isolated \\
NGC 1851B &	2.8162	 & 52.07	 & Isolated \\
NGC 1851C &	5.5648	 & 52.05	 & Isolated \\
NGC 1851D &	4.5543	 & 52.17	 & Binary\\
NGC 1851E &	5.5952	 & 51.95	 & Binary \\
NGC 1851F &	4.3294	 & 51.63	 & Binary \\
NGC 1851G &	3.8028	 & 51.01	 & Binary \\
NGC 1851H &	5.5061	 & 52.26	 & Binary \\
NGC 6752H &	2.0155	 & 33.25	 & Isolated\\
NGC 6752I &	2.6582	 & 33.34	 & Isolated\\
NGC 6624I &	4.3195	 & 87.35	 & Isolated\\
NGC 6624J &	20.8995	 & 86.85	 & Isolated\\
NGC 6441E &	251.1587 & 221.16 & 	Isolated\\
NGC 6441F &	6.0006	 & 228.22 & 	Binary\\
NGC 6441G &	5.3522	 & 229.20 & 	Isolated\\
NGC 1851I &	32.6538	 & 52.42	 & Binary\\
NGC 1851J &	6.6329	 & 52.06	 & Isolated\\
NGC 1851K &	4.6920	 & 51.93	 & Isolated\\
NGC 1851L &	2.9586	 & 51.23	 & Binary\\
NGC 1851M &	4.7977	 & 51.66	 & Isolated\\
NGC 1851N &	5.5679	 & 51.11	 & Isolated\\
J1823-3022 & 2497.72  & 96.9      & Isolated\\
NGC 6624K & 2.7686  & 87.46      & Isolated\\
\botrule
\end{tabular}
}
\label{table2}
\end{table}

\subsection{MeerTIME discoveries}
The first batch of discoveries made with observations taken by MeerTIME is already published\cite{Ridolfi2021}. This publication contains the description and the properties of 8 new millisecond pulsars found in six different globular clusters. Three of them are isolated pulsars while the rest are found in binary systems. The profiles of these pulsars are shown in Figure \ref{figure1}.

\begin{figure}[h]
\begin{center}
\includegraphics[width=5in]{Figures/Discoveries_integrated_profiles_fullres_n128.pdf}
\end{center}
\caption{Profiles of the eight new pulsar discoveries published in Ridolfi et al. 2021\cite{Ridolfi2021}.}
\label{figure1}
\end{figure}

These discoveries include two new eclipsing millisecond pulsars 47 Tuc ac and 47 Tuc ad. Fitting for the orbits of these pulsars suggest that 47 Tuc ac is a ‘black widow' pulsar while 47 Tuc ad is a ‘red-back'.  Both of these discoveries were then re detected in previous Parkes observations but due to scintillation effects and eclipses the detections were too few to allow for deriving phase-connected timing solutions.

Possibly the most interesting pulsar in this set of discoveries is NGC 6624G. This pulsar is found in a binary with orbital period of 1.54 days and eccentricity of 0.38. Because of the very high eccentricity it was possible to use only 11 months of data to measure the rate of advance of periastron $\dot \omega =0.217 \pm 0.004$ deg yr$^{-1}$. Assuming that this advance is caused entirely by relativistic effects and using formulas derived in General Relativity (GR) it is possible to derive the total mass of the system, $M_{tot}= 2.65 \pm 0.07 \rm M_{\odot}$. This parameter alone is not enough to measure precisely the mass of the pulsar itself, see Figure \ref{figure2}, but suggests that the neutron star is particularly heavy. Measurements of further post-keplerian parameters will be needed to accurately determine the mass. The companion mass is also particularly heavy suggesting that it was not the original companion that recycled the pulsar but a new companion gained in an exchange encounter. This is the type of secondary exchange products that we have been looking for.

\begin{figure}[h]
\begin{center}
\includegraphics[width=5in]{Figures/mass_mass_NGC6624G.pdf}
\end{center}
\caption{Mass-inclination and mass-mass diagrams for NGC 6624G. The main square panels depict the $\cos i$-$M_{\rm c}$ and $M_{\rm c}$-$M_{\rm p}$ planes. The grey areas in the former are excluded by the requirement that the mass of the pulsar is positive, the grey areas in the latter are excluded by the mass function and the requirement that $\sin i \leq  1$. The red lines depict the masses consistent with the measurement of $\dot{\omega}$ and its $\pm  1$-$\sigma$  uncertainties, under the assumption that this effect is dominated by the GR contribution and that GR is the correct theory of gravity.  The contours include 68.3 and 95.4\% of a 2-D probability distribution function (pdf) derived from the $\chi^2$ of \texttt{TEMPO}\ fits that assumed all GR effects to be according to the masses and orbital inclination at each point. The side panels show the probability density functions for the $\cos i$ (top left), $M_{\rm p}$ (top right) and $M_{\rm c}$ (right) derived by marginalizing the aforementioned 2-D pdfs. we obtain a median for the pulsar mass of $2.1$~$\rm M_{\odot}$, but with a tail of probability  that extends to lower masses: there is a 31\% probability of $M_{\rm p} < 2$~$\rm M_{\odot}$\ and a 3.8\% probability of $M_{\rm p} < 1.4$~$\rm M_{\odot}$. Thus, either the pulsar is very massive, or it has a massive companion; the system was likely formed in an exchange encounter. Credits: Ridolfi et al. 2021\cite{Ridolfi2021}.}
\label{figure2}
\end{figure}

Another newly discovered pulsar that has the potential of being particularly massive is Ter 5 an. This binary has an eccentricity of 0.0066, two orders of magnitude higher than the average millisecond pulsar binaries found in the Galactic field. It is likely caused by distant stellar encounters in the dense cluster environment. Also here, owing to the recovery of the system in data from the Green Banks Telescope and the resulting long baseline, it was possible to measure an advance of the periastron with a rate of $\dot \omega =0.009 \pm 0.001 $ deg yr$^{-1}$. With the assumption that it fully caused by GR effects, the total mass becomes $M_{tot}= 2.95 \pm 0.52 \rm M_{\odot}$. The mass-inclination and mass-mass diagrams are shown in Figure \ref{figure3}. The uncertainties in this case are higher allowing for a greater range in pulsar mass.

\begin{figure}[h]
\begin{center}
\includegraphics[width=5in]{Figures/mass_mass_Ter5an.pdf}
\end{center}
\caption{Same mass-inclination and mass-mass diagrams as in Figure \ref{figure2} but for Ter 5 an. The estimated median pulsar mass is $2.13$~$\rm M_{\odot}$\ and median companion mass is $0.75$~$\rm M_{\odot}$. Credits: Ridolfi et al. 2021\cite{Ridolfi2021}.}
\label{figure3}
\end{figure}

The other pulsars discovered by MeerTIME are in the process of being published. Of these, one is particularly interesting: NGC 6440H. It is a binary pulsar with a very light companion. Assuming a pulsar mass of 1.4 $\rm M_{\odot}$, the companion has a minimum mass of 0.0062 $\rm M_{\odot}$. This very low mass of the companion suggests that it could be brown dwarf or a planet. If the planetary nature of the companion is confirmed, this would be only the second pulsar with such a companion \cite{Thorsett1999,Sigurdsson2003}.

\subsection{TRAPUM discoveries}

Only a small fraction of the observations of globular clusters with the APSUSE machine with TRAPUM have currently been processed. The large volumes of data produced in these observations require a large processing power to be properly analyzed. Despite this, the new pulsar discoveries are already more than double than those made using the PTUSE machines and only the central 44 antennas. The list is shown in Table \ref{table2}.

More than half of these discoveries come from just one globular cluster NGC 1851. This is a core-collapsed cluster that was known to contain only one pulsar NGC 1851A. TRAPUM observations in just the central beam have confirmed the presence of 13 new pulsars. This large population makes NGC 1851 the third globular cluster for number of known pulsars together with M28. Once phase connected ephemerides are determined for all these pulsars it will be possible to probe the potential well and test the presence of gas and magnetic fields within the cluster.

Two of these new discoveries, NGC 1851D and NGC 1851E, are binary systems with very high eccentricities. NGC 1851E has an orbital period of around 7 days and an eccentricity of $e=0.71$. This high value means that in the near future it will be possible to measure the rate of precession of the periastron to high accuracy and therefore the total mass of the system. Early fits to the orbit suggest that, assuming a pulsar mass of 1.4 $\rm M_{\odot}$, the companion has a minimum mass of 1.5 $\rm M_{\odot}$. This suggests that there is a high probability of this being another double neutron star system.

The vast majority of the discoveries are fast spinning MSPs but there are some exceptions: NGC 6624J and NGC 1851I have periods respectively of 20.89 and 32.65 ms, that suggest that they are mildly recycled pulsars. Another exceptions is NGC 6441E, a slow pulsar with a period of 251.15 ms. Particularly interesting is also J1823-3022 with a period of 2497.72 ms. This potentially record breaking pulsar has been found in the outskirts of the globular cluster NGC 6624 and has a DM that is higher than the rest of the pulsars in the cluster. Because of these reasons it is not clear if it is really a member of the cluster or if it is a Galactic disk pulsar positioned along the same line of sight. All of these slow pulsars are in core-collapsed clusters validating the correlation between spin periods and dynamical state of the cluster suggested previously \cite{Verbunt2014}.

These early results from a fraction of the observations suggest that many new discoveries are possible by analyzing the data that has already been captured.

\section{Future plans}\label{sec:future}

Early observations have showed that the future of the two projects MeerTIME and TRAPUM can prove revolutionary for the science of globular cluster pulsars.

MeerTIME will focus on timing the already discovered pulsars in order to measure the masses of the most interesting binaries and characterise the eclipses of the ecplising ‘black-widows' and ‘redbacks'. Furthermore, by improving the timing models of all pulsars it will be possible the study the properties of the clusters like the gas content, magnetic field and the presence of intermediate mass black holes.

TRAPUM, on the other hand, will continue to search for new pulsars broadening the scope in order to include clusters with no previously known pulsars. The searches will involve all of the available bands, going from the UHF band, where the pulsars are brighter but the interstellar medium is more noisy to the L band and the S band where the effects of scattering and the interstellar medium are minimized.

In closing, the last few years have shown an almost unprecedented growth in the number of pulsars in globular clusters going from 150 in 2018 to over 230 in 2021. The largest contributor is the MeerKAT telescope with 36 but a very close second is FAST with 32 new discoveries. These numbers are likely to increase in the following years and with them also the possibility of finding even more exotic systems like a double millisecond pulsar binary or a pulsar-black hole binary.

\bibliographystyle{ws-procs961x669}