\section*{Acknowledgments}

\section*{Supplementary materials}

\subsection{Optical setup}
\label{MM:optical_setup}

The experimental setup is displayed in Fig.~\ref{fig:setup}. The beam of a continuous wave diode laser (Eagleyard EYP-DFB-0785), operating at 785 nm, is sent through a single-mode fibre (Thorlabs P3-780PM-FC-1) in order to obtain Gaussian spatial profile. A half-wave plate (HWP) and a polarizing beam splitter (PBS) further split it into two paths, which will correspond to the local oscillator (LO) and signal of the heterodyne detector (HD).

In the LO path, the laser beam is magnified by a telescope ($f_1$ = 50 mm, $f_2$ = 200 mm) in order to fully and (almost) uniformly illuminate the screen of a reflective phase-only liquid-crystal-on-silicon SLM. The incident beam hits the SLM screen almost perpendicularly, so that the angle between the incoming and reflected waves is smaller than 5 degrees, and polarized parallel to the SLM's LC director axis, so that all the incident light is modulated. The SLM holograms output the desired phase and amplitude profiles in the first order of diffraction, which is selected by a telescope ($f_1$ = 250 mm, $f_2$ = 100 mm) and an iris at its focal plane.

In the signal path, the beam is sent through an acousto-optics modulator (AOM, Isomet 1205C-1), driven at 92.05 MHz. The produced first diffraction order mode is incident upon a DMD (DLP LightCrafter 6500), which modulates its amplitude with the binary bitmap of the object to be imaged. The imaging system aperture is set by an iris whose diameter is set to $3.5\pm0.03$ mm and placed at a distance of $245.5\pm1$ cm from the DMD, corresponding to an optical system of numerical aperture (NA) of $(7.1 \pm 0.07) \times$ 10$^{-4}$. The iris diameter was measured by placing a metal ruler next to it and imaging both these objects with a high-resolution camera.

After the iris, a set of three telescopes magnify and collimate the signal beam in order to match the waist and wave-vector of the LO beam in the 0th order HG mode. The first telescope ($f_1$ = 100 mm, $f_2$ = 50 mm) collects the light transmitted through the iris. In the second ($f_1$ = 50 mm, $f_2$ = 100 mm) and third telescopes ($f_1$ = 75 mm, $f_2$ = 75 mm), the ``eyepiece" lenses are mounted on translation stages to independently control the signal beam’s diameter and divergence, hence enabling mode matching to the LO. The signal and LO paths are reunited at a PBS, whose output beams feed the photodetectors of a homemade balanced detector \cite{kumar2012homodyne,masalov2017}.

\subsection{Spatial light modulator}
\label{MM:SLM}

The SLM (Hamamatsu X13138-02, 1272$\times$1024 pixel resolution and 12.5 $\mu \textup{m}$ pixel pitch) modulates the phase of an incoming optical wave by controlling the effective refractive index of the liquid crystal layer in each of its pixels. In order to generate  Hermite-Gaussian modes, we display a phase grating (``hologram") on the SLM screen. Setting the grating’s depth and offset at each point of the SLM surface allows us to generate any desired complex spatial profile in the first order of diffraction in the reflected wave \cite{bolduc2013exact}.

Fig.~\ref{fig:holograms} illustrates a few examples of the resulting holograms. When computing the hologram, we compensate for the incoming beam's non-uniform profile and the SLM backplane curvature using the procedure described in Ref.~\cite{pushkina2020slm}. The holograms can be modified further to  alter the waist, displacement and orientation of the generated modes.

\begin{figure}[h]
	\includegraphics[width=0.8\columnwidth]{figures/holograms.png}
	\caption{Examples of SLM holograms used to produce HG modes.
		(A) $HG_{0,0}$, (B) $HG_{0,1}$, (C) $HG_{1,0}$,  (D) $HG_{5,5}$, (E) $HG_{10,10}$, (F) $HG_{20,20}$.
	}
	\label{fig:holograms}
\end{figure}

\subsection{Digital micromirror device and generated sources}
\label{MM:DMD}
\indent
The DMD consists of 1920$\times$1080 micromirrors. Each of them can be set to the ON or OFF state, corresponding, respectively, to a tilt by $+12^\circ$ or $-12^\circ$ along their diagonal axis and relative to their flat position. The signal beam illuminates the central area of the DMD, where 210$\times$210 micromirrors display the objects to be imaged.
The DMD micromirrors outside that area are set to the OFF state. The micromirror pitch is 7.56 $\mu \textup{m}$, so the total working area is 1.588 mm wide.

The DMD is driven in ‘pattern on-the-fly’ mode, in which a sequence of 400 binary bitmap images is loaded into the internal memory of the DMD controller board via USB. Each loaded sequence includes 398 objects to be imaged, a phase reference object and an alignment
square, whose purpose is described below.

\begin{figure}[h]
	\includegraphics[width=0.6\columnwidth]{figures/DMD_frames.png}
	\caption{Examples of DMD frames used in the experiment: (A) alignment square; (B) phase reference frame; (C) random matrix (train dataset); (D) superposed random lines (train dataset); (E) random parts of an ellipse (train dataset); (F) pair of lines (test dataset); (G) frame from Oxford University logo (test dataset); (H) text symbol (test dataset).
	}
	\label{fig:DMD_frames}
\end{figure}

\subsection{Acquisition}

We set the acquisition order to minimize the overall measurement time, taking into account the SLM and DMD refresh times ($\sim0.5$ s and $\sim0.1$ ms, respectively), and the time required to load new frames into the DMD internal memory ($\sim 20$ s for a batch of 400 binary frames). After each such batch is loaded into the DMD, it is displayed sequentially while keeping the LO mode constant. The corresponding 400 photocurrents are acquired with a digital oscilloscope in a single trace, along with the AOM driving signal to keep track of the phase. While the trace is transferred from the
oscilloscope to the PC in the binary format via USB (it takes approximately 1 s to transfer 55000 data points and 4 channels), the SLM hologram is updated to produce the next HG mode. The acquisition for one such batch and all 441 LO modes  lasts about 9 minutes. Acquiring the training set ($26501$ DMD frames) takes about 10 hours, while the test set (Oxford logos, text symbols and line pairs) is limited to $302$ DMD frames and is hence acquired as a single batch.

The acquired traces are processed to extract the amplitude and phase of the heterodyne detection photocurrent for each object. The relative phase between the signal and LO arms of the optical setup drifts with time. To keep track of this drift, we use the phase reference object [Fig.~\ref{fig:DMD_frames}(B)] that is the union of the bottom and top-right quadrants of the DMD screen. Its asymmetric shape has nonzero overlap with all LO modes, so this object always generates a measurable photocurrent. The phase $\phi_{j,m,n}$ associated with the LO mode HG$_{m,n}$ and object $j$ can then be calculated according to
\begin{equation}
	\phi_{j,m,n} = \left(\phi_{j,m,n}^{\rm HD} - \phi_{j,m,n}^{\rm AOM} \right) - \left(\phi_{{\rm ref},m,n}^{\rm HD} - \phi_{{\rm ref},m,n}^{\rm AOM}\right),
\end{equation}
where each term is parentheses is the difference between the phases of the photocurrent and the AOM driving signal acquired by the oscilloscope, both oscillating at 92.05 MHz.

The Python code used for the acquisition and pre-processing can be seen in our GitHub repository \cite{HGM_github}.

\subsection{Automatic interferometer alignment procedure}
\indent
Because the data acquisition takes several hours, it is necessary to regularly realign the LO and signal beams with respect to each other in both transverse position and direction. To this end, we have the DMD display at its centre the alignment square of size 100$\times$100 pixels [Fig.~\ref{fig:DMD_frames}(A)]. When the signal and LO are perfectly aligned, we expect null photocurrents when the LO is in modes HG$_{1,0}$ and HG$_{0,1}$. The vertical or horizontal misalignment will give rise to a nonzero photocurrent for the corresponding LO mode.
This error signal is reduced by modifying the displacement and phase gradient of the SLM hologram. The procedure lasts for about 3 minutes and is performed every 3 hours. The corresponding code can be found in Ref.~\cite{HGM_github}.

Note that the alignment square is also displayed in every batch of objects during the regular data acquisition for monitoring purposes.

\subsection{Neural Network}

The NN architecture is shown in Fig.~\ref{fig:NN} and consists of two hidden layers with 6000 units each. Since the normalized input ranges from $-1$ to $1$, we use the hyperbolic tangent as activation function for the first layer. For the second layer, we use a ReLU activation for train speed purposes and, after the third, a sigmoid, in order to adapt to the $[0,1]$ range of the labels. The optimization method is Adaptive Moment Estimation (Adam), with a learning rate of $10^{-4}$, exponential decay rate moving parameters for the first and second moment estimates $(0.9, 0.999)$ and weight decay $0$.
We set the batch size to $512$. We notice that a large batch size, apart from reducing the training time, allows achieving lower training and cross-validation losses.
The loss function is the mean squared error (MSE) loss, which makes the NN effectively behave as a nonlinear regressor. The NN is trained for 900 epochs, achieving training and validation log$_{10}$-loss of $-6.674$ and $-6.498$ per sample, respectively, as illustrated in the inset of Fig.~\ref{fig:NN}.

\subsection{Train, validation and test datasets}
The training/cross-validation dataset consists of
\begin{itemize}
	\item 20000 frames that are tiled with black and white rectangles. The tile size is randomly chosen from 10 to 50 pixels in each dimension and the ``color" of each tile is also random, sampled from a uniform distribution between 0.2 and 0.8; the entire matrix is then randomly rolled along both axes and randomly oriented, Fig.~\ref{fig:DMD_frames}C;
	\item 3000 random sets of 1 to 5 lines with random orientation, position and width ($10$ to $70$ pixel), Fig.~\ref{fig:DMD_frames}D;
	\item 3000 random sets of 1 to 5 segments of elliptic rings with randomly variable eccentricity and width (minor and major axis from $40$ to $120$ pixel, width from $10$ to $15$ pixel), Fig.~\ref{fig:DMD_frames}E;
	\item 58 single squares of sizes $30$ pixel and $70$ pixel at various positions within the  $210 \times 210$ frame;
	\item 441 squares of $30$ pixel size at variable positions within the $210 \times 210$ frame, from which a smaller square of $10$ pixel size is subtracted, with variable positions of the smaller square with respect to the larger one;
	\item a blank image;
	\item a square occupying the entire $210 \times 210$ pixel frame.
\end{itemize}
During the training time, the train/cross-validation dataset is randomly split into 90\% for training and 10\% for validation.

The test dataset consists of
\begin{itemize}
    \item University of Oxford logo, divided into 10$\times$3 frames [Fig.~\ref{fig:reconstructions}(A), top].
    \item University of Oxford coat of arms, divided into 10$\times$12 frames [Fig.~\ref{fig:reconstructions}(A), bottom].
    \item 96 text symbols from the Latin alphabet and special ASCII characters [Fig.~\ref{fig:alphabet_prediction}(A)].
    \item 56 centred pairs of vertical lines of $10$ pixel thickness and varied spacing, used to evaluate the Rayleigh distance.
\end{itemize}

These datasets can be downloaded at \cite{HGM_github}.

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{figures/alphabet_sup.pdf}
	\caption{Original (A), HGM-reconstructed (B) and camera (C) images of the \textit{text symbols} test set.}
	\label{fig:alphabet_prediction}
\end{figure}

\subsection{Direct imaging experiment and deconvolution algorithm}

In order to assess the resolution improvement of HGM with respect to conventional microscopy, we perform direct intensity measurement of the image field with a CMOS camera. For the parallel lines and alphanumerical symbols, a single direct image of the $210\times210$ pixel object is acquired. For the two larger objects, we simulate the scanning as described in the main text.

To further evaluate the performance of HGM, we post-process the experimental camera images using one of the commonly used deconvolution algorithms,
Richardson-Lucy (RL) deconvolution. RL is an iterative algorithm that aims to estimate the original object $O(x,y)$, given the point spread function $P(x,y)$ and the experimental camera image $I(x,y)$, via
\begin{equation}
\hat{O}_{t+1} =\hat{O}_{t} \cdot \Bigg( \frac{I}{\hat{O}_{t} * P} * P \Bigg)
\end{equation}
with $\hat{O}_t$ being the estimate of the object at the  iteration $t$ and asterisks denote convolution.

\begin{figure}[h!]
	\includegraphics[width=1\columnwidth]{figures/deconvolution.png}
	\caption{RL applied to the camera images of the test datasets using the Airy point spread function (A). Estimates of the two University of Oxford logos (B,C) and text symbols (D).}
	\label{fig:RL}
\end{figure}

Fig.~\ref{fig:RL}(B-D) shows the result of RL applied to the test datasets (two University of Oxford logos and the text symbols test), after $1000$ iterations and using the Airy point spread function (PSF) reported in Fig.~\ref{fig:RL}A. The size of the Airy PSF is chosen to match the experimental PSF, with 45 camera pixels from its centre to the radius of the first zero. The program used to perform RL can be found in \cite{HGM_github}.

By qualitatively comparing the HGM predictions (Fig.~\ref{fig:reconstructions} and Fig.~\ref{fig:alphabet_prediction}) with the results of RL (Fig.~\ref{fig:RL}), we conclude that HGM allows for better reconstruction of the original object, with a higher resolution and fewer artifacts.

\subsection{Alternative resolution improvement benchmark}

\begin{figure}
	\includegraphics[width=1\columnwidth]{figures/Full_loss.pdf}
	\caption{Resolution improvement for Oxford logo, Oxford coat of arms and and the text symbols test set, evaluated through the MSE loss as described in the text. For the 20$\times$20 reconstructions, we achieved resolution improvements of $1.76$, $1.78$ and $2.22$, respectively.}
	\label{fig:MSE}
\end{figure}

An alternative way of estimating the resolution improvement is by evaluating the mean squared error (MSE) between the reconstructed images and the original objects. We calculate the MSE $M_{\rm DI}({\rm NA})$ between theoretically calculated direct camera images with varying numerical apertures and the original objects, and then compare these values with the MSE $M_{\rm HGM}({\rm NA}_0)$ for the HGM reconstruction with the numerical aperture ${\rm NA}_0$. Because HGM is generally superior to direct imaging, $M_{\rm DI}({\rm NA}_0)>M_{\rm HGM}({\rm NA}_0)$. However, the function $M_{\rm DI}({\rm NA})$ monotonically decreases with ${\rm NA}$ and becomes equal to $M_{\rm HGM}({\rm NA}_0)$ at some ${\rm NA}_1>{\rm NA}_0$: in other words, the image quality obtained via direct imaging with the numerical aperture ${\rm NA}_1$ is similar to that with HGM and numerical aperture ${\rm NA}_0$. The quantity ${\rm NA}_1/{\rm NA}_0$ corresponds to the resolution enhancement and is plotted, as a function of the number of modes used in HGM, in Fig.~\ref{fig:MSE}.

\rednew{\subsection{Hermite-Gaussian microscopy and shot noise}}

In this section we recapitulate the derivation of the HGM equations for coherent light as in Ref.~\cite{yang2016far} and extend it for two-dimensional imaging. We also discuss how HGM behaves in noisy environments.

Consider a plane object illuminated by a uniform coherent light beam and imaged by an objective lens. At the image plane of the lens, we have the field distribution
\begin{equation}
	E'(x,y) = \int \int_{-\infty}^{+\infty} E(x',y') \psi(x-x',y-y') dx'dy',
\end{equation}
where $(x,y)$ and $(x',y')$ are the image and object plane coordinates, respectively, $E(x,y)$ is the object plane field distribution and $\psi(x,y)$ the point-spread function (PSF) of the lens. This PSF can be well approximated by a ``soft-aperture'' Gaussian PSF $\psi(x,y) = e^{-(x^2+y^2)/4\sigma^2}/\sqrt{2 \pi \sigma^2}$, where $\sigma \approx 0.21 \lambda / \textup{NA}$ \cite{zhang2007gaussian}. Instead of directly detecting the light's intensity $|E'(x,y)|^2$, we perform heterodyne detection using as LO the $(n,m)$-th order HG mode
\begin{equation}
	\psi_{mn}(x,y) = \frac{1}{2^n n!}  H_n\left(\frac{x}{\sqrt{2}\sigma}\right) H_m\left(\frac{y}{\sqrt{2}\sigma}\right)\psi(x,y),
\end{equation}
with the same waist $\sigma$ as the PSF. The resulting heterodyne photocurrents are then proportional to the overlap between the signal and LO beams,
\begin{equation}\label{eq:HGM_J_n}
	J_{mn} \propto \int \int_{-\infty}^{+\infty}E'(x,y)\psi_{mn}(x,y) dx dy = \frac{1}{\sqrt{m! n!}}\int \int_{-\infty}^{+\infty}E(x,y)\left(\frac{x}{2 \sigma}\right)^m \left(\frac{y}{2 \sigma}\right)^n e^{-(x^2+y^2)/8\sigma^2} dx dy.
\end{equation}

The original object plane field can be recovered from the measured photocurrents by calculating the coefficients
\begin{equation}\label{eq:HGM_beta}
	\beta_{mn} = \sum_{k=0}^{m}\sum_{l=0}^{n}\sqrt{m! n!} \alpha_{mk}\alpha_{nl} J_{kl},
\end{equation}
where the $\alpha_{mk}$ are the coefficients of Hermite polynomials, $H_m(x) = \sum_{k=0}^{m} \alpha_{mk} x^k$, and then reconstruct the object field $E(x,y)$ via
\begin{equation}\label{eq:HGM_reconstruction}
	E_{\textup{rec}}(x,y) = \sum_{m,n=0}^{\infty}  \frac{\beta_{mn} H_m\left(\frac{x}{2\sigma}\right)H_n\left(\frac{y}{2\sigma}\right)e^{-(x^2+y^2)/8\sigma^2}}{2^{m+n+2}m!n!\pi \sigma^2}.
\end{equation}

\begin{figure}
	\includegraphics[width=1\columnwidth]{figures/new_HGM.pdf}
	\caption{Oxford logo reconstruction using the closed-form expression \eqref{eq:HGM_reconstruction} in the presence of shot noise. Simulations for two detection methods, different sizes of HG mode basis and different number of photons available per image section (DMD frame) are displayed. When fewer photons are available, the relative shot noise is more significant and reasonable reconstruction is possible only with smaller HG bases. }
	\label{fig:new_HGM}
\end{figure}

These equations are derived in a noiseless, ideal scenario. In a real experiment, even if all sources of systematic errors are eliminated, the unavoidable shot noise will still be present. As can be seen from Eq.~\eqref{eq:HGM_J_n}, the photocurrents get exponentially smaller with the HG mode order, and therefore become heavily affected by shot noise. These small photocurrents, however, are  amplified via the factorial terms in Eq.~\eqref{eq:HGM_beta}. The result is that even the slightest noise leads to a complete breakdown of the reconstruction in Eq.~\eqref{eq:HGM_reconstruction}

In Fig.~\ref{fig:new_HGM} we simulated an idealized scenario where there are no experimental imperfections and the only source of error is shot noise (which is the inverse square root of the number of photons available). The results show that $10^{17}$--$10^{18}$ photons per image section are necessary for HGM to work when HG modes up to (20,20) modes are used, which is an enormous photon overhead. Even when limiting the mode order to  (5,5), still $10^4$--$10^5$ photons are necessary to produce any meaningful result.

Therefore, practical HGM requires an algorithm that regularizes the noise from the higher-order photocurrents while at the same time being able to extract information from them. In this work, this task is addressed by the neural network, and in Fig.~\ref{fig:new_predictions} we assessed its efficiency against shot noise. We observe that although the reconstruction quality degrades with the falling number of photons, we are still able to obtain meaningful results even with a very limited number of photons. Remarkably, HGM with as low as $10^3$ photons per image section yields the same reconstruction quality as direct imaging in the absence of the shot noise.

In Fig.~\ref{fig:new_predictions}, we also analyze how the reconstruction quality is affected when the NN trained on noisy data with certain  number of photons available (and hence a certain shot noise level) is used to reconstruct an image from data acquired with a different number of photons. For a given number of photons in the test set, the optimal reconstruction quality is obtained when training set has been acquired with the same number of photons (i.e. has the same noise level). If difference is present, the reconstruction quality is affected less significantly if the test set has a better signal-to-noise ratio than the training set, rather than the other way around.

\begin{figure}
	\includegraphics[width=1\columnwidth]{figures/loss_photon_n.pdf}
	\caption{Simulated Oxford logo reconstruction via direct imaging and HGM using heterodyne detection with different photon numbers available per image section during training and test phases. a) Examples. b) (Logarithm of) mean squared error between the noisy reconstructed logo and the original object. Horizontal lines on the right correspond to the error obtained for the direct imaging reconstructions with the same photon numbers.}
	\label{fig:new_predictions}
\end{figure}

\end{document}