# Controlling solid-liquid interfacial energy anisotropy through the isotropic
liquid

Lei Wang<EMAIL_ADDRESS><EMAIL_ADDRESS>Jeffrey J. Hoyt  Nan Wang  Nikolas
Provatas  Chad W. Sinclair  Department of Materials Engineering, The
University of British Columbia, 309-6350 Stores Road, Vancouver, British
Columbia, V6T 1Z4, Canada  Key Laboratory of Space Applied Physics and
Chemistry, Ministry of Education, School of Science, Northwestern
Polytechnical University, Xi’an, 710072, China  Department of Materials
Science and Engineering, McMaster University, Hamilton, ON, L8S-4L7, Canada
Department of Physics and Centre for the Physics of Materials, McGill
University, Montreal, H3A 2T8 ,Canada

###### Abstract

Although the anisotropy of the solid-liquid interfacial free energy for most
alloy systems is very small, it plays a crucial role in the growth rate,
morphology and crystallographic growth direction of dendrites. Previous work
posited a dendrite orientation transition via compositional additions. In this
work we examine experimentally the change in dendrite growth behaviour in the
Al-Sm (Samarium) system as a function of solute concentration and study its
interfacial properties using molecular dynamics simulations. We observe a
dendrite growth direction which changes from $\langle 100\rangle$ to $\langle
110\rangle$ as Sm content increases. The observed change in dendrite
orientation is consistent with the simulation results for the variation of the
interfacial free energy anisotropy and thus provides definitive confirmation
of conjecture in previous works. In addition, our results provide physical
insight into the atomic structural origin of the concentration dependent
anisotropy, and deepens our fundamental understanding of solid-liquid
interfaces in binary alloys.

††journal: Nature Communications

## Introduction

Though interfaces constitute a vanishingly small volume fraction of bulk
materials, they play an essential role in determining bulk properties. Among
the myriad of ways that interfaces impact on properties, one of the most
important is their use to control microstructures resulting from phase
transitions. In solidification, it is the intrinsic properties of the solid-
liquid interface that determines the morphology of the selected product phase
and the composition distribution. The interfacial free energy also determines
the characteristic scale and morphology of the microstructure of the solid.
This connection arises directly from the tendency for pattern formation via
dendritic growth [1, 2] owing to the solid dendrite arms wanting to grow
parallel to crystallographic directions whose interfacial energy is highest.
The practical significance of such dendritic growth and its consequences for
the properties of bulk products are broad, from the strength of additively
manufactured parts [3] to the performance of metal-ion batteries [4]. The
significance of this has driven an enormous body of work in the past 10-15
years aimed at identifying the key structural parameters that correlate to the
properties of interfaces at the atomistic scale [5].

According to the generally accepted microscopic solvability theory [6, 7],
dendritic growth velocity and tip radius are sensitively controlled by the
(small) anisotropy in the solid-liquid interfacial free energy, $\gamma$.
Moreover, the anisotropy also determines the crystallographic growth direction
of the dendrites. Pure metals of cubic symmetry will exhibit dendrite arms
extending along the $\langle 100\rangle$ directions so as to minimize the
content of high energy (100) surfaces. Due to its importance, experimental and
modelling techniques have been devised to estimate the anisotropy of $\gamma$
in several systems [8], yet the atomistic origin of this anisotropy remains
elusive [9, 10, 11].

With well-controlled experiments and simulations, Haxhimali et al. show that
the dendrite growth direction in Al-Zn (Zinc) alloys changes with an increase
in the solute content [12]. Starting from normal $\langle 100\rangle$
dendrites in Al-rich alloys, a transition to hyper-branched dendrites and
further to $\langle 110\rangle$ dendrites happens when the content of Zn is
increased. The hyper-branched structure is characterized by no preferred
growth direction while $\langle 110\rangle$ dendrites have distinctly
crystallographic morphologies but with arms oriented parallel to $\langle
110\rangle$. The authors showed that the observed dendrite orientation
transition (DOT) could be explained by assuming that the solid-liquid
anisotropy changes in some fashion with Zn concentration.

Two parameters are often used to characterize the anisotropy of $\gamma$ when
the solid phase has underlying cubic symmetry. The first parameter
($\epsilon_{1}$) represents a four fold symmetric contribution, the second
($\epsilon_{2}$) a six fold contribution. Using phase-field simulations,
Haxhimali et al. found that the space spanned by various values of
$\left(\epsilon_{1},-\epsilon_{2}\right)$ can be separated into three regions
where $\langle 100\rangle$, hyper-branched, or $\langle 110\rangle$ dendrites
dominate the microstructure. Compiling existing estimates of $\epsilon_{1}$
and $\epsilon_{2}$ for a variety of pure cubic metals, Haxhimali et al.[12]
found that most lie inside the $\langle 100\rangle$ region or near the
$\langle 100\rangle$/hyper-branched boundary. This is consistent with the
dominance of $\langle 100\rangle$ dendrites observed experimentally. To
explain the DOT in Al-Zn it was speculated that the anisotropy must vary with
Zn concentration such that the combination of $\epsilon_{1}$ and
$\epsilon_{2}$ move from the $\langle 100\rangle$ to $\langle 110\rangle$
regime. However, it must be stressed that there is no direct evidence why such
an anisotropy variation exists in Al-Zn, nor is there any explanation as to
why Al-Zn is unique with respect to this solid-liquid interfacial free energy
behaviour.

Subsequent molecular dynamics (MD) simulations revealed that the $\gamma$
anisotropy can indeed vary with composition within the
($\epsilon_{1},-\epsilon_{2}$) space. However, for systems examined so far,
such as hard-spheres [13], Lennard-Jones [14] and Cu-Ni [15], there is no
example where the system can vary sufficiently to predict a complete
transition from $\langle 100\rangle$ to $\langle 110\rangle$, which is in
contrast with Al-Zn experiments. This illustrates that the DOT phenomenon, and
the associated anisotropy variations are quite unusual. Interfacial properties
are intrinsically determined by the atomic-scale structure of interfaces. For
solid-liquid interfaces, a well-accepted feature is that the liquid extending
a few interatomic distances from the the crystal exhibits partial ordering.
This interfacial ordering has been investigated extensively (see review [9]
and references therein) and are related to different interface-control
phenomena including nucleation and glass formation. But even though liquid
ordering is generally observed, the DOT has only been reported for several
systems. This indicates that the DOT must rely on yet undescribed features of
interfaces.

In this work, we demonstrate from laser melting experiments on the Al-Sm
system that a complete DOT takes place with increasing concentration of Sm.
Through MD simulations, we verify that its anisotropy behaviour is consistent
with the Haxhimali et al. conjecture of the observed dendrite morphology [12].
We also provide evidence of an unusual atomistic structure of the solid-liquid
interface that may provide insights to the origin of the DOT.

## Results

Figure 1: Dendrite orientation transition in Al-Sm alloys. a, A columnar
dendritic microstructure obtained in a laser remelted Al-2.2 at.% Sm alloy.
The microstructure shown here is consistent with a seaweed structure, closely
related to a hyper-branched dendritic structure described by Haxihimali et al.
[12]. The scale bar indicates 10 $\mathrm{\mu m}$. b, Equiaxed $\langle
110\rangle$ dendrites observed in a laser remelted Al-5.6 at.% Sm alloy (scale
bar 2 $\mathrm{\mu m}$). c, A three dimensional dendrite from a phase-field
simulation, with a section showing the highlighted six $\langle 110\rangle$
branches. The microstructure in (b) is composed of $\langle 110\rangle$
dendrites is verified by the observed 6-fold symmetry.

### DOT in Al-Sm alloys

Figure 1 presents two abnormal microstructures observed in laser melted Al-Sm
(Samarium) alloys. Figure 1(a) shows a columnar seaweed microstructure
obtained from a laser remelted Al-2.2 at.% Sm alloy. The key feature of a
seaweed microstructure is that it exhibits no shape-preserving steady-state
tip morphology [16, 17]. Such a seaweed microstructure is similar to the
hyper-branched structure discussed by Haxihimali _et al._ [12, 18]. Figure
1(b) shows equiaxed dendrites obtained from a laser remelted sample of Al-5.6
at.% Sm. The equiaxed dendrites shown in Figure 1(b) notably display six
branches in one section implying that these cannot be $\langle 100\rangle$
dendrites since a $\langle 100\rangle$ dendrite can have six secondary arms
but those can be only simultaneously viewed in 3D, for example see Figure 2 in
Ref [12]. Six branches in one plane can however be observed with $\langle
110\rangle$ dendrites. As a comparison, a three dimensional $\langle
110\rangle$ dendrite calculated via phase-field simulation is shown in Figure
1(c) to illustrate how six branches can lie in one section. Since dendrites
grow from pure Al melts are $\langle 100\rangle$ oriented, Figure 1 provides
compelling evidence of a transition with increasing Sm content in the liquid,
from $\langle 100\rangle$ to a hyper-branched/seaweed structure at
approximately$\sim$2.2 at.% Sm, and to $\langle 110\rangle$ at $\sim$5.6 at.%
Sm. To make this interpretation fully quantitative, modern high resolution
diffraction techniques could be employed in future studies. It should also be
mentioned that the small size of the observed dendrites is due to the small
excluded volume available. This may also limit the time/space available for
the final arm-selection mechanism to operate. While morphologies can be
influenced when dendrites are confined in space/time, quantitative phase field
simulations do not provide evidence of a significant change in the number or
orientation of dendrite side-arms once the initial crystallographic pattern
emerges.

We attribute the above DOT to solute-induced changes in interfacial free
energy anisotropy rather than to possible solute effects on the interfacial
kinetics. Although laser melting is often taken as an example of rapid
solidification, the solidification conditions vary with position inside the
melt pool. The local solidification velocity, $v_{\mathrm{n}}$, is related to
the laser beam scanning velocity, $v_{\mathrm{B}}$, by
$v_{n}=\mathbf{v}_{\mathrm{B}}\cdot\mathbf{n}$ where $\mathbf{n}$ is the
normal vector to the melt pool boundary after the remelting process reached
the steady state. From the bottom of melt pool $v_{\mathrm{n}}$ increases from
zero to the possible order of $v_{\mathrm{B}}$ at the top surface, depending
on the specific shape [19]. In our case, the seaweeds and $\langle 110\rangle$
dendrites are observed near the bottom of melt pools. For the seaweed
microstructure, we can estimate the growth velocity to be roughly the order of
$\mathrm{mm\ s^{-1}}$ since the $v_{\mathrm{B}}=10\ \mathrm{mm\ s^{-1}}$. For
equiaxed dendrites, the controlling parameter is the undercooling which can be
estimated using the columnar to equiaxed transition (CET) model [20].
Specifically, the equiaxed grains under constrained growth conditions form by
heterogeneous nucleation driven by the constitutional undercooling ahead of
the columnar grains. The undercooling of columnar fronts can be approximated
as the bath undercooling of equiaxed grains. Unfortunately, the equiaxed
grains near the melt pool boundary as in our case (also observed in [21])
cannot be explained satisfactorily by CET theory since the local
solidification velocity and thus undercooling approaches zero there. We
believe those equiaxed grains are due to nucleation of pre-existing inoculants
with very high potency, perhaps from the first scan. Under these speculated
conditions, even a negligible undercooling should be enough to trigger
nucleation and growth near the melt pool bottom/boundary although it is
difficult for us to put an accurate value on this undercooling based on
available data.

At first sight, a direct comparison between experiments under constrained
growth conditions and phase-field simulations in free growth conditions is not
obvious. We emphasize that as long as kinetic effects on anisotropy are not
too strong the transition is expected to happen whether the growth condition
is constrained or not. That is the same strategy used by Haxhimali et al.
[12]. The main purpose of phase-field simulations is to illustrate how
dendritic growth direction changes with the anisotropy and to illustrate six
side-arms in 2D is only possible for $\langle 110\rangle$ dendrites. A
quantitative comparison between simulations and experiments in terms of
supersaturation and the operating state of dendrite is not attempted here.
From previous work [22, 23], it is expected that at a speed of
$\sim\mathrm{mm\ s^{-1}}$, kinetic effects (although present) are not expected
to change the selection mechanisms of dendritic arms for atomically rough
interfaces. If the growth rates become sufficiently high such that the
interface deviates from local equilibrium, then one needs to consider other
factors (including latent heat production and dissipation and the kinetic
effect) into a discussion. For the purposes of the following, we deal only
with conditions where the local equilibrium approximation is appropriate such
that the anisotropy of the interfacial free energy contributes to the dendrite
arm orientation primarily.

Figure 2:  Al-Sm phase diagram. The phase boundaries between liquid and FCC-Al
are calculated from semi-grand canonical MC simulations (squares) using the
Mendelev interatomic potential [24] and compared to the CALPHAD assessment of
Zhou et al. [27] (solid lines). Note that both assessments report $\leq$ 100
ppm solubility of Sm in FCC-Al.

### Solute-dependent interfacial anisotropy

To better understand the atomistic origins of the DOT in Al-Sm alloys, we have
employed molecular dynamics (MD) simulations. An empirical Finnis-Sinclair
type potential developed by Mendelev et al. [24], obtained by fitting the ab-
initio calculated properties of a range of Al-Sm compounds with compositions
close to Al90Sm10, was used here. Notably this potential included fitting of
the structure of liquid Al90Sm10 in its construction. The Al-rich equilibrium
phase diagram predicted by this potential for equilibrium between the liquid
and Al rich FCC phase was computed using semi-grand canonical Monte Carlo
simulations (see [25, 26] and the Methodology Section for details of the
method employed). A comparison between the CALPHAD (CALculation of PHAse
Diagrams) computed [27] Al-Sm phase diagram and that obtained from the
Mendelev interatomic potential are shown in Figure 2. As can be seen, the
agreement between the two estimations are good, though the phase diagram
obtained from the Mendelev potential exhibits a lower melting point for pure
Al and thus a lower liquidus line. Most notable for our purposes is that the
Mendelev potential reproduces the negligible solubility of Sm in FCC-Al solid
solution.

From the results presented in Figure 2 it was possible to construct
equilibrium solid-liquid interfaces at various temperatures/alloy compositions
in MD simulations. Equilibrated solid-liquid interfaces were obtained with
(100), and (110) crystal orientations. Once equilibrated, the capillary
fluctuation method (CFM) [28, 29] was used to obtain the interfacial
stiffness, from which $\epsilon_{1}$ and $\epsilon_{2}$ can be calculated (see
Methods) [12]. Figure 3 illustrates the results of calculated
$\left(\epsilon_{1},-\epsilon_{2}\right)$ values made for interfaces held at
temperatures between 800 K and 900 K in 25 K increments. Overlaid on this plot
are the expected regimes where $\langle 100\rangle$, hyper-branched and
$\langle 110\rangle$ oriented dendrites would be found [12]. Within the
uncertainty of the calculated values of $\epsilon_{1}$ and $\epsilon_{2}$ one
can see a clear variation in anisotropy predicted for the Al-Sm alloys, namely
$\langle 100\rangle$ oriented dendrites predicted for low Sm contents ($<$ 4
at.% Sm, T $>$ 875 K), $\langle 110\rangle$ oriented dendrites for high Sm
contents ($>$ 8 at.% Sm, T $<$ 800 K) and hyper-branched microstructures
predicted in between.

Figure 3: Composition-dependent anisotropy parameters $\epsilon_{1}$ and
$\epsilon_{2}$. The boundaries (dashed lines) separating different dendritic
growth regimes are reproduced from [12]. The coloured points indicate the
specific temperatures at which equilibrium MD simulations were performed. As
the temperature drops, the Sm composition in the liquid phase increases along
the liquidus line while the solid phase remains pure Al as shown by the phase
diagram in Figure 2. The liquid composition increases from $\sim$1.7 at.% at
900 K to $\sim$3.9 at.% (875 K) to $\sim$6.1 at.% (850 K) to $\sim$7.4 at.%
(825 K) to $\sim$8.3 at.% at 800 K. As the Sm content in the liquid increases
(temperature drops), the interfacial anisotropy changes such that the expected
dendritic growth mode shifts from $\langle 100\rangle$ to hyper-branched to
$\langle 110\rangle$. The error bars is calculated based on propagation of
stiffness uncertainty, as described in Methods.

Both our experiments and atomistic simulations show the full transition from
$\langle 100\rangle$ to $\langle 110\rangle$ in Al-Sm alloys. Therefore, the
above results represent the definitive confirmation of the Haxhimali et al.
DOT conjecture [12]. As discussed above, several MD simulations on various
binary alloy systems have shown a variation in anisotropy within the
($\epsilon_{1},-\epsilon_{2}$) space, but no previous case has explicitly
demonstrated a change from $\langle 100\rangle$ oriented dendrites to the
$\langle 110\rangle$ oriented dendrites. Ideally, one would like to compare
with experimentally measured anisotropy parameters. In principle,
($\epsilon_{1},-\epsilon_{2}$) can be derived from the Wulff shape by
measuring the shape of a liquid droplet embedded in a solid matrix. As the
anisotropies of solid-liquid interfaces are extremely weak, experimental
attempts have been only very limited [30, 31]. In addition, the results either
suffer from experimental uncertainty or are not consistent with the observed
dendrite growth directions.

Figure 4: Composition-dependent anisotropy strength and interfacial free
energy. The data used are the same as those in Figure 3. When plotted this
way, it shows that the effect of Sm is to drastically reduce the anisotropy
strength, meaning that the preference for $\langle 110\rangle$ vs $\langle
100\rangle$ dendrites discussed in relation to Figure 3 is small. The error
bars is calculated based on propagation of stiffness uncertainty.

### Atomic structures of interfaces

While the above results give clear evidence that changing the Sm content of
the liquid can change the anisotropy of the interfacial free energy, the
mechanistic origins of this change remain unclear. Another way to view the
data shown in Figure 3 is directly in terms of the difference in interfacial
free energy between the $\langle 100\rangle$ and $\langle 110\rangle$
directions (Figure 4). As can be seen, the anisotropy strength (left axis)
decreases while the average interfacial energy (right ordinate) increases with
Sm content in the liquid. Both of these trends are unusual in alloy systems
and the increase of $\gamma_{0}$ with decreasing temperature runs counter to
the well known result of Spaepen valid for pure metals [10]. However, it is
worth mentioning that in a recent MD study by Hoyt et al. [32], the results
shown in Figure 4 are also observed in Cu-Zr alloys. It is also noteworthy
that both Al-Sm and Cu-Zr exhibit a very small solid solubility and a large
size mismatch. Although it is unclear how these factors directly affect the
interfacial energy and its anisotropy, we discuss below how the presence of Sm
in Al-Sm alloys creates an unusual solid-liquid interfacial structure, which
may provide a possible atomistic origin of the DOT.

Figure 5(a) shows the fine grained number density as a function of distance
perpendicular to the liquid-solid interface for both the (100) and (110)
interfaces at the highest (900 K) and lowest (800 K) temperatures studied.
Both Al atomic density and the atomic density of the Sm (orange) are shown in
the plots. The density profile of Sm is multiplied by a factor of 10 to more
clearly compare with the Al density. Focusing on the Al density we see the
expected behaviour, with partial ordering extending 3-4 layers into the liquid
with a characteristic spacing between peaks similar to that of the solid (see
e.g. [13]). However, a striking feature of the density profiles in Figure 5(a)
is the peaks of the Sm profile do not align with the Al density peaks. To the
best of our knowledge this misalignment differs from any previous MD study of
the structure of solid-liquid interfaces in other systems, see e.g., Cu-Ni
[15] and Lennard-Jones [14]. A solute density profile that is displaced from
the total density (dominated by Al) will result in excess number of solute
atoms and excess interfacial energy contributions that are not found in the
more commonly observed aligned structure and these additional excess
quantities will affect the interfacial free energy [32]. From the point of
view of the DOT, however, a more significant observation is the fact that the
displacement between profiles appears to be greater for the (110) orientation
than for the (100) interface for all temperatures studied. Therefore, this
unique atomistic structure of the Al-Sm solid-liquid interface provides an
additional source of anisotropy, which may play a crucial role in identifying
those alloys likely to undergo a DOT.

A possible source of the temperature dependence of the interfacial energy and
anisotropy is provided in Figure 5(b). These panels show the total in-plane
ordering obtained by time averaging the atomic density in a layer $\sim$3
$\AA$ wide (grey shaded box in Figure 5(a) centered on the interface
position). The insets to Figure 5(b) show the corresponding in-plane density
for the Sm atoms alone. At the higher temperature (900 K), where the Sm
concentration is low, there is little in-plane disordering due to the absence
of Sm. On the other hand, at a temperature of 800 K, shown by the right two
panels, there is considerable disorder induced by the Sm atoms. The
disordering is more pronounced for the (100) interface, a result consistent
with the observation discussed above that the two density profiles are more
closely aligned for the case of (100). Since the disordered planes will
contribute to the excess energy of the interface, the misaligned atomic
structure seen in Al-Sm will provide a contribution, not found in most other
alloys, to the temperature dependence of $\gamma$ and its anisotropy.

Figure 5: Ordering near the solid-liquid interfaces. a, Time averaged atomic
density profiles calculated from layers parallel to the solid-liquid interface
for the highest (900 K) and lowest (800 K) temperatures studied. The solid-
liquid interface has been set at $y=0$. The green curve shows the solvent (Al)
atomic density while the orange curve shows the average density obtained from
Sm only (the Sm density profile is multiplied by 10 so as to plot on the same
scale). b, The total (Al+Sm) in-plane average atomic density. The averages are
performed in $y$ direction over the regions identified by the shaded grey
boxes in (a) (extending $\sim$3 $\AA$ in length perpendicular to the
interface). The inset to each figure shows the average atomic density obtained
from the Sm atoms alone. The in-plane ordering can illustrate a structure
commensurate with the underlying solid. All the profiles are averaged over 500
snapshots (0.5 ns).

## Discussion

The above observations suggest that in alloys which exhibit an interface
structure with misaligned total density (dominated by the solvent density) and
the solute density peaks may be candidates for identifying a DOT. It is of
course difficult to assess whether the same structure observed in Figure 5 is
also present in the Al-Zn system studied by Haxhimali et al., but it is worth
noting that in Al-Zn the equilibrium crystal structure of Zn is HCP and not
FCC. Therefore, one can imagine that at high enough Zn content an out of phase
density structure may be possible, at least in some crystallographic
directions. In addition to the structural properties of the solid-liquid
interface, the results of this study offer other insights into the DOT and
hence the variation in interface anisotropy. Because the crystal induces local
ordering in the liquid at the interface, it is tempting to conclude that any
variation in anisotropy arises from changes in the crystalline structure with
composition (e.g. lattice parameter). But in the case of Al-Sm there is no
solubility in the solid and the origin of the anisotropy is necessarily more
complicated. Also, while the DOT’s in Al-Sm and Al-Zn alloys appear
superficially similar, it must be recalled that Zn, unlike Sm, can be
dissolved significantly (up to 67 at.%) into solid Al. This suggests that the
solution thermodynamics of the bulk liquid and solid phases do not play a
significant role in the explanation of the DOT behaviour.

The results presented above open new challenges and opportunities for
furthering our understanding of the role of the liquid-solid interface
structure and composition on interfacial free energy and consequently
solidification microstructure evolution. If the above interpretation is
correct, and alloying elements that disrupt the partial ordering in the liquid
at the liquid-solid interface can have a strong effect on both the absolute
magnitude and anisotropy of the interfacial free energy, it then suggests that
efforts aimed at quantitative prediction of this effect would be desirable.
Even without a fully quantitative model, however, this correlation suggests a
route to identify alloying elements that would be suitable for inducing a DOT.
This in turn would pave the way for additional alloy design possibilities.
That is, solute additions whose main role is to alter the dendrite morphology
and produced more desirable properties in the solidified material. More
broadly, these results also point to the role that solute could play in
modifying the partial ordering of the liquid at the liquid-solid interface
[33, 34, 35]. Aside from changing interfacial free energy anisotropy, this
ordering has been cited as an important ingredient for promoting or
suppressing crystallization [36, 37], this being an important consideration
both for grain refinement in commercial alloys and in the production of bulk
metallic glasses.

## Methods

### Experiments

The laser remelting experiments presented in Figure 1 were performed on
material prepared from master alloys prepared as part of a previous study
[38]. These alloys were prepared by chill casting 99.999 wt.% pure Al and
99.99 wt.% Sm in a cylindrical mold water cooled from the bottom. The cast
ingots were then cut into plates of thickness of 3 mm. The laser remelting
experiments were conducted with a laser operated at 1280 W using a scanning
speed of 10 $\mathrm{mm~{}s^{-1}}$ for the Al-2.2 at.% Sm (11 wt.% Sm) alloys.
Due to the formation of large primary intermetallic compounds during chill
casting in the Al-5.6 at.% Sm (25 wt.% Sm) alloy, a double scan strategy was
used. A first scan at a low speed (21 $\mathrm{mm~{}s^{-1}}$) was used to
remelt the intermetallic compound and to generate a fine microstructure. The
second scan was them performed inside the remelted region of the first scan at
a high scan speed (667 $\mathrm{mm~{}s^{-1}}$).

### Phase-field simulations

The quantitative phase-field model presented in Figure 1 was used to simulate
equiaxed dendritic growth [39]. This model has been implemented in a 3D
Adaptive mesh refinement Message Passing Interface C++ code as described by
Greenwood et al. [40]. The simulations were performed using dimensionless
parameters, i.e., $W_{\mathrm{0}}=\tau=1$, where $W_{\mathrm{0}}$ and $\tau$
are the characteristic interface width and the kinetic attachment time
respectively. The interfacial free energy anisotropy to produce the results
shown were taken as $\epsilon_{1}=0,\epsilon_{2}=-0.02$. In addition the
supersaturation was taken to be $\omega=0.55$, solute distribution coefficient
$k=0.1$ and dimensionless diffusion coefficient $D=3$.

### Atomistic simulations

All atomistic simulations (molecular dynamics and Monte Carlo simulations)
were performed using the LAMMPS (Large-scale Atomic/Molecular Massively
Parallel Simulator) software [41]. All calculations were performed using the
semi-empirical potential developed by Mendelev et al. [24] for the properties
of Al-Sm alloys. The liquidus and solidus boundaries shown in Figure 2 were
obtained using the thermodynamic integration technique described by Ramalingam
et al. [42], which computes the relative free energy between solid and liquid
phases. The melting point of pure Al which serves as the reference for all
free energies is determined first. Starting from the melting point (914 K
[24]), the free energy difference between liquid and solid pure Al was found
as a function of temperature by a Gibbs-Helmholtz integration. This requires
knowledge of the enthalpy per atom of pure Al in the solid and liquid phases
at a range of temperatures, which can be determined from MD simulations. The
chemical potential difference between Al and Sm in the liquid and FCC phases
was then found by performing semi-grand canonical Monte Carlo simulations at
the several temperatures between 800 K to 900 K with 25 K as the interval. In
this case, simulations were performed on boxes containing 6912 atoms in either
the liquid or solid FCC state. The chemical potential difference applied to
liquid phase ranges from 2.0-2.5 eV/atom and 1.9-2.1 eV/atom for solid phase
approximately. The systems were equilibrated for a total of 2,360,000 MD time
steps and statistics were generated over runs of at least 1,000,000 steps. The
equilibrium composition for liquid ranges from 0-15 at.% while the solid
composition is of the order of 0.1 at.%. The obtained relationship between
chemical potential difference and composition was then integrated to obtain
the functional form of the Gibbs free energy as function of Sm content in the
liquid and FCC solid phases. From this relationship, the liquidus and solidus
lines were determined from a common tangent construction.

To determine the interfacial free energy, the Capillary Fluctuation Method
(CFM) [28, 29] has been used. This method provides the interface stiffness,
from which interface anisotropy parameters can be computed. For each alloy
composition/temperature considered, simulation cells were set up with (100)
and (110) oriented crystal-melt interfaces containing 64,000 ($20\times
20\times 40$ units cells) and 72,000 atoms ($15\times 20\times 30$ unit
cells), respectively. For (100) interface, only one stiffness can be obtained
due to the four-fold symmetry:

$\langle A(k)\rangle=\frac{k_{\mathrm{B}}T}{aSk^{2}},$ (1)

where $A(k)$ is the Fourier amplitude of the interface height profile at a
given value of the magnitude of the wave vector $k$, $S$ stiffness, $a$ the
interfacial area, $k_{\mathrm{B}}$ the Boltzmann’s constant, and $T$ absolute
temperature. For the (110) interface, the other two stiffnesses can be
obtained:

$\langle
A(k_{\mathrm{x}},k_{\mathrm{y}})\rangle=\frac{k_{\mathrm{B}}T}{a(S_{1\bar{1}0}k_{\mathrm{x}}^{2}+S_{001}k_{\mathrm{y}}^{2})}.$
(2)

From those three stiffnesses, the interfacial free energy as a function of
orientation can be obtained:

$\frac{\gamma(n)}{\gamma_{0}}=1+\epsilon_{1}\left(Q-3/5\right)+\epsilon_{2}\left(3Q+66S-17/7\right),$
(3)

where $Q=n_{\mathrm{x}}^{4}+n_{\mathrm{y}}^{4}+n_{\mathrm{z}}^{4}$ and
$S=n_{\mathrm{x}}^{2}n_{\mathrm{y}}^{2}n_{\mathrm{z}}^{2}$ are two cubic
harmonics. $\gamma_{0}$ is the average value and $\epsilon_{i}(i=1,2)$ are two
adjustable parameters. The error bars are calculated based on the propagation
of stiffness uncertainty $\sigma$, which is estimated according to
$\sigma^{2}=\langle|A(k)|^{2}\rangle^{2}\tau(k)/t_{\mathrm{run}}$ where
$t_{\mathrm{run}}$ is total running time and $\tau(k)$ the relax times.
$\tau(k)$ can be calculated by integrating time correlation functions $\langle
A(k,t)A(-k,0)\rangle^{2}/\langle|A(k)|^{2}\rangle^{2}$ derived from the MD
data [29]. For each orientation, 4 replica were used from which to gather
statistics. To achieve two-phase equilibrium, the middle half of the system
was loaded randomly with the equilibrium concentration of the liquid phase,
whereas the two ends of the periodic cell contained the equilibrium
concentration of the solid. The middle portion was then melted at high
temperatures while the remaining atoms were held fixed. Equilibration was
performed for 1 ns in the $\mathrm{NP_{\mathrm{y}}AT}$ ensemble with
$\mathrm{y}$ denoting the interface normal direction. During the equilibration
runs the solute atoms were periodically swapped with solvent atoms according
to a Monte Carlo procedure to achieve different configurations. The frequency
dependence of the averaged power spectrum can then be used to directly obtain
the stiffness [28]. 6,000 snapshots which last 6 ns were collected for each
replica to generate the spectrum. In order to compute the interface position,
a local order parameter was defined as
$\phi=\frac{1}{12}\sum_{i}|r_{i}-r_{\mathrm{fcc}}|^{2}$ where
$r_{\mathrm{fcc}}$ corresponds to ideal fcc positions in the crystal, $r_{i}$
denotes the instantaneous atom position, $i$ refers to the nearest neighbour
atoms and 12 represents the total number of nearest neighbour atoms over which
the summation is taken. The interface location ($y_{\mathrm{loc}}$) and width
($W$) were determined by fitting the local order parameter profile across the
interface to the hyperbolic tangent function,
$c_{1}+c_{2}\tanh(\left(y-y_{\mathrm{loc}}\right)/{W})$.

## Data Availability

The data and codes that support the findings of this study are available from
the corresponding author upon request.

## References

## References

  * [1] Ben-Jacob, E. & Garik, P.  The formation of patterns in non-equilibrium growth.  _Nature_ 343, 523 (1990).
  * [2] Trivedi, R. & Kurz, W.  Dendritic growth.  _International Materials Reviews_ 39, 49–74 (1994).
  * [3] Wang, Y. M. _et al._ Additively manufactured hierarchical stainless steels with high strength and ductility.  _Nature materials_ 17, 63 (2018).
  * [4] Wang, K. _et al._ Dendrite growth in the recharging process of zinc–air batteries.  _Journal of Materials Chemistry A_ 3, 22648–22655 (2015).
  * [5] Mishin, Y., Asta, M. & Li, J.  Atomistic modeling of interfaces and their impact on microstructure and properties.  _Acta Materialia_ 58, 1117–1151 (2010).
  * [6] Amar, M. B. & Brener, E.  Theory of pattern selection in three-dimensional nonaxisymmetric dendritic growth.  _Physical review letters_ 71, 589 (1993).
  * [7] Langer, J. S.  Instabilities and pattern formation in crystal growth.  _Reviews of Modern Physics_ 52, 1 (1980).
  * [8] Hoyt, J. _et al._ Crystal–melt interfaces and solidification morphologies in metals and alloys.  _MRS bulletin_ 29, 935–939 (2004).
  * [9] Kaplan, W. D. & Kauffmann, Y.  Structural order in liquids induced by interfaces with crystals.  _Annu. Rev. Mater. Res._ 36, 1–48 (2006).
  * [10] Spaepen, F.  A structural model for the solid-liquid interface in monatomic systems.  _Acta Metallurgica_ 23, 729–743 (1975).
  * [11] Oh, S. H., Kauffmann, Y., Scheu, C., Kaplan, W. D. & Rühle, M.  Ordered liquid aluminum at the interface with sapphire.  _Science_ 310, 661–663 (2005).
  * [12] Haxhimali, T., Karma, A., Gonzales, F. & Rappaz, M.  Orientation selection in dendritic evolution.  _Nature Materials_ 5, 660–664 (2006).
  * [13] Davidchack, R. L. & Laird, B. B.  Simulation of the hard-sphere crystal–melt interface.  _The Journal of chemical physics_ 108, 9452–9462 (1998).
  * [14] Becker, C. A., Olmsted, D., Asta, M., Hoyt, J. J. & Foiles, S. M.  Atomistic underpinnings for orientation selection in alloy dendritic growth.  _Physical Review Letters_ 98, 125701 (2007).
  * [15] Qi, C., Li, J., Xu, B., Kong, L. & Zhao, S.  Atomistic characterization of solid-liquid interfaces in the Cu-Ni binary alloy system.  _Computational Materials Science_ 125, 72–81 (2016).
  * [16] Amoorezaei, M., Gurevich, S. & Provatas, N.  Orientation selection in solidification patterning.  _Acta Materialia_ 60, 657–663 (2012).
  * [17] Xing, H. _et al._ Degenerate seaweed to tilted dendrite transition and their growth dynamics in directional solidification of non-axially oriented crystals: a phase-field study.  _Sci Rep_ 6, 26625 (2016).
  * [18] Dantzig, J., Di Napoli, P., Friedli, J. & Rappaz, M.  Dendritic growth morphologies in Al-Zn alloys——Part II: phase-field computations.  _Metallurgical and Materials Transactions A_ 44, 5532–5543 (2013).
  * [19] Zimmermann, M., Carrard, M. & Kurz, W.  Rapid solidification of Al-Cu eutectic alloy by laser remelting.  _Acta Metallurgica_ 37, 3305 – 3313 (1989).
  * [20] Hunt, J.  Steady state columnar and equiaxed growth of dendrites and eutectic.  _Materials science and engineering_ 65, 75–83 (1984).
  * [21] Mokedem, S.  _Epitaxial Laser Treatment of Single Crystal Nickel-base Superalloys_.  Ph.D. thesis, Ecole Polytechnique Fédérale de Lausanne, Lausanne (2004).
  * [22] Bragard, J., Karma, A., Lee, Y. H. & Plapp, M.  Linking phase-field and atomistic simulations to model dendritic solidification in highly undercooled melts.  _Interface Science_ 10, 121–136 (2002).
  * [23] Pinomaa, T. & Provatas, N.  Quantitative phase field modeling of solute trapping and continuous growth kinetics in quasi-rapid solidification.  _Acta Materialia_ 168, 167–177 (2019).
  * [24] Mendelev, M. I. _et al._ Development of interatomic potentials appropriate for simulation of devitrification of Al90Sm10 alloy.  _Modelling & Simulation in Materials Science & Engineering_ 23, 045013 (2015).
  * [25] Hoyt, J. J., Asta, M. & Karma, A.  Atomistic simulation methods for computing the kinetic coefficient in solid-liquid systems.  _Interface Science_ 10, 181–189 (2002).
  * [26] Hoyt, J. J., Asta, M. & Karma, A.  Atomistic and continuum modeling of dendritic solidification.  _Materials Science & Engineering R Reports_ 41, 121–163 (2003).
  * [27] Zhou, S. & Napolitano, R. E.  Modeling of thermodynamic properties and phase equilibria for the Al-Sm binary system.  _Metallurgical and Materials Transactions A_ 39, 502–512 (2008).
  * [28] Hoyt, J. J., Asta, M. & Karma, A.  Method for computing the anisotropy of the solid-liquid interfacial free energy.  _Physical Review Letters_ 86, 5530 (2001).
  * [29] Asta, M., Hoyt, J. J. & Karma, A.  Calculation of alloy solid-liquid interfacial free energies from atomic-scale simulations.  _Physical Review B_ 66, 351–363 (2002).
  * [30] Liu, S., Napolitano, R. & Trivedi, R.  Measurement of anisotropy of crystal-melt interfacial energy for a binary Al-Cu alloy.  _Acta Materialia_ 49, 4271 – 4276 (2001).
  * [31] Napolitano, R. E. & Liu, S.  Three-dimensional crystal-melt Wulff-shape and interfacial stiffness in the Al-Sn binary system.  _Physical Review B_ 70, 214103 (2004).
  * [32] Hoyt, J., Raman, S., Ma, N. & Asta, M.  Unusual temperature dependence of the solid-liquid interfacial free energy in the Cu-Zr system.  _Computational Materials Science_ 154, 303 – 308 (2018).
  * [33] Hashibon, A., Adler, J., Finnis, M. W. & Kaplan, W. D.  Ordering at solid-liquid interfaces between dissimilar materials.  _Interface Science_ 9, 175–181 (2001).
  * [34] Hashibon, A., Adler, J., Finnis, M. W. & Kaplan, W. D.  Atomistic study of structural correlations at a liquid–solid interface.  _Computational materials science_ 24, 443–452 (2002).
  * [35] Sun, Y. _et al._ Structural ordering at solid-liquid interfaces in al-sm system: a molecular-dynamics study.  _Materials Letters_ 186, 26–29 (2016).
  * [36] Greer, A. L.  Liquid metals: Supercool order.  _Nature Materials_ 5, 13 (2006).
  * [37] Geysermans, P., Gorse, D. & Pontikis, V.  Molecular dynamics study of the solid–liquid interface.  _The Journal of Chemical Physics_ 113, 6382–6389 (2000).
  * [38] Wang, N., Kalay, Y. E. & Trivedi, R.  Eutectic-to-metallic glass transition in the Al–Sm system.  _Acta Materialia_ 59, 6604–6619 (2011).
  * [39] Karma, A.  Phase-field formulation for quantitative modeling of alloy solidification.  _Physical Review Letters_ 87, 115701 (2001).
  * [40] Greenwood, M. _et al._ Quantitative 3D phase field modelling of solidification using next-generation adaptive mesh refinement.  _Computational Materials Science_ 142, 153–171 (2018).
  * [41] Plimpton, S.  Fast parallel algorithms for short-range molecular dynamics.  _Journal of computational physics_ 117, 1–19 (1995).
  * [42] Ramalingam, H., Asta, M., Van de Walle, A. & Hoyt, J.  Atomic-scale simulation study of equilibrium solute adsorption at alloy solid-liquid interfaces.  _Interface science_ 10, 149–158 (2002).

## Acknowledgement

We thank Compute Canada for computing resources. We also thank Mike Greenwood
who developed, with N. P., the phase-field code used herein. N. P. would like
to thank the Canada Research Chairs (CRC) program and N. P. and C. S. the
National Science and Engineering Research Council of Canada (NSERC) for
funding.

## Authors contribution

The first idea started from discussions about the seaweed formation between L.
W. and N. P. at McGill University. L. W. observed the abnormal dendritic
morphologies in Al-Sm alloys and initiated the collaboration with the other
authors. L. W. and J. J. H calculated the interface energy anisotropy. N. W.
did the laser remelting experiments. L. W. and N. P. conducted the phase-field
simulations. C. S. contributed his expertise in atomic-scale simulations and
interfacial thermodynamics. All authors co-edited the manuscript.

## Competing Interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to L. W.


# Passive superresolution imaging of incoherent objects

Jernej Frank  Department of Physics, University of Oxford, Oxford, OX1 3PU, UK
Alexander Duplinskiy  Department of Physics, University of Oxford, Oxford, OX1
3PU, UK  Kaden Bearne  Department of Physics, University of Oxford, Oxford,
OX1 3PU, UK  A. I. Lvovsky  Department of Physics, University of Oxford,
Oxford, OX1 3PU, UK<EMAIL_ADDRESS>
###### Abstract

We investigate Hermite Gaussian Imaging (HGI) — a novel passive super-
resolution technique — for complex 2D incoherent objects in the sub-Rayleigh
regime. The method consists of measuring the field’s spatial mode components
in the image plane in the overcomplete basis of Hermite-Gaussian modes and
their superpositions and subsequently using a deep neural network to
reconstruct the object from these measurements. We show a three-fold
resolution improvement over direct imaging. Our HGI reconstruction retains its
superiority even if the same neural network is applied to improve the
resolution of direct imaging. This superiority is also preserved in the
presence of shot noise. Our findings are the first step towards passive super-
resolution imaging protocols in fluorescent microscopy and astronomy.

††journal: opticajournal

##  1 Introduction

With the human eye, we can only see objects as small as a thin strand of hair,
which is about 100 microns in width [1]. Modern life and materials sciences
are keen to study structures well beyond this limit, and microscopy is an
indispensable tool for this purpose. The resolution of conventional optical
microscopy is, however, limited due to diffraction on the objective lens to
approximately a half of the optical wavelength ($\sim$200 nm) [2, 3, 4]. In
the last few decades, super-resolution (SR) imaging techniques emerged
focusing on either manipulating the illumination of the sample or using
stochastic properties of fluorescent emitters to achieve resolution well
beyond the diffraction limit [5, 6, 7, 8]. Some of these methods were marked
by the 2014 Nobel Prize in chemistry [9]. Even higher resolution — down to
single nanometers — can be achieved by atomic force microscopy [10] and
electron microscopy [11]. However, all existing superresolution methods
involve active interaction with the sample, which imposes significant
limitations. First, they involve the risk of altering or damaging the samples.
Second, these techniques are not universally applicable to all imaging
scenarios since they rely on specific underlying properties of the samples —
such as, e.g., nonlinear susceptibility. Even if the sample is amenable to a
particular method, SR microscopes require expert knowledge, and various
parameters must be adjusted on a case-by-case basis to render a truthful image
reconstruction.

These limitations can be overcome if the imaging method involves only far-
field, linear-optical, and passive measurement of the optical field arriving
from the object. In 2016 Tsang et al. [12, 13, 14, 15, 16] theoretically
predicted that the diffraction limit can be overcome for estimating the
distance between two incoherent point sources in the far-field regime by
decomposing the field in the image plane into an orthonormal basis of spatial
modes (typically, Hermite-Gaussian) and measuring the intensity of each
component. This approach, dubbed SpaDe for “spatial demultiplexing” was soon
experimentally demonstrated by various groups [17, 18, 19, 20, 21, 22, 23].

On this basis, Yang et al. proposed an extension into reconstructing complete
2D objects dubbed Hermite-Gaussian Imaging (HGI) [17]. Pushkina et al.
experimentally implemented HGI for coherent 2D objects and showed a two-fold
resolution improvement over the diffraction limit [24].

However, the latter experimental work used coherent light to illuminate the
sample and detect the reflected light. On the other hand, fluorescence emitted
by biological samples is incoherent. Also, if SpaDe is to be applied outside
microscopy — say, in astronomy or remote imaging, it needs to be adapted to
incoherent input. Here we address this requirement by extending HGI to imaging
incoherent 2D objects. Since the incoherent light emitted from the object
lacks phase information, it is necessary to measure in an overcomplete spatial
mode basis, involving not only HG modes but also their superpositions. We
measure the intensity in each basis mode using heterodyne detection and then
utilize a deep neural network (DNN) to reconstruct the objects. We observe a
significant resolution improvement over direct imaging (DI) in various
scenarios and certified by quantitative benchmarks.

##  2 Theoretical Background

###  2.1 Image reconstruction theory

We briefly recap the reconstruction method of Yang et al. [17]. The object
being imaged is defined by the incoherent intensity distribution $I(x,y)$. We
reconstruct this distribution by expanding it into the HG basis

$\displaystyle
I_{\text{rec}}(x,y)=\sum_{m,n=0}^{N}\frac{\beta_{mn}\phi_{mn}(x,y)}{2^{m+n+2}m!n!\pi\sigma^{2}},$
(1)

where
$\phi_{mn}(x,y)=H_{m}\left(\frac{x}{\sqrt{2}\sigma}\right)H_{n}\left(\frac{y}{\sqrt{2}\sigma}\right)e^{-\frac{x^{2}+y^{2}}{4\sigma^{2}}}$
are the HG polynomials, $N$ is the cutoff for the number of HG modes to be
used, $\sigma\approx 0.21\lambda/\text{NA}$ is the radius of the PSF at
$e^{-1/2}$ intensity [25], and $\beta_{mn}$ are the expansion coefficients.
The PSF mode is assumed to match $\phi_{00}$.

These coefficients can be obtained experimentally by measuring the powers of
the heterodyne detector photocurrents with the local oscillator prepared in a
set of modes (see Appendix for a full derivation). This set of modes can be
the HG basis [Fig. 1(a)]; however, in this case, we are only able to recover
the $\beta$’s with even indices, which, in turn, allow us to reconstruct only
the even component of the object intensity distribution: $I_{\rm
even}(x,y)=[I(x,y)+(-x,y)+I(x,-y)+I(-x,-y)]/4$, which results in ghost image
artefacts [17, 13]. To address this problem, Yang et al.[17] proposed
displacing the object into a single quadrant of the x-y plane, thereby
separating the image from its “ghosts”. However, the displaced object together
with the “ghosts” is larger in size and requires higher-order modes for its
expansion into the HG basis. Furthermore, this approach is unsuitable for
objects with a wide spatial extent.

Figure 1: (a) Subset of holograms displayed on the SLM to generate the LO in
Hermite-Gaussian modes. (b) Expanded iHG basis, which includes 8 additional
superpositions of neighbouring HG modes to obtain the odd geometric moments
for image reconstruction.

Tsang et al. proposed an alternative strategy: recovering the odd components
of the image by measuring the powers in an overcomplete mode basis known as
interferometric HG (iHG) [13, 26]. The version of the iHG basis used in our
work includes the following superpositions of individual HG modes [Fig. 1(b)]:

$\phi_{mn}^{0-8}:=\\{\phi_{mn},\phi^{\pm}_{m,m+1,n,n},\phi^{\pm}_{m,m,n,n+1},\phi^{\pm}_{m,m+1,n,n+1},\phi^{\pm}_{m,m+1,n,n-1}\\},$

where

$\displaystyle\phi^{\pm}_{m,m^{\prime},n,n^{\prime}}:=\frac{1}{\sqrt{2}}[\phi_{m^{\prime}n^{\prime}}\pm\phi_{mn}].$
(2)

For each mode $\phi_{mn}^{i}$, the power $P_{mn}^{i}$ of the associated
heterodyne photocurrent is measured. These measurements then evaluate the
coefficients $\beta_{mn}$, from which the image is recovered (see Appendix).
Note that theoretically, iHG modes $\phi_{mn}^{i}$ with $i\in\\{0,\ldots,6\\}$
are sufficient for complete reconstruction. However, the experiment uses the
entire iHG mode set with $i\in\\{0,\ldots,8\\}$, as defined above.

###  2.2 Simulations

To test the viability of our approach we first simulate image reconstruction
computationally using (1) and compare it to simulated DI. We set the numerical
aperture (NA) to $5.5\times 10^{-4}$ and $\lambda=$ 785 nm (the same as used
in the experiment).

As our objects, we use symbols of a bitmap ASCII font. Each symbol is a black-
and-white bitmap $8\times 6$ pixels, with the size of each pixel being 0.36
$\sigma$ [Fig. 2(a), upper left panel]. To simulate DI, we convolve the PSF
with the objects [Fig. 2(a), upper middle panel]. The upper right panel of
Fig. 2(a) shows the result of applying the conventional Richardson-Lucy (RL)
deconvolution algorithm [27] to the DI data. Next, we simulate the image
reconstruction for various cut-off numbers $N=2,5,30$ [Fig. 2(a), bottom row].
As can be seen, the DI of most of the letters is unreadable, whether or not RL
is applied. But HGI allows us to achieve a satisfactory level of readability,
the quality improving with $N$.

Figure 2: (a) Letters of the alphabet with varying resolution. Top row: GT is
the ground truth, DI represents the direct image whose resolution is limited
by diffraction, RL is the Richardson-Lucy deconvolution of the DI. Bottom row:
simulated HGI reconstructions in an idealized setting with up to
$m=n=N=2,5,30$ modes. (b) MS-SSIM calculated between (1) the ground truth
(infinite resolution) and (2) the simulated HGI, DI, and RL reconstructions.
MS-SSIM score of $\gtrsim 80\%$ achieves letter readability.

We quantify the image quality with respect to the ground truth objects by
computing the Multi-Scale Structural Similarity (MS-SSIM) [28] [Fig. 2(b)],
which is a popular metric in the image processing community, designed to mimic
how the human eye perceives images [29]. We find that the resolution
improvement is proportional to $N$ and saturates for $N>25$. Comparing the
images with the MS-SSIM graph in Fig. 2, we can see that satisfactory image
quality is achieved at about $80\%$.

##  3 Experiment

###  3.1 Setup

We adopt the experimental setup from Ref. [24] (Fig. 3). We use a continuous
wave diode laser (Eagleyard EYP-DFB-0785), operating at 785 nm. The powers of
the image field in individual HG modes are measured by means of heterodyne
detection with the local oscillator (LO) prepared in the corresponding mode.

Figure 3: Experimental setup to image incoherent objects with HGI via homodyne
detection and DI.

At the outset, we split the laser into the signal LO paths. We use the signal
path to illuminate the object displayed on a digital micromirror device or DMD
(DLP LightCrafter 6500). The DMD pixel pitch is 7.56 $\mu$m, and our objects
are displayed in a square area of $160\times 160$ pixels. This area is divided
into $10\times 10$ “logical pixels", each being a square of $16\times 16$
physical DMD pixels. Our binary objects are bitmaps of these logical pixels.
To simulate incoherent objects, the DMD displays logical pixels one at a time,
and for each pixel, photocurrent powers $P_{mn}^{i}$ for all iHG modes are
measured. In order to determine $P_{mn}^{i}$ corresponding to a given multi-
pixel object, we digitally add the measured values of $P_{mn}^{i}$ for all of
the object’s pixels. This approach is chosen instead of the usual method of
obtaining spatial incoherence via a rotating ground glass because the latter
method would require averaging of the speckle pattern for each object and each
combination $(m,n,i)$, resulting in prohibitively slow data acquisition.

After the light is reflected from the object we collect it by an objective
lens. In front of the objective lens, at a distance 242.5 cm away from the
DMD, we place an iris of diameter 2.66 mm to reduce the NA to $5.5\times
10^{-4}$. This configuration sets the Rayleigh diffraction limit to a spacing
of 7 logical DMD pixels (870 $\mu$m).

To generate a particular iHG mode in the LO path, we display a phase grating
(“hologram”) on a reflective phase-only liquid-crystal-on-silicon SLM
(Hamamatsu X13138-02) and select the first diffraction order of the reflected
light [30, 31]. We switch through holograms corresponding to different iHG
modes sequentially.

The signal and LO are recombined on a beam splitter and sent onto a balanced
detector. By placing a flip mirror after the objective lens we are able to
switch between recording HGI with the heterodyne detector and DI using a high-
resolution CMOS camera.

###  3.2 Deep neural networks

As realized in Ref. [24], the experimental imperfections such as an imperfect
Gaussian PSF, optical aberrations, shot noise, electronic white noise, optical
misalignment, iHG mode fidelity in the LO, imperfect overlap between signal
and LO, and general air fluctuation or vibrations render an image
reconstruction using (1) impossible because the measurement errors of lower
order iHG modes propagate into higher order $\beta_{mn}$ coefficients.

We, therefore, recast our reconstruction problem into a machine-learning task
and let a DNN find an optimal mapping between imperfect measurements and the
underlying basis needed for image reconstruction. In this way, we can benefit
from both the new physics of HGI measurement and machine learning to extract
the image with the highest possible resolution. The neural network takes the
measured photocurrent powers as input and outputs the reconstructed image.

The DNN architecture is shown in Fig. 4(a). We have four fully-connected
hidden layers with 400-800-800-400 neurons. The activation functions are
hyperbolic tangent in the first layer ReLU in the next three layers; the last
layer adapts the [0,1] range by using sigmoid.

Figure 4: (a) Architecture of the DNN used for image reconstruction from iHG
modes. The input layer (blue) dimension varies with the number of iHG modes
utilized. The hidden layers (green) reduce the latent space to 400 neurons.
The output layer (red) is fixed to output images $50\times 50$ pixels. (b)
Examples of bitmaps used for training the DNN.

To acquire the DNN training and cross-validation sets, we obtain the
photocurrent powers in iHG modes $\phi_{mn}^{l}$, with
$m,n\in\\{0,\ldots,20\\}$ and $l\in\\{0,\ldots,8\\}$ for 20000 bitmaps
representing random matrices, ellipses, and lines[Fig. 4(b)]. In order to
investigate the dependence of the reconstruction quality on the number of
modes acquired, several subsets are selected from the above set, corresponding
to $m,n\leq M$ with $M=1,4,9,14,20$. These subsets are split into training and
cross-validation sets in proportion 90:10.

Using the ground truth images for labels would lead to severe overfitting
[24]. We, therefore, produce our labels according to (1). That is, our label
is the image of the object as if it were reconstructed via HGI in the absence
of imperfections. We produced several label sets of varied resolution by
choosing the mode cut-off number $N$ between 5 and 35.

For the DNN training, we set our loss function to be the mean squared error
(MSE) and use a minibatch size of $32$. We use the adaptive moment estimation
(Adam) optimiser, learning rate $10^{-4}$, exponential moving averaging
parameters for the first and second moment estimates $(0.9,0.999)$
respectively, and weight decay $0$. The training is run on an 11th Generation
Intel Core i9-11900K CPU with 64 GB memory and an NVIDIA® RTX™ A5000 GPU with
24 GB memory. Training the DNN up to $2000$ epochs takes 40–60 minutes on
average.

##  4 Results

We evaluate the performance of our method using the alphabet letters [Fig.
2(a)] as well as pairs of parallel lines with varying separation. None of
these objects belongs to the training set. The experimental data for these
objects have been acquired separately from the training set data.

We train several DNNs where we vary both the number $M$ of iHG modes in the
input and the number $N$ of HG modes used in simulating the labels. To show
that the resolution improvement arises due to the new physics of HGI
measurements and not only thanks to machine learning we also trained the same
DNN using DI intensity measurements as inputs. These inputs are $63\times
63=3969$ pixel bitmaps which correspond to roughly the same dimension of the
input layer as for iHG with $M=20$.

We compare the results qualitatively by looking at the image reconstruction
results [Fig. 5(a)] and quantitatively by computing the MS-SSIM between the
ground truth (infinite resolution) and the respective outputs of the DNNs
[Fig. 5(c)].

We see that the resolution improvement is maximized for $N=30$ label modes and
$M=10$ to $15$ measured iHG modes. We attribute the lack of resolution
improvement for higher $M$ to the photocurrents in higher-order iHG modes
being very low and hence prone to noise.

We can also observe some resolution improvement beyond the diffraction limit
with DI purely through DNN software enhancement. However, this enhancement is
smaller than that of HGI. We can see this directly from the reconstructed
images [Fig. 5(a)] and from the MS-SSIM benchmarks [Fig. 5(b)]. We also
observe that training the DNN for DI with the number $N$ of label modes as low
as 10 quickly leads to overfitting, which reduces the reconstruction quality
on the test set [Fig. 5(c)].

To determine the resolution limit, we reconstruct the images for pairs of
parallel lines and measure the intensity profile by integrating over the
vertical dimension [Fig. 6]. We see that the resolution of direct imaging is
limited by the diffraction limit, corresponding to a line separation of 7
logical DMD pixels ($847\,\mu$m). DI+DNN gives some improvement over that
limit by resolving lines 4 logical pixels ($484\,\mu$m) apart, but HGI is able
to resolve the separation of 2 logical pixels ($242\,\mu$m), which compared to
DI+DNN still yields a factor of 2 resolution enhancement.

Figure 5: (a) Alphabet letters imaged using DI (left), DNN-enhanced DI
(middle), and HGI using $M=15$ iHG modes (right). (b) Multi-scale structural
similarity values for different pairs of input and label mode numbers. (c)
Training and cross-validation MSE as a function of training epoch. DI+DNN
overfits and we stop training at 1000 epochs, while for HGI we train for the
complete 2000 epochs.  Figure 6: Reconstructed images of pairs of parallel
lines via direct imaging (top), DNN-enhanced direct imaging (middle) and HGI
(bottom). Coloured frames correspond to the resolution limits in the three
settings.

Lastly, we investigate the effect of simulated shot noise on the performance
of HGI and DI. To this end, we assume that all measurements (iHG mode
intensities in HGI and pixel intensities in DI) are implemented via photon
counting, and the total number of photons detected is $N_{\rm ph}$. Then the
number of photons $\lambda_{i}$ in each measurement is $N_{i}=N_{\rm
ph}\lambda_{i}/\sum_{i}\lambda_{i}$ and the corresponding rms shot noise is
$\Delta N_{i}/N_{i}=\Delta\lambda_{i}/\lambda_{i}=1/\sqrt{N_{i}}$. We simulate
the shot noise by adding a normally distributed random value with the rms
deviation $\Delta\lambda_{i}$ to every $\lambda_{i}$.

We set $N_{\rm ph}$ to 100, 400 and 10,000 and add the photon shot noise
numerically before training the DNNs. As can be seen from Fig. 7, the
advantage of HGI over DI for $N_{\rm ph}=10^{4}$ is largely similar to the
classical case. Particularly strong is the impact of the shot noise on the
loss function of DI+DNN [Fig. 7(c)]: the neural network begins overfitting as
early as at epoch 300.

On the other hand, the edge of HGI over DI+DNN degrades for $N_{\rm ph}=400$
and especially $100$ [Fig. 8]. This is because the total number of photons
available is no longer sufficient to reliably reconstruct the relatively
intricate image even with the help of HGI.

Figure 7: (a) Alphabet letter images with simulated shot noise assuming
$N_{\rm ph}=10000$ photons using DI (left), DNN-enhanced DI (middle), and HGI
using $M=15$ iHG modes (right). (b) MS-SSIM values for different pairs of
inputs and labels. The performance is similar to the classical case. (c)
Training and cross-validation MSE as a function of training epochs for
selected input and label pairs. DI+DNN overfits at 300 epochs, at which point
the training is stopped, while the DNN for HGI is trained for 2000 epochs.
Figure 8: Reconstructed images using HGI (left) and DI+DNN (right) of pairs of
parallel lines with added photon shot noise for $N_{\rm ph}=100,400,10.000$.

##  5 Conclusion

Here we presented the first results in superresolving reconstruction of
complex incoherent 2D objects by completely passive measurements using linear
optics in the far-field regime. We have shown that a combination of
information-rich measurements with a simple DNN produces images with $\sim
3$-fold resolution improvement over DI. We can enhance the DI resolution by
applying the same DNN; however, the resulting resolution is still a factor of
$\sim 2$ below that achieved with HGI. This shows the key role of the new
physics of HGI in resolution enhancement.

Due to its passive nature, HGI can become beneficial in various scenarios far
beyond microscopy. This includes astronomy, remote sensing, satellite imaging
and many other settings where active interaction with the object is
impossible. The main obstacles to the practical application of our scheme are
the training of the DNN required for the reconstruction and the bandwidth
limitations of heterodyne detection. To address the latter concern, the
homodyne detector needs to be replaced by a mode sorter — a device that would
spatially separate individual modes of the incoming light field, after which
the intensities of these modes can be measured by direct photodetection [32].

The training set samples can be fabricated using commercially available
fluorescent beads in various sizes. Alternatively, calibration slides can be
fabricated with a resolution of a few tens of nanometers through lithography
or laser writing. In both cases, the training samples will need to be imaged
employing electron microscopy or another method to obtain ground truth data.
With the rapid advancement of machine learning technology, we can also
envision that better DNN architectures and loss functions will emerge that
would drastically reduce the size of the required training set. In the
meantime, we can expand the training set using known data augmentation
techniques [33].

It would be interesting to combine HGI with existing passive SR microscopy
methods based on engineered sample illumination, such as confocal and
structured-illumination microscopy [34]. Because these methods and HGI rely on
different physics to enhance the resolution, we expect the enhancement effect
achieved through both methods to be cumulative. In this way, we hope passive
microscopy resolution with engineered illumination to eventually reach the
scale of 50 nm.

Funding The project is funded by BBSRC grant BB/X004317/1. JF is supported by
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 956071.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement Data underlying the results presented in this
paper are not publicly available at this time but may be obtained from the
authors upon reasonable request.

Supplemental document See Supplement 1 for supporting content.

## References

  * [1] M. Yanoff and J. S. Duker, _Ophthalmology_ (Elsevier Health Sciences, 2008).
  * [2] L. Rayleigh, “Xxxi. investigations in optics, with special reference to the spectroscope,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 8, 261–274 (1879).
  * [3] E. Abbe, “Beiträge zur theorie des mikroskops und der mikroskopischen wahrnehmung,” Archiv für mikroskopische Anatomie 9, 413–468 (1873).
  * [4] E. Abbe, “Xv.—the relation of aperture and power in the microscope (continued),” Journal of the royal microscopical society 3, 790–812 (1883).
  * [5] S. Pujals, N. Feiner-Gracia, P. Delcanale, I. Voets, and L. Albertazzi, “Super-resolution microscopy as a powerful tool to study complex synthetic materials,” Nature Reviews Chemistry 3, 68–84 (2019).
  * [6] G. Huszka and M. A. Gijs, “Super-resolution optical imaging: A comparison,” Micro and Nano Engineering 2, 7–28 (2019).
  * [7] J. Valli, A. Garcia-Burgos, L. M. Rooney, B. V. d. M. e Oliveira, R. R. Duncan, and C. Rickman, “Seeing beyond the limit: A guide to choosing the right super-resolution microscopy technique,” Journal of Biological Chemistry 297 (2021).
  * [8] S. Liu, P. Hoess, and J. Ries, “Super-resolution microscopy for structural cell biology,” Annual Review of Biophysics 51, 301–326 (2022). PMID: 35119945.
  * [9] L. Möckl, D. C. Lamb, and C. Bräuchle, “Super-resolved fluorescence microscopy: nobel prize in chemistry 2014 for eric betzig, stefan hell, and william e. moerner,” Angewandte Chemie International Edition 53, 13972–13977 (2014).
  * [10] G. Binnig, C. F. Quate, and C. Gerber, “Atomic force microscope,” Physical review letters 56, 930 (1986).
  * [11] D. Mathys, “Die entwicklung der elektronenmikroskopie vom bild über die analyse zum nanolabor,” Zentrum für Mikroskopie, University of Basel p. 8 (2004).
  * [12] M. Tsang, R. Nair, and X.-M. Lu, “Quantum theory of superresolution for two incoherent optical point sources,” Physical Review X 6, 031033 (2016).
  * [13] M. Tsang, “Subdiffraction incoherent optical imaging via spatial-mode demultiplexing,” New Journal of Physics 19, 023054 (2017).
  * [14] M. Tsang, “Resolving starlight: a quantum perspective,” Contemporary Physics 60, 279–298 (2019).
  * [15] M. Tsang and R. Nair, “Resurgence of rayleigh’s curse in the presence of partial coherence: comment,” Optica 6, 400–401 (2019).
  * [16] M. Tsang, “Quantum limit to subdiffraction incoherent optical imaging,” Physical Review A 99, 012305 (2019).
  * [17] F. Yang, A. Tashchilina, E. S. Moiseev, C. Simon, and A. I. Lvovsky, “Far-field linear optical superresolution via heterodyne detection in a higher-order local oscillator mode,” Optica 3, 1148–1152 (2016).
  * [18] Z. S. Tang, K. Durak, and A. Ling, “Fault-tolerant and finite-error localization for point emitters within the diffraction limit,” Optics express 24, 22004–22012 (2016).
  * [19] M. Paúr, B. Stoklasa, Z. Hradil, L. L. Sánchez-Soto, and J. Rehacek, “Achieving the ultimate optical resolution,” Optica 3, 1144–1147 (2016).
  * [20] W.-K. Tham, H. Ferretti, and A. M. Steinberg, “Beating rayleigh’s curse by imaging using phase information,” Physical review letters 118, 070801 (2017).
  * [21] M. Parniak, S. Borówka, K. Boroszko, W. Wasilewski, K. Banaszek, and R. Demkowicz-Dobrzański, “Beating the rayleigh limit using two-photon interference,” Physical review letters 121, 250503 (2018).
  * [22] U. Zanforlin, C. Lupo, P. W. Connolly, P. Kok, G. S. Buller, and Z. Huang, “Optical quantum super-resolution imaging and hypothesis testing,” Nature Communications 13, 5373 (2022).
  * [23] I. Ozer, M. R. Grace, and S. Guha, “Reconfigurable spatial-mode sorter for super-resolution imaging,” in _2022 Conference on Lasers and Electro-Optics (CLEO),_ (IEEE, 2022), pp. 1–2.
  * [24] A. Pushkina, G. Maltese, J. Costa-Filho, P. Patel, and A. Lvovsky, “Superresolution linear optical imaging in the far field,” Physical Review Letters 127, 253602 (2021).
  * [25] B. Zhang, J. Zerubia, and J.-C. Olivo-Marin, “Gaussian approximations of fluorescence microscope point-spread function models,” Applied optics 46, 1819–1829 (2007).
  * [26] M. Tsang, “Subdiffraction incoherent optical imaging via spatial-mode demultiplexing: Semiclassical treatment,” Physical Review A 97, 023830 (2018).
  * [27] W. H. Richardson, “Bayesian-based iterative method of image restoration,” JoSA 62, 55–59 (1972).
  * [28] Z. Wang, E. Simoncelli, and A. Bovik, “Multiscale structural similarity for image quality assessment,” in _The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003,_ vol. 2 (2003), pp. 1398–1402 Vol.2.
  * [29] H. Zhao, O. Gallo, I. Frosio, and J. Kautz, “Loss functions for image restoration with neural networks,” IEEE Transactions on computational imaging 3, 47–57 (2016).
  * [30] E. Bolduc, N. Bent, E. Santamato, E. Karimi, and R. W. Boyd, “Exact solution to simultaneous intensity and phase encryption with a single phase-only hologram,” Optics letters 38, 3546–3549 (2013).
  * [31] A. Pushkina, J. Costa-Filho, G. Maltese, and A. Lvovsky, “Comprehensive model and performance optimization of phase-only spatial light modulators,” Measurement Science and Technology 31, 125202 (2020).
  * [32] N. K. Fontaine, R. Ryf, H. Chen, D. T. Neilson, K. Kim, and J. Carpenter, “Laguerre-gaussian mode sorter,” Nature communications 10, 1865 (2019).
  * [33] C. Shorten and T. M. Khoshgoftaar, “A survey on image data augmentation for deep learning,” Journal of big data 6, 1–48 (2019).
  * [34] K. K. Bearne, Y. Zhou, B. Braverman, J. Yang, S. Wadood, A. N. Jordan, A. Vamivakas, Z. Shi, and R. W. Boyd, “Confocal super-resolution microscopy based on a spatial mode sorter,” Optics Express 29, 11784–11792 (2021).
  * [35] https://en.wikipedia.org/wiki/Weierstrass_transform#Transforms_of_some_important_functions.

Appendix: Derivation of Eqs. (1) and (2)

Consider an object with the electric field distribution in the object plane
given by $E(\vec{\rho})$. Imaging the object with the objective lens will
result in the electric field distribution in the image plane that is the
convolution of $E$ and the point spread function $\psi$:

$\displaystyle\tilde{E}(\vec{r})=\int
E(\vec{\rho})\psi(\vec{r}-\vec{\rho}){\rm d}\vec{\rho}.$ (S1)

Heterodyne detection of the image field with the local oscillator prepared in
a Hermite-Gaussian mode $\phi_{mn}(\vec{r})$ will produce the photocurrent

$\displaystyle J_{mn}=\int\phi_{mn}(\vec{r})\tilde{E}(\vec{r})=\iint
E(\vec{\rho})\psi(\vec{r}-\vec{\rho})\phi_{mn}(\vec{r}){\rm d}\vec{\rho}.$
(S2)

For incoherent objects, we measure the time-averaged power (absolute square)
of the heterodyne photocurrent, which is given by

$\displaystyle P_{mn}=\langle|J_{mn}|^{2}\rangle$
$\displaystyle=\iiiint\langle
E(\vec{\rho})E^{*}(\vec{\rho}^{\prime})\rangle\psi(\vec{\rho}-\vec{r})\psi^{*}(\vec{\rho}^{\prime}-\vec{r}^{\prime})\phi_{mn}^{*}(\vec{r})\phi_{mn}(\vec{r}^{\prime})\
{\rm d}\vec{\rho}{\rm d}\vec{\rho}^{\prime}{\rm d}\vec{r}{\rm
d}\vec{r}^{\prime}.$ (S3)

In the above result, we can recognise $\langle
E(\vec{\rho})E^{*}(\vec{\rho}^{\prime})\rangle$ as the mutual intensity
function. For a fully incoherent object with the intensity distribution
$I(\vec{\rho})$, the mutual intensity is

$\langle
E(\vec{\rho})E^{*}(\vec{\rho}^{\prime})\rangle=I(\vec{\rho})\delta(\vec{\rho}-\vec{\rho}^{\prime}).$

Hence

$\displaystyle P_{mn}=\iiint
I(\vec{\rho})\psi(\vec{\rho}-\vec{r})\psi^{*}(\vec{\rho}-\vec{r}^{\prime})\phi_{mn}^{*}(\vec{r})\phi_{mn}(\vec{r}^{\prime})\
{\rm d}\vec{\rho}{\rm d}\vec{r}{\rm d}\vec{r}^{\prime}.$ (S4)

The above integral can be written as

$\displaystyle P_{mn}=\int
I(\vec{\rho})\left|[\phi_{mn}*\psi](\vec{\rho})\right|^{2},$ (S5)

where $\phi_{mn}*\psi$ is the convolution of the Hermite-Gaussian mode

$\phi_{mn}(\vec{r})=\frac{H_{k}(\frac{x}{\sqrt{2}\sigma})H_{k}(\frac{y}{\sqrt{2}\sigma})}{\sqrt{2^{m+n}m!n!}}\psi(\vec{r})$

and the Gaussian PSF

$\displaystyle\psi(\vec{\rho})=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{x^{2}+y^{2}}{4\sigma^{2}}},$
(S6)

which is assumed to match $\phi_{00}$. This convolution is known to equal [35]

$\displaystyle\left[\phi_{mn}*\psi\right](x,y)=\frac{1}{\sqrt{m!n!}}\left(\frac{x}{2\sigma}\right)^{m}\left(\frac{y}{2\sigma}\right)^{n}e^{-\frac{x^{2}+y^{2}}{8\sigma^{2}}},$
(S7)

so we finally get

$\displaystyle P_{mn}$ $\displaystyle=\frac{1}{m!n!}\int
I(x,y)\left(\frac{x}{2\sigma}\right)^{2m}\left(\frac{y}{2\sigma}\right)^{2n}e^{-\frac{x^{2}+y^{2}}{4\sigma^{2}}}\
dxdy.$ (S8)

From these measurements, we can calculate the following sums with even indices
$p$ and $q$,

$\displaystyle\beta_{pq}=\sum_{\text{even
}k,l=0}^{p,q}\left(\frac{k}{2}\right)!\left(\frac{l}{2}\right)!\alpha_{pk}\alpha_{ql}P_{\frac{k}{2}\frac{l}{2}},$
(S9)

with the $\alpha$’s being the coefficients of Hermite polynomials
$H_{i}(x)=\sum_{j}^{i}\alpha_{ij}x^{j}$. These sums are equal to scalar
products between the object intensity distribution and the corresponding
Hermite-Gaussian polynomials:

$\displaystyle\beta_{pq}$ $\displaystyle=\sum_{\text{even
}k,l=0}^{p,q}\alpha_{pk}\alpha_{ql}\int
I(x,y)\left(\frac{x}{2\sigma}\right)^{k}\left(\frac{y}{2\sigma}\right)^{l}e^{-\frac{x^{2}+y^{2}}{4\sigma^{2}}}\
dxdy$ (S10) $\displaystyle=\sum_{\text{even }k,l=0}^{p,q}\int
I(x,y)\phi_{kl}(x,y)\ dxdy.$ (S11)

Because 2-dimensional Hermite-Gaussian polynomials $\phi_{pq}$ are orthonormal
with weight $2^{p+q+2}p!q!\pi\sigma^{2}$, we have that $\beta$’s are the
expansion coefficients of $I(x,y)$ into the Hermite-Gaussian basis:

$\displaystyle I_{\rm
even}(x,y)=\sum_{p,q=0}^{N}\frac{\beta_{pq}\phi_{pq}(x,y)}{2^{p+q+2}p!q!\pi\sigma^{2}},$
(S12)

Because the measurements only recover the $\beta$’s with even indices, we can
reconstruct the even component of the object intensity distribution $I(x,y)$.
In other words, if we calculate Eq. (1) from the main text including only even
$m$ and $n$, we will reconstruct the distribution
$[I(x,y)+I(-x,y)+I(x,-y)+I(-x,-y)]/4$.

Alternatively, we can recover the odd components by measuring in the
overcomplete measurement basis, which includes HG mode superpositions

$\phi_{mn}^{0-6}:=\\{\phi_{m,n},\phi^{\pm}_{m,m+1,n,n},\phi^{\pm}_{m,m,n,n+1},\phi^{\pm}_{m,m+1,n,n+1}\\},$

where

$\displaystyle\phi^{\pm}_{m,m^{\prime},n,n^{\prime}}:=\frac{1}{\sqrt{2}}[\phi_{m^{\prime},n^{\prime}}\pm\phi_{mn}].$
(S13)

By measuring the heterodyne photocurrent powers $P_{mn}^{(0)-(6)}$
corresponding to these mode superpositions, we can compute the integrals

$\displaystyle Q_{mn}$ $\displaystyle=\int
I(x,y)\left(\frac{x}{2\sigma}\right)^{m}\left(\frac{y}{2\sigma}\right)^{n}e^{-\frac{x^{2}+y^{2}}{4\sigma^{2}}}\
dxdy.$ (S14)

according to

$\displaystyle Q_{mn}$
$\displaystyle=\begin{cases}\left(\frac{m}{2}\right)!\left(\frac{n}{2}\right)!P^{(0)}_{\frac{m}{2}\frac{n}{2}}&m,n\text{
even}\\\
\left(\frac{m-1}{2}\right)!\left(\frac{n}{2}\right)!\sqrt{\frac{m+1}{8}}\left(P^{(1)}_{\frac{m-1}{2}\frac{n}{2}}-P^{(2)}_{\frac{m-1}{2}\frac{n}{2}}\right)&\text{$m$
odd, $n$ even}\\\
\left(\frac{m}{2}\right)!\left(\frac{n-1}{2}\right)!\sqrt{\frac{n+1}{8}}\left(P^{(3)}_{\frac{m}{2}\frac{n-1}{2}}-P^{(4)}_{\frac{m}{2}\frac{n-1}{2}}\right)&\text{$m$
even, $n$ odd}\\\
\left(\frac{m-1}{2}\right)!\left(\frac{n-1}{2}\right)!\sqrt{\frac{m+1}{8}}\sqrt{\frac{n+1}{8}}\left(P^{(5)}_{\frac{m-1}{2}\frac{n-1}{2}}-P^{(6)}_{\frac{m-1}{2}\frac{n-1}{2}}\right)&m,n\text{
odd.}\end{cases}$ (S15)

Knowing $Q_{mn}$, we can evaluate

$\displaystyle\beta_{pq}=\sum_{m,n=0}^{pq}\alpha_{pm}\alpha_{qn}Q_{mn}$ (S16)

for all $p$ and $q$, and hence reconstruct $I(x,y)$.