# Nonlinear interferometry-based metrology of magneto-optical properties at
infrared wavelengths

Tanmoy Chakraborty1 Thomas Produit1 Harish N S Krishnamoorthy2,3 Present
address: Tata Institute of Fundamental Research Hyderabad, Gopanpally,
Hyderabad 500046, Telangana, India  Cesare Soci2,3 Anna V. Paterova1 1A*STAR
Quantum Innovation Centre (Q.InC), Institute of Materials Research and
Engineering (IMRE), Agency for Science, Technology and Research (A∗STAR), 2
Fusionopolis Way, Innovis #08-03, Singapore 138634, Republic of Singapore
2Centre for Disruptive Photonic Technologies, The Photonics Institute, Nanyang
Technological University, 21 Nanyang Link, Singapore 637371, Republic of
Singapore  3Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link,
Singapore 637371, Republic of Singapore

###### Abstract

Magneto-optical properties of materials are utilized in numerous applications
both in scientific research and industries. The novel properties of these
materials can be further investigated by performing metrology in the infrared
wavelength range, thereby enriching their potential applications. However,
current infrared metrology techniques can be challenging and resource-
intensive due to the unavailability of suitable components. To address these
challenges, we propose and demonstrate a set of measurements based on
nonlinear interferometry, which allows us investigating magneto-optical
properties of materials at infrared wavelength range by performing optical
detection at the visible range. For a proof-of-principle study, we measure the
Verdet constant of a bismuth-iron-garnet, over a spectral bandwidth of 600 nm
in the near-IR range.

##  I Introduction

Quantum metrology exploits the principles of quantum mechanics and allows
measuring physical quantities with higher accuracy, better sensitivity and
precision compared to conventional classical methods [1, 2, 3, 4, 5]. Research
in quantum metrology has emerged substantially in recent years encompassing
measurements of a broad variety of physical parameters. This includes
extremely sensitive detection of fundamental constants of atoms and molecules
[6, 7], probing magnetic and electric fields using ensembles of cold atoms [8,
9, 10, 11], quantities associated with relativistic quantum field [12],
gravitational waves [13, 14, 15] and dark matter [16, 17]. The present work is
related to using quantum optics for metrology applications at infrared (IR)
wavelength range, where a number of interesting physical properties of solid
state and biological systems can be captured [18, 19]. Several applications
involve Faraday rotation spectroscopy in the infrared part of the spectrum as
well [20, 21].

Conventional IR metrology techniques require light sources and detectors
operating at IR wavelengths [22, 23]. Hence, the performance of these
techniques is limited by the significant inefficiency of existing IR detectors
and the need for their operation at cryogenic temperatures, which reduces
their accessibility and ease of use. The efficiency of the mercury cadmium
telluride detectors is limited due to their high dark counts, low sensitivity
even at cryogenic temperatures, and susceptibility to surrounding temperature
fluctuations [24, 25]. Alternatively, IR spectroscopy can be performed by up-
conversion of the IR signal, where an additional high power pump laser is
employed to convert the wavelength of the IR photons to visible (or near-IR),
where efficient detectors are easily available [26]. However, the detection
performance in such systems is strongly constrained due to low conversion
efficiency of nonlinear crystals [27].

Another method relies on the nonlinear interferometry of correlated photon
pairs generated via spontaneous parametric down-conversion (SPDC) [28] [29],
[30], where one of the photons is generated at visible wavelength (signal),
and the other at infrared (idler). Due to correlations between signal and
idler SPDC photons, interference for signal photons is sensitive to the
amplitude, phase and polarization of the electric field of idler photons.
Thus, light-matter interactions at IR wavelength can be studied by probing the
sample with idler photons and analyzing the interference for the signal
photons at visible wavelengths. The main advantage of this method is that it
does not require a detector for IR, the detection can be performed at visible
wavelength, where efficient Si-based detectors are available. Moreover, a wide
spectral tunability of the photons can be achieved by changing the phase
matching within the same nonlinear crystal. The notion of ”induced coherence
without induced emission” was introduced by Wang _el. al._ [30], and later
applied for various applications, including spectroscopy [28, 31, 32, 33, 34,
35, 36, 37], imaging [38, 39, 40, 41, 42, 43], as well as optical coherence
tomography [44, 45, 46, 47]. Polarization degree-of-freedom of the photons in
such scheme is exploited for metrology applications [48, 49], tomography [50]
and quantum erasers [51, 52].

Through the present work we demonstrate a polarimetry method for measuring the
magneto-optical properties of materials at IR wavelengths. Our technique is
based on nonlinear interferometry of non-degenerate SPDC photon pairs that are
generated by optically pumping a nonlinear crystal, which is embedded into the
Michelson interferometer. The idler photon generated at IR wavelength serves
as the probing photon, while the signal photon generated at a visible/near-IR
wavelength is detected. In such a system, any change in the polarization of
the idler photons introduced by the magneto-optical properties of the sample
at IR wavelength is inferred by studying the interference of the signal
photons. For a proof-of-principle demonstration of our method, we measure the
Vedet constant and saturation Faraday rotation of a Bi3Fe5O12 (BIG) crystal in
the wavelength range of 1540 - 2141 nm by performing detection at 820 - 708
nm, respectively.

##  II Theoretical models

Interference pattern of signal photons $I_{s}$ in a nonlinear interferometer
with optically active media at the idler arm of the interferometer is
expressed as follows [49] :

$I_{\mathrm{s}}\propto
1+{|\tau_{i}|^{2}}|t||\mu|\cos(\phi_{p}-\phi_{i}-\phi_{s}+\Phi),$ (1)

where $\absolutevalue{\tau_{i}}$ is amplitude transmission coefficient for the
idler photons through the optical components ($\absolutevalue{\tau_{i}}^{2}$
accounts for double pass of infrared photons though the specimen), $|t|$
determines polarization transfer, $\absolutevalue{\mu}$ is the normalized
correlation function for the SPDC photons, $\phi_{p}$, $\phi_{s}$, $\phi_{i}$
are acquired phases by pump, signal and idler photons, respectively [29],
$\Phi$ is the initial phase in the interferometer.

The visibility $\mathcal{V}$ of the interference pattern is defined from the
maxima and minima of interference fringes and is proportional to the losses in
the interferometer $|\tau_{i}|^{2}$, polarization transfer $|t|$ and
correlation function $\absolutevalue{\mu}$:

$\mathcal{V}=\frac{{I_{\mathrm{s}}^{\mathrm{max}}-I_{\mathrm{s}}^{\mathrm{min}}}}{{I_{\mathrm{s}}^{\mathrm{max}}+I_{\mathrm{s}}^{\mathrm{min}}}}\propto{|\tau_{i}|^{2}}|t||\mu|.$
(2)

The parameter $|t|$ is calculated employing Jones matrix formalism [49]. Here
we use combinations of half- (HWP), quarter- (QWP) waveplates and a polarizer
to observe the Faraday rotation introduced by a sample and construct two
schemes for the polarimetry measurements. For case I, the idler photons
propagate through the HWP, the sample, the polarizer, then reflect from the
mirror and travel back along the same path [see Fig. 1(a)]. For the case II,
the idler photons travel through the QWP and the sample, before reflecting
back from the mirror [see Fig. 1(b)]. Thus, the parameter $|t|$ is defined by
transformation matrices
$\mathcal{M}_{\mathrm{HWP}}({\theta}_{{}_{\mathrm{HWP}}})$,
$\mathcal{M}_{\mathrm{QWP}}({\theta}_{{}_{\mathrm{QWP}}})$,
$\mathcal{M}_{\mathrm{P}}({\theta}_{{}_{\mathrm{P}}})$,
$\mathcal{M}_{\mathrm{F}}({\theta}_{{}_{\mathrm{F}}})$ and
$\mathcal{M}_{\mathrm{M}}$ for HWP, QWP, polarizer, the sample with Faraday
rotation and mirror, respectively (see Appendix A). Here
$\theta_{{}_{\mathrm{HWP}}}$, $\theta_{{}_{\mathrm{QWP}}}$,
$\theta_{{}_{\mathrm{F}}}$, and ${\theta}_{{}_{\mathrm{P}}}$ are the
polarization rotation angles of the idler photons introduced by HWP, QWP, the
sample and the polarizer, respectively. Thus, in the scheme for case I the
combined Jones matrix after double pass through the optical components is
given by:

$J_{{}_{\mathrm{I}}}=\mathcal{M}_{\mathrm{HWP}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{P}}\,\mathcal{M}_{\mathrm{M}}\,\mathcal{M}_{\mathrm{P}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{HWP}}.$
(3)

The initial state $|\sigma\rangle$ of vertically polarized idler photons can
be written as $|\sigma\rangle=\begin{pmatrix}0&1\end{pmatrix}^{\intercal}$.
Hence, the resulting polarization state of the idler photons is described as:

$|\sigma\rangle_{{}_{\mathrm{I}}}=J_{{}_{\mathrm{I}}}|\sigma\rangle$ (4)

Since interference visibility for the signal photons is proportional to $|t|$
from Eq. 1, the visibility for the interference of signal photons in this
scheme is defined as a function of both $\theta_{{}_{\mathrm{HWP}}}$ and
$\theta_{\mathrm{F}}$, see Appendix A for more details:

$\mathcal{V}^{\mathrm{I}}=|\cos(2\theta_{{}_{\mathrm{HWP}}}+\theta_{\mathrm{F}})\cos(2\theta_{{}_{\mathrm{HWP}}}-\theta_{\mathrm{F}})|.$
(5)

Figure 1: Graphical representation of Jones matrix formalism and simulated
2-dimensional mapping showing variation of normalized visibility as a function
of Faraday rotation angle introduced by the sample and the rotation angle of
the (a) half-waveplate (with polarizer position fixed) and (b) quarter-
waveplate.

Next, the combined Jones matrix for case II can be written as:

$J_{{}_{\mathrm{II}}}=\mathcal{M_{\mathrm{QWP}}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{M}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{QWP}}.$
(6)

Then, in this case the normalized visibility as a function of
$\theta_{{}_{\mathrm{QWP}}}$ and $\theta_{\mathcal{\mathrm{F}}}$ is given by:

$\mathcal{V}^{\mathrm{II}}=|\cos(\theta_{{}_{\mathrm{QWP}}})\cos(\theta_{\mathrm{F}})|.$
(7)

Figure 1(a) [(b)] shows the mapping of normalized interference visibility as a
function of the Faraday rotation angle by the sample and the angle of rotation
of the HWP [QWP], calculated using Eq.(5) [Eq.(7)]. In both visibility maps,
the higher values correspond to better indistinguishability of the idler
photons, which is discussed more in details in the following sections. In the
experiments the visibilities can be measured depending on the orientation of
the waveplates for a given applied magnetic field (or the introduced Faraday
rotation). These measurements correspond to the horizontal cross-sections of
the graphs in Fig. 1.

##  III Experimental setup and measurement protocols

A simplified schematic of our experimental setup is shown in Fig. 2. Signal-
idler SPDC photon pairs are generated by focusing the continuous wave green
laser (532 nm wavelength, 40 mW power, Laser Quantum) pump beam into a
periodically poled lithium niobate (PPLN) nonlinear crystal (Covesion). The
PPLN crystal consists of multiple periodic poling channels, each designed to
meet different quasi-phase matching (QPM) conditions. QPM via different poling
and temperature tuning allows generating idler photons over 1500 nm - 2300 nm
wavelength range, while conjugate signal photons are generated at 824 nm to
692 nm. After the PPLN crystal the co-propagating signal and idler beams are
separated using a dichroic mirror (DM1) and collimated using two lenses, each
with $f$=125 mm focal length. The signal and idler photons propagate along the
two arms of the Michelson interferometer and reflect back from two metallic
mirrors. The mirror Mi at the idler arm is mounted on a motorized translation
stage, which allows to introduce the additional phase for the idler photons.
The pump beam following the signal path is reflected back by the mirror Ms
into the PPLN crystal, which generates the signal and idler SPDC photons for
the second time. After separating the signal and pump photons by a second
dichroic mirror (DM2), the signal photons are focused onto a sCMOS camera
(Thorlabs, Inc.). At the input of the sCMOS camera we use a notch at 532 nm
and a band-pass filters to block any residual pump photons.

Figure 2: A simplified schematic of the experimental setup based on Michelson
interferometer. PPLN: periodically poled lithium niobate; DM1: dichroic mirror
that transmits idler photons and reflects signal and pump photons; DM2:
dichroic mirror that reflects pump and transmits signal; Ms and Mi: signal and
idler mirrors, respectively. Interference of the signal photons are captured
using the sCMOS camera.

The BIG sample under study is placed at the path of idler photons (see the
crystal structure of the sample and its magnetization properties given in
Appendix B). An externally applied magnetic field $B$ to the sample along the
optical axis results in Faraday rotation of the polarization of idler photons.
Below the saturation field, the Faraday rotation angle ${\Theta}_{\mathrm{F}}$
increases linearly with the strength of the applied magnetic field as follows:

${\Theta}_{\mathrm{F}}=V(\lambda)Bl,$ (8)

where $l$ is length of the crystal along the light propagation direction (in
our case $l=$0.38 mm). Thus, at the linear regime of the Faraday rotation it
is possible to measure the Verdet constant $V(\lambda)$.

To detect the Faraday rotation introduced by the sample we use two different
configurations of the idler arm of the interferometer. In case I, after the
DM1 idler photons pass through the HWP, BIG sample and polarizer. The HWP
introduces a known polarization transformation to idler photons, which is
followed by Faraday rotation from the BIG crystal. Subsequently, the polarizer
retrieves the polarization back to its initial state, as it was before the
HWP. After reflecting back from Mi, the idler photons travel back along the
same path through the polarizer, BIG and HWP. The polarization state is not
changed by the polarizer, though it is altered by the Faraday rotation in the
BIG crystal and the retardation introduced by the HWP. In case II, a known
amount of retardation is introduced to the idler photons by QWP, which follows
the Faraday rotation by the BIG crystal. These two operations are then applied
in reverse order after the photons are reflected back by mirror Mi.

We analyze the resulting polarization state of the idler photons by measuring
the visibilities of the interference fringes for the signal photons. These
fringes are observed via introducing a relative phase by translating the
mirror Mi along the optical axis of the idler arm. Thus, by capturing the
variation of visibility pattern as a function of Faraday rotation angles at
different applied magnetic fields $B$, we study the magneto-optical properties
of BIG at IR wavelengths, according to Eqs. 5, 7 and 8.

##  IV Experimental results

First, we pump the PPLN crystal along a polling period that allows generating
the signal and idler photons at 813 nm and 1540 nm, respectively (see the
signal spectrum in Appendix C). The full width at half maximum of the spectrum
is 1.4$\pm$0.02 nm, which gives us $l_{{}_{\mathrm{coh}}}=0.48\pm 0.01$ mm
coherence length. We measure the intensity of signal photons across $\approx$
0.8 mm path length difference without and with the sample at the path of idler
photons (see Figs. 8 and 9 in Appendix C) and confirm the calculated coherence
length experimentally. The inclusion of the BIG crystal into the path of idler
arm changes the optical path length to $\approx$ 0.53 mm, which depends on the
refractive index of BIG and its thickness $l$. In this case the optical path
length difference can be compensated by shifting the position of the mirror
Mi.

First, we perform interference measurements for case I and II without
introducing the BIG crystal into the interferometer (see Figs. 10 and 11 in
Appendix D for the case I and II, respectively). The visibility of the
interference fringes depends on the orientation of the waveplates
(${\theta}_{\mathrm{HWP}}$ for case I and ${\theta}_{\mathrm{QWP}}$ for case
II) as calculated by our theoretical model in Eq. 5 for case I and Eq. 5 for
case II. Next, we introduce the BIG sample, balance the interferometer arms
and perform the polarimetry measurements.

###  IV.1 Results for case I: with HWP and polarizer

We measure the visibility of the interference fringes for the signal photons
with the orientation of the HWP at zero magnetic field $B=0$ mT, which is
shown in Fig. 3(a). The normalized visibility $\mathcal{V}$ is maximum when
${\theta}_{\mathrm{HWP}}=0^{\circ}$. It reaches minimum value when
${\theta}_{\mathrm{HWP}}=45^{\circ}$, and restores its maximum again when
${\theta}_{\mathrm{HWP}}=90^{\circ}$, which is in accordance with Eq. 5 for
${\theta}_{\mathrm{F}}=0^{\circ}$.

Figure 3: (a) Visibility of interference fringes of signal photons depending
on the rotation angle of the HWP at different magnetic field strengths applied
to the BIG crystal. (b) Spatial intensity distribution of the signal photons
revealing constructive interference patterns captured using the sCMOS camera.
A few of such patterns are shown for different magnetic field strengths and
HWP rotation angles.

The visibility as a function of ${\theta}_{\mathrm{HWP}}$ for different
strengths of the magnetic field is shown in Fig. 3(a) as well. Interestingly,
initially one minimum in visibility curve at $B=0$ mT splits into two for
$B\neq 0$ mT. The separation of these minima increases with the strength of
the applied magnetic field and reaches maximum for the saturation field. These
results are consistent with the theoretical model given by Eq. 5. We fit the
experimental data for $B=75$ mT and 141 mT with keeping
${\theta}_{\mathrm{F}}$ as the fitting parameter, which yields
${\theta}_{F}=34.08^{\circ}\pm 0.22^{\circ}$ and $44.34^{\circ}\pm
0.46^{\circ}$, respectively.

Change in visibility of the interference pattern is associated with
indistinguishability of SPDC photons generated in two passes of the pump beam
through the nonlinear crystal. The waveplate and Faraday rotator change the
polarization state of idler photons, which reduces the degree of
indistinguishability and hence, the interference visibility of signal photons
[49]. The stronger magnetic field causes larger Faraday rotation, reducing the
degree of indistinguishability further. As an example, Fig. 3(b) shows a few
measured spatial intensity distributions of signal photons for different
values of ${\theta}_{\mathrm{HWP}}$ and $B$.

###  IV.2 Results for case II: with QWP

In this scheme, we replace the HWP shown in Fig. 2 with a QWP, remove the
polarizer and perform manipulation of the input polarization of idler photons
by rotating the QWP. The dependence of the visibility of the interference
fringes on the orientation of the QWP for different values of $B$ is shown in
Fig. 4. For $B=0$ mT and $\theta_{\mathrm{QWP}}=0^{\circ}$, the photons
generated during the forward and backward pass of the pump are
indistinguishable, being in the same polarization state. However, the degree
of indistinguishability is reduced with increasing $\theta_{QWP}$, and the
visibility of the interference fringes reaches minimum when
$\theta_{QWP}=45^{\circ}$. Additional polarization rotation introduced by the
BIG varies with the strength of the magnetic field $B$, and the maximum
visibility of the interference fringes are reduced with increasing the
magnetic field strength. Visibility curves depending on the orientation of the
QWP for each of the values of $B$ are fitted to Eq. 7, wheres, the Faraday
rotation angles caused by the BIG crystal are determined (see the summarized
Faraday rotation angles depending on the applied magnetic field in table 1 in
Appendix E). For the applied magnetic fields $B=75$ mT and 141 mT the
resulting Faraday rotation angles are ${\theta}_{F}=31.24^{\circ}\pm
0.15^{\circ}$ and $44.9^{\circ}\pm 0.45^{\circ}$, respectively, which agree
with the measured values in case I.

Figure 4: Visibility of the interference fringes for the signal photons
measured as a function of the rotation angle of QWP for different magnetic
field strengths, causing different amount of Faraday rotations.

###  IV.3 Wavelength dependent Faraday rotation

Next, we generate the idler photons at IR wavelength range from
$\lambda_{\mathrm{idler}}$ = 1516 nm to 2141 nm. The coarse tuning of
$\lambda_{\mathrm{idler}}$ is done by changing the poling period and fine-
tuning is done by changing the applied temperature to the PPLN crystal while
being in the same poling channel. Following the energy conservation law, the
detected signal photons are generated in the range of 708 nm to 820 nm (see
table 1 in Appendix F for the details). To capture the dependence of the
Verdet constant $V(\lambda)$ on wavelength we choose the scheme in case I of
the interferometer, as both the maxima and the minima of the visibility
undergoes changes with the applied magnetic field. For each wavelength we
capture Faraday rotations for different strengths of the applied magnetic
field $B$. Following same protocol as for the data captured at $\lambda$=1540
nm, we measure $V(\lambda)$ at IR wavelengths over 600 nm bandwidth range [see
Fig. 5(a)].

Figure 5: (a) The Verdet constant depending on wavelength. Points correspond
to data measured employing the interferometry scheme. Inset shows Faraday
rotations caused by the BIG crystal with applied magnetic field and its linear
fit for $\lambda_{{}_{\mathrm{\mathrm{idler}}}}$=1540 nm. (b) Maximum Faraday
rotation (at saturation magnetic field) introduced by the BIG sample depending
on wavelength: the hexagons represent the experimental data measured in
interferometry scheme, and the solid line show the values provided by
manufacturer.

The estimated Faraday rotation angles for case I are plotted depending on $B$
in the inset of Fig. 5. The data that are below the saturation field $B_{sat}$
is fitted as a linear curve. The slope of the fit allows estimating the Verdet
constant $V(\lambda)$ for BIG crystal at $\lambda$=1540 nm. Note that beyond
the saturation field at $B=150$ mT the Faraday rotation angle is
${\theta}_{F}=45.0^{\circ}\pm 0.5^{\circ}$, which corresponds to the supposed
specifications. The plot shows that $V(\lambda)$ decreases with the wavelength
and becomes significantly small at 2141 nm. From the acquired data we also
estimate the Faraday rotation angles at saturation magnetic field, which is
shown the in Fig. 5(b). The data is well consistent with their classically
measured counterparts.

##  V Discussion and conclusions

To summarize, we have demonstrated a metrology scheme based on a nonlinear
interferometry to investigate the magneto-optical properties of a specimen at
IR wavelengths. The method is implemented by using a Michelson interferometer,
where the paths of the signal and idler SPDC photons construct the two arms of
the interferometer. Through our experiments, a polarization rotation of the
idler photons introduced by the sample at IR wavelength is inferred by
studying the interference of the signal photons. This allows performing
measurement of magneto-optical properties of the sample at IR wavelengths by
performing detection at visible or near-IR wavelengths.

For a proof-of-concept demonstration of our method, we performed metrology of
magneto-optical properties of Bi3Fe5O12 crystal, which shows detectable
Faraday rotation at IR wavelength range. We formulated theoretical analysis
for calculating the visibility of the interference for signal photons as a
function of the polarization rotation introduced by the sample at the idler
photons wavelength. We captured the Verdet constant and Faraday rotation
angles over a bandwidth from 1516 to 2141 nm, while the detection is performed
at the range from 820 nm to 708 nm. The results obtained through our method
show consistency with the classically measured data, which proves the
applicability of our method. By choosing different QPM conditions in PPLN
crystal used in the experiments, the metrology experiments can be performed at
even longer wavelengths, up to 5 $\mu$m, which is given by the transparency
range of the PPLN crystal. The method can be extended even further towards
measuring materials which show detectable Faraday rotation at longer
wavelengths (5-10 $\mu$m) using suitable nonlinear crystals like silver
gallium sulfide [53] and gallium phosphide [54].

The proposed IR metrology protocol and its experimental demonstration in the
present paper open up possibilities for unleashing magneto-optical properties
of a large class of materials, that is not yet well explored due to technical
challenges. Such investigations can reveal novel phenomena of these materials
in the IR wavelength range, which can significantly boost their applications
in scientific research and industries. IR polarimetry can have applications in
defence and facial recognition at night [55] and in characterization of
painting materials [56]. These rely on polarimetry of reflection and thermal
emission from surfaces, which is rather complex to analyse. Such measurements
can be benefited from our IR polarimetry technique. Our method also provides a
new path towards investigating light-matter interaction phenomena in novel
material systems such as topological insulators, which have exotic surface
states that are sensitive to the helicity of light [57, 58].

##  VI Acknowledgements

We thank Jian Rui Soh and Vytautas Valuckas for their help with characterizing
the sample. AP, TP and TC acknowledge support to this research A*STAR through
the grant no. C230917004. and for support by NRF Singapore through the grant
NRF2021-QEP2-03-P08.

## Appendix A: Calculation of the interference visibilities for schemes of the
interferometer in _case I_ and _case II_

To calculate the parameter $t$ expressed in Eq. 1, we first determine the
transfer matrix for the optical components described in section III for cases
I and II. Change in the polarization of the idler photons by each of the
optical components is determined by employing Jones matrix formalism.
According to the schemes of the interferometer we express Jones matrices for
half- and quarter- wave plates, Faraday rotator, polarizer, and a mirror as:
$\mathcal{M}({\theta}_{{}_{\mathrm{HWP}}})$,
$\mathcal{M}({\theta}_{{}_{\mathrm{QWP}}})$,
$\mathcal{M}({\theta}_{{}_{\mathrm{F}}})$, $\mathcal{M}_{\mathrm{P}}$ and
$\mathcal{M}_{\mathrm{M}}$, respectively:

$\mathcal{M}_{\mathrm{HWP}}=\begin{pmatrix}\cos
2\theta_{{}_{\mathrm{HWP}}}&\sin 2\theta_{{}_{\mathrm{HWP}}}\\\ \sin
2\theta_{{}_{\mathrm{HWP}}}&-\cos 2\theta_{{}_{\mathrm{HWP}}}\end{pmatrix},$
(9)

$\mathcal{M}_{\mathrm{QWP}}=\frac{1}{\sqrt{2}}\begin{pmatrix}1+i\cos
2\theta_{{}_{\mathrm{QWP}}}&i\sin 2\theta_{{}_{\mathrm{QWP}}}\\\ i\sin
2\theta_{{}_{\mathrm{QWP}}}&1-i\cos 2\theta_{{}_{\mathrm{QWP}}}\end{pmatrix},$
(10)

$\mathcal{M}_{\mathrm{F}}=\begin{pmatrix}\cos\theta_{{}_{\mathrm{F}}}&\sin\theta_{{}_{\mathrm{F}}}\\\
-\sin\theta_{{}_{\mathrm{F}}}&\cos\theta_{{}_{\mathrm{F}}}\end{pmatrix},$ (11)

$\mathcal{M}_{\mathrm{P}}=\begin{pmatrix}0&0\\\
0&1\end{pmatrix},\mathcal{M}_{\mathrm{M}}=\begin{pmatrix}1&0\\\
0&1\end{pmatrix},$ (12)

where angles $\theta_{{}_{\mathrm{HWP,QWP,F}}}$ denote angles between initial
polarization state of the idler photons and an orientation of the waveplates,
or Faraday rotation. In Eq. 12 the matrix for polarizer
$\mathcal{M}_{\mathrm{P}}$ is given for the vertical orientation, which is
aligned with the initial polarization for the idler photons. The matrix
$\mathcal{M}_{\mathrm{M}}$ is given by the unit matrix since the mirror does
not alter the polarization of the propagating beams.

In case I, the polarization is changed by the HWP, BIG sample and polarizer at
the forward path along the arm of the Michelson interferometer:

$\mathcal{M}_{{}_{1}}^{\mathrm{I}}=\mathcal{M}_{\mathrm{HWP}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{P}}=\begin{pmatrix}0&\sin{(2\theta_{{}_{\mathrm{F}}}+\theta_{{}_{\mathrm{HWP}}})}\\\
0&\cos{(2\theta_{{}_{\mathrm{F}}}+\theta_{{}_{\mathrm{HWP}}})}\end{pmatrix}.$
(13)

After the reflection from Mi, the combined matrix is given in reversed order
$\mathcal{M}_{{}_{2}}^{\mathrm{I}}=\mathcal{M}_{\mathrm{P}}\,\mathcal{M}_{\mathrm{F}}\,\mathcal{M}_{\mathrm{HWP}}$.
Hence the total transfer matrix is defined as:

$\displaystyle
J_{\mathrm{I}}=\mathcal{M}_{{}_{1}}^{\mathrm{I}}\,\mathcal{M}_{\mathrm{M}}\,\mathcal{M}_{{}_{2}}^{\mathrm{I}}=$
(14)
$\displaystyle\frac{1}{2}\begin{pmatrix}\cos{2\theta_{{}_{\mathrm{F}}}}-\cos{4\theta_{{}_{\mathrm{HWP}}}}&-\sin{2\theta_{{}_{\mathrm{F}}}}-\sin{4\theta_{{}_{\mathrm{HWP}}}}\\\
\sin{2\theta_{{}_{\mathrm{F}}}}-\sin{4\theta_{{}_{\mathrm{HWP}}}}&\cos{2\theta_{{}_{\mathrm{F}}}}+\cos{4\theta_{{}_{\mathrm{HWP}}}}\end{pmatrix}.$

The initial state $|\sigma\rangle$ of vertically polarized idler photons can
be written as $|\sigma\rangle=\begin{pmatrix}0&1\end{pmatrix}^{\intercal}$.
The final state of the idler photons after using the Jones matrix and
simplification of the equations is given by:

$|\sigma\rangle_{{}_{\mathrm{I}}}=J_{{}_{\mathrm{I}}}\begin{pmatrix}0\\\
1\end{pmatrix}=\frac{1}{2}\begin{pmatrix}-\sin{2\theta_{{}_{\mathrm{F}}}}-\sin{4\theta_{{}_{\mathrm{HWP}}}}\\\
\cos{2\theta_{{}_{\mathrm{F}}}}+\cos{4\theta_{{}_{\mathrm{HWP}}}}\end{pmatrix}.$
(15)

The resulting interference depends on the vertical polarization component.
Thus, the visibility of the interference fringes depends on the orientation of
the HWP and the Faraday rotation as follows:

$\displaystyle\mathcal{V}^{\mathrm{I}}=\frac{1}{2}|(\cos{2\theta_{{}_{\mathrm{F}}}}+\cos{4\theta_{{}_{\mathrm{HWP}}}})|=$
(16)
$\displaystyle|\cos{(2\theta_{{}_{\mathrm{HWP}}}+\theta_{{}_{\mathrm{F}}})}\cos{(2\theta_{{}_{\mathrm{HWP}}}-\theta_{{}_{\mathrm{F}}})}|.$

Similarly, in case II, the Jones matrices for combined QWP and BIG crystal
towards Mi and backwards from Mi are given respectively by:

$\displaystyle\mathcal{M}_{{}_{1}}^{\mathrm{II}}=\frac{\sqrt{2}}{2}\cdot$ (17)
$\displaystyle\begin{pmatrix}\cos{\theta_{{}_{\mathrm{F}}}}+i\cos{(2\theta_{{}_{\mathrm{QWP}}}+\theta_{{}_{\mathrm{F}}})}&\sin{\theta_{{}_{\mathrm{F}}}}+i\sin{(2\theta_{{}_{\mathrm{QWP}}}+\theta_{{}_{\mathrm{F}}})}\\\
-\sin{\theta_{{}_{\mathrm{F}}}}+i\sin{(2\theta_{{}_{\mathrm{QWP}}}+\theta_{{}_{\mathrm{F}}})}&\cos{\theta_{{}_{\mathrm{F}}}}-i\cos{(2\theta_{{}_{\mathrm{QWP}}}+\theta_{{}_{\mathrm{F}}})}\end{pmatrix}.$

Hence, the Jones matrix for double-pass configuration is
$J_{\mathrm{II}}=\mathcal{M}_{{}_{1}}^{\mathrm{II}}\,\mathcal{M}_{\mathrm{M}}\,\mathcal{M}_{{}_{2}}^{\mathrm{II}}$.
The final visibility depending on the orientation of the QWP and the Faraday
rotation by the BIG sample is given by:

$\mathcal{V}^{\mathrm{II}}=|\cos(\theta_{{}_{\mathrm{QWP}}})\cos(\theta_{\mathrm{F}})|.$
(18)

## Appendix B: Characterization of the Bi3Fe5O12 sample

We perform magnetization and X-ray diffraction (XRD) measurements for
characterizing of B3Fe5O12 crystal used in our measurements. Fig. 6 shows the
isothermal magnetization curve measured (using a superconducting quantum
interference-based magnetometer by Quantum Design) as a function of externally
applied magnetic field at 300 K. The data shows that magnetization saturates
when the magnetic field is approximately 100 mT, which is consistent with the
magnetic field dependent Faraday rotation results shown in the inset of Fig.
5(a). The upper inset of Fig. 6 shows the $\theta$-2$\theta$ XRD scanning
depicting only the (444) peak, that confirms the single crystallinity of the
BIG sample [59, 60]. The unit cell is cubic with $a$ = 12.493 $\AA$ (see the
lower inset of Fig. 6).

Figure 6: Isothermal magnetization data as a function of applied magnetic
field. Upper inset shows single-crystal $\theta$-2$\theta$ X-ray diffraction
scan data and lower inset shows the cubic lattice structure of BIG.

## Appendix C: Coherence length of the idler photons at 1540 nm

The spectrum of the signal SPDC photons generated from PPLN crystal with 7.4
$\mu$m poling periodicity and T=126∘C applied temperature is shown in Fig. 7.
The central wavelength of the spectrum is 813$\pm$0.2 nm, the bandwidth is
1.40$\pm$0.02 nm, which gives the 0.48$\pm$0.01 mm coherence length.

Figure 7: Optical spectrum of the signal photons at 813 nm. The data is fitted
to a Sinc2 curve. The linewidth is 1.40$\pm$0.02 nm, which is the QPM
bandwidth of the photons at near-IR range.

Intensity of the signal photons captured as a function of the translation of
M${}_{\emph{i}}$ (MTS25/M, Thorlabs Inc.) shows the interference of the signal
photons in Fig. 8. The interference is observable within the coherence length
of $\approx$ 0.44 of the SPDC photons, which is twice the width of the
envelope (due to double pass configuration) taken at the $1/e$ amplitude decay
of the interference fringes [61], see Fig. 8. In the balanced position of the
interferometer, the visibility of the interference is maximum. The inset at
the bottom of Fig. 8 shows an enlarged view of the interference fringes close
to the balanced position of the interferometer. Due to the double-pass along
the arms of the Michelson-type interferometer, the periodicity of the
interference fringes corresponds to half of the idler wavelength. Fitting the
data to a cosine curve gives the periodicity of the interference fringes,
which is 759 nm. We attribute the error in this value to the irreversibility
caused by the backlash error of the transnational stage and other mechanical
components.

Figure 8: Interference fringes depending on the additional phase introduced by
translating mirror $M_{i}$ of the interferometer when the sample is absent.
The inset at the bottom shows the zoomed in interference data around the
balanced position of the interferometer and its fit to cosine function.
Figure 9: Interference fringes depending on the translation of the mirror
$M_{i}$ when the sample is present. The interference fringes at the optical
path length difference of 0.53 mm are shown in the inset at the bottom. The
fitting is performed using cosine function.

The interference fringes after the insertion of the BIG crystal into the idler
arm of the interferometer are shown in Fig. 9. Due to change in the optical
path length, the balanced position of the interferometer is shifted by 0.53
mm. This change of the optical path length corresponds to the 0.38 mm
thickness of the BIG sample and depends on the wavelength of the idler photons
(as refractive index of a medium varies with wavelength of the light according
to the Sellmeier relation [62]). The inset at the bottom in Fig. 9 shows an
enlarged view of the interference fringes and their fit to cosine curve.

## Appendix D: Visibility of the interference without the sample

In case I, in the absence of BIG sample, the visibility of interference
fringes for signal photons depends on the orientation of the HWP and polarizer
(see the inset in Fig. 10). Since the polarizer is aligned with the initial
polarization of the idler photons, the visibility of the interference is given
by the orientation of the HWP only. For each rotation angle of the HWP we
introduce additional phase into the interferometer by translating mirror
M${}_{\emph{i}}$ over a few periods of $\lambda_{\mathrm{idler}}$, capture the
interference fringes and estimate the visibilities, see Fig. 10. The
constructive interference pattern of the signal photons captured by the sCMOS
camera at angles of HWP at 0∘, 45∘ and 90∘ are shown in Fig. 10(b).

Figure 10: (a) The visibility of the interferene at the absence of the sample
as a function of HWP orientation. (b) The examples of the constructive
interference pattern captured by sCMOS camera at few orientation angles of the
HWP.  Figure 11: (a) The visibility of the interference at the absence of the
sample as a function of QWP orientation. (b) The examples of the constructive
interference pattern captured by sCMOS camera at few orientation angles of the
QWP.

In case II, when we replace the HWP by a QWP and remove the polarizer (see the
inset), the balance condition for the interferometer changes. We balance the
interferometer by translating the mirror M${}_{\emph{i}}$. Next, we record the
change in visibility of the interference depending on the orientation of QWP
in the absence of the BIG sample (see Fig. 11). Likewise, the constructive
interference pattern of the signal photons at orientation angles of QWP at 0∘,
45∘ and 90∘ are shown Fig. 11(b). Both of the obtained visibility data is in a
good agreement with the theoretical analysis given by Eq. 5 (case I) and Eq.
10 (case II).

## Appendix E: Faraday rotation angles measured via nonlinear interferometry

Table 1: Applied magnetic field strengths and measured Faraday rotations. Magnetic field (case I) |  Faraday rotation(case I)
---|---
$75\pm 8$ mT | $34.08^{\circ}\pm 0.22^{\circ}$
$141\pm 15$ mT | $44.34^{\circ}\pm 0.46^{\circ}$
Magnetic field (case II) |  Faraday rotation(case II)
$36.5\pm 3.9$ mT | $16.46^{\circ}\pm 0.47^{\circ}$
$75\pm 8$ mT | $31.24^{\circ}\pm 0.15^{\circ}$
$106\pm 4.1$ mT | $42.31^{\circ}\pm 0.08^{\circ}$
$141\pm 15$ mT | $44.9^{\circ}\pm 0.45^{\circ}$

Here in table 1 we show the experimentally measured Faraday rotations
depending on the applied magnetic fields in case I and II.

## Appendix F: Wavelength dependent Faraday rotation

To perform metrology of the Verdet constant and the Faraday rotation at the
saturation magnetic field over a broad IR wavelength range, we generate the
idler photons in the range of 1540 nm to 2141 nm by pumping the PPLN through
different quasi phase-matched channels and setting the corresponding
temperatures. We measure the signal (visible/near-IR) photon spectra using a
spectrometer (HR4000, Ocean Optics) and determine the signal wavelengths.
Following the energy conservation, we determine the wavelengths of idler
photons. Table 2 summarizes signal$\lambda_{\mathrm{idler}}$ and idler
wavelengths $\lambda_{\mathrm{idler}}$, with corresponding values of the
polling periods and PPLN temperatures.

Table 2: Poling periods of the PPLN and temperatures for generating signal and
idler photons at different wavelengths. The wavelength of the pump laser is
532 nm in all these cases.
Polling period | Temperature |  Signal wavelength (detected) |  Idler wavelength (Probed)
---|---|---|---
7.4 $\mu$m | 378 K | 820 nm | 1516 nm
7.4 $\mu$m | 399 K | 813 nm | 1540 nm
7.4 $\mu$m | 423 K | 805 nm | 1568 nm
7.71 $\mu$m | 303 K | 793 nm | 1617 nm
8.36 $\mu$m | 352 K | 725 nm | 2000 nm
8.36 $\mu$m | 571 K | 708 nm | 2141 nm
Figure 12: Variation of visibilities of interference of the signal photons
with the rotation angle of the HWP for different idler wavelengths.

For each of these cases, we measure visibility of the interfering signal
photons as a function of the input polarization of the idler photons. Figure
12 shows the visibilities of the interference pattern depending on HWP
rotation angle (protocol for the case I) and their fit using Eq. 5. We
determine $\theta_{\mathrm{F}}$ from the fitting and estimate the Verdet
constant. For the estimation of the Verdet constant we choose the magnetic
fields below the saturation field. At the saturation magnetic field we
calculate the Faraday rotation angles as a function of wavelength.

## References

  * Göbel and Siegner [2015] E. O. Göbel and U. Siegner,  _Quantum Metrology: Foundation of Units and Measurements_ (John Wiley & Sons, 2015).
  * Barbieri [2022] M. Barbieri, PRX Quantum 3, 010202 (2022).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Nature photonics 5, 222 (2011).
  * Tóth and Apellaniz [2014] G. Tóth and I. Apellaniz, Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Mukamel  _et al._ [2020] S. Mukamel, M. Freyberger, W. Schleich, M. Bellini, A. Zavatta, G. Leuchs, C. Silberhorn, R. W. Boyd, L. L. Sánchez-Soto, A. Stefanov, _et al._ , Journal of physics B: Atomic, molecular and optical physics 53, 072002 (2020).
  * Rosenband  _et al._ [2008] T. Rosenband, D. Hume, P. Schmidt, C.-W. Chou, A. Brusch, L. Lorini, W. Oskay, R. E. Drullinger, T. M. Fortier, J. E. Stalnaker, _et al._ , Science 319, 1808 (2008).
  * Morel  _et al._ [2020] L. Morel, Z. Yao, P. Cladé, and S. Guellati-Khélifa, Nature 588, 61 (2020).
  * Vengalattore  _et al._ [2007] M. Vengalattore, J. Higbie, S. Leslie, J. Guzman, L. Sadler, and D. Stamper-Kurn, Physical review letters 98, 200801 (2007).
  * Wildermuth  _et al._ [2006] S. Wildermuth, S. Hofferberth, I. Lesanovsky, S. Groth, P. Krüger, J. Schmiedmayer, and I. Bar-Joseph, Applied physics letters 88 (2006).
  * Behbood  _et al._ [2013] N. Behbood, F. Martin Ciurana, G. Colangelo, M. Napolitano, M. W. Mitchell, and R. J. Sewell, Applied physics letters 102 (2013).
  * Pezze  _et al._ [2018] L. Pezze, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Reviews of Modern Physics 90, 035005 (2018).
  * Ahmadi  _et al._ [2014] M. Ahmadi, D. E. Bruschi, and I. Fuentes, Physical Review D 89, 065028 (2014).
  * Schnabel  _et al._ [2010] R. Schnabel, N. Mavalvala, D. E. McClelland, and P. K. Lam, Nature communications 1, 121 (2010).
  * Tse  _et al._ [2019] M. e. Tse, H. Yu, N. Kijbunchoo, A. Fernandez-Galiana, P. Dupej, L. Barsotti, C. Blair, D. Brown, S. Dwyer, A. Effler, _et al._ , Physical Review Letters 123, 231107 (2019).
  * Danilishin  _et al._ [2019] S. L. Danilishin, F. Y. Khalili, and H. Miao, Living Reviews in Relativity 22, 1 (2019).
  * Backes  _et al._ [2021] K. M. Backes, D. A. Palken, S. A. Kenany, B. M. Brubaker, S. Cahn, A. Droster, G. C. Hilton, S. Ghosh, H. Jackson, S. K. Lamoreaux, _et al._ , Nature 590, 238 (2021).
  * Dixit  _et al._ [2021] A. V. Dixit, S. Chakram, K. He, A. Agrawal, R. K. Naik, D. I. Schuster, and A. Chou, Physical review letters 126, 141302 (2021).
  * Colarusso  _et al._ [1998] P. Colarusso, L. H. Kidder, I. W. Levin, J. C. Fraser, J. F. Arens, and E. N. Lewis, Applied Spectroscopy 52, 106A (1998).
  * Bhargava [2012] R. Bhargava, Applied spectroscopy 66, 1091 (2012).
  * Lewicki  _et al._ [2009] R. Lewicki, J. H. Doty III, R. F. Curl, F. K. Tittel, and G. Wysocki, Proceedings of the national academy of sciences 106, 12587 (2009).
  * Liu  _et al._ [2023] C. Liu, Y. Hai, and A. Dang, Physical Review A 107, 053105 (2023).
  * Smith [2011] B. C. Smith,  _Fundamentals of Fourier transform infrared spectroscopy_ (CRC press, 2011).
  * Fard  _et al._ [2013] A. M. Fard, P. Vacas-Jacques, E. Hamidi, H. Wang, R. W. Carruth, J. A. Gardecki, and G. J. Tearney, Optics express 21, 30849 (2013).
  * Rogalski [2012] A. Rogalski, Progress in Quantum Electronics 36, 342 (2012).
  * Lei  _et al._ [2015] W. Lei, J. Antoszewski, and L. Faraone, Applied Physics Reviews 2 (2015).
  * Junaid  _et al._ [2019] S. Junaid, S. C. Kumar, M. Mathez, M. Hermes, N. Stone, N. Shepherd, M. Ebrahim-Zadeh, P. Tidemand-Lichtenberg, and C. Pedersen, Optica 6, 702 (2019).
  * Israelsen  _et al._ [2019] N. M. Israelsen, C. R. Petersen, A. Barh, D. Jain, M. Jensen, G. Hannesschläger, P. Tidemand-Lichtenberg, C. Pedersen, A. Podoleanu, and O. Bang, Light: Science & Applications 8, 11 (2019).
  * Kalashnikov  _et al._ [2016] D. A. Kalashnikov, A. V. Paterova, S. P. Kulik, and L. A. Krivitsky, Nature Photonics 10, 98 (2016).
  * Paterova  _et al._ [2018a] A. Paterova, H. Yang, C. An, D. Kalashnikov, and L. Krivitsky, New Journal of Physics 20, 043015 (2018a).
  * Wang  _et al._ [1991] L. Wang, X. Zou, and L. Mandel, Physical Review A 44, 4614 (1991).
  * Paterova and Krivitsky [2020] A. V. Paterova and L. A. Krivitsky, Light: Science & Applications 9, 82 (2020).
  * Kuznetsov  _et al._ [2020] K. A. Kuznetsov, E. I. Malkova, R. V. Zakharov, O. V. Tikhonova, and G. K. Kitaeva, Physical Review A 101, 053843 (2020).
  * Paterova  _et al._ [2017] A. Paterova, S. Lung, D. Kalashnikov, and L. Krivitsky, Rep 7, 42608 (2017).
  * Kutas  _et al._ [2021] M. Kutas, B. Haase, J. Klier, D. Molter, and G. von Freymann, Optica 8, 438 (2021).
  * Burlakov  _et al._ [1998] A. Burlakov, S. Kulik, A. Penin, and M. Chekhova, Journal of Experimental and Theoretical Physics 86, 1090 (1998).
  * Kulik  _et al._ [2004] S. Kulik, G. Maslennikov, S. Merkulova, A. Penin, L. Radchenko, and V. Krasheninnikov, Journal of Experimental and Theoretical Physics 98, 31 (2004).
  * Klyshko [1993] D. Klyshko, JETP 104, 2676 (1993).
  * Lemos  _et al._ [2014] G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, and A. Zeilinger, Nature 512, 409 (2014).
  * Pearce  _et al._ [2023] E. Pearce, N. R. Gemmell, J. Flórez, J. Ding, R. F. Oulton, A. S. Clark, and C. C. Phillips, Optics Continuum 2, 2386 (2023).
  * Paterova  _et al._ [2020a] A. V. Paterova, H. Yang, Z. S. Toa, and L. A. Krivitsky, Applied Physics Letters 117 (2020a).
  * Paterova  _et al._ [2020b] A. V. Paterova, S. M. Maniam, H. Yang, G. Grenci, and L. A. Krivitsky, Science Advances 6, eabd0460 (2020b).
  * Yang  _et al._ [2023] Y. Yang, H. Liang, X. Xu, L. Zhang, S. Zhu, and X.-s. Ma, npj Quantum Information 9, 2 (2023).
  * Santos  _et al._ [2022] E. A. Santos, T. Pertsch, F. Setzpfandt, and S. Saravi, Physical review letters 128, 173601 (2022).
  * Paterova  _et al._ [2018b] A. V. Paterova, H. Yang, C. An, D. A. Kalashnikov, and L. A. Krivitsky, Quantum Science and Technology 3, 025008 (2018b).
  * Vallés  _et al._ [2018] A. Vallés, G. Jiménez, L. J. Salazar-Serrano, and J. P. Torres, Physical Review A 97, 023824 (2018).
  * Vanselow  _et al._ [2020] A. Vanselow, P. Kaufmann, I. Zorin, B. Heise, H. M. Chrzanowski, and S. Ramelow, Optica 7, 1729 (2020).
  * Machado  _et al._ [2020] G. J. Machado, G. Frascella, J. P. Torres, and M. V. Chekhova, Applied Physics Letters 117 (2020).
  * Grayson  _et al._ [1994] T. Grayson, J. Torgerson, and G. Barbosa, Physical Review A 49, 626 (1994).
  * Paterova  _et al._ [2019] A. Paterova, H. Yang, C. An, D. Kalashnikov, and L. Krivitsky, Optics express 27, 2589 (2019).
  * Fuenzalida  _et al._ [2024] J. Fuenzalida, J. Kysela, K. Dovzhik, G. B. Lemos, A. Hochrainer, M. Lahiri, and A. Zeilinger, Physical Review A 109, 022413 (2024).
  * Herzog  _et al._ [1995] T. J. Herzog, P. G. Kwiat, H. Weinfurter, and A. Zeilinger, Physical Review Letters 75, 3034 (1995).
  * Gemmell  _et al._ [2023] N. R. Gemmell, Y. Ma, E. Pearce, J. Florez, O. Czerwinski, M. Kim, R. F. Oulton, A. S. Clark, and C. C. Phillips, arXiv preprint arXiv:2311.13421 (2023).
  * Kumar  _et al._ [2021] M. Kumar, P. Kumar, A. Vega, M. A. Weissflog, T. Pertsch, and F. Setzpfandt, Applied Physics Letters 119 (2021).
  * Wilson  _et al._ [2020] D. J. Wilson, K. Schneider, S. Hönl, M. Anderson, Y. Baumgartner, L. Czornomaz, T. J. Kippenberg, and P. Seidler, Nature Photonics 14, 57 (2020).
  * Bieszczad  _et al._ [2021] G. Bieszczad, S. Gogler, and J. Świderski, Opto-Electronics Review 29 (2021).
  * Hagen [2022] N. Hagen, Optical Engineering 61, 070902 (2022).
  * Krishnamoorthy  _et al._ [2023] H. N. Krishnamoorthy, A. M. Dubrovkin, G. Adamo, and C. Soci, Chemical Reviews 123, 4416 (2023).
  * Spektor  _et al._ [2015] G. Spektor, A. David, G. Bartal, M. Orenstein, and A. Hayat, Optics Express 23, 32759 (2015).
  * Flament  _et al._ [2002] S. Flament, L. Warsito, L. Mechin, C. Gunther, K. Kawano, G. Kong, C. RA, S. Kahl, K. SI, and G. AM, Transactions of the Magnetics Society of Japan 2, 214 (2002).
  * Adachi  _et al._ [2002] N. Adachi, T. Okuda, V. Denysenkov, A. Jalali-Roudsar, and A. Grishin, Journal of magnetism and magnetic materials 242, 775 (2002).
  * Lopez-Mago and Novotny [2012] D. Lopez-Mago and L. Novotny, Physical Review A 86, 023820 (2012).
  * Sellmeier [1872] W. Sellmeier, Annalen der Physik 223, 525 (1872).


# Spatial deterministic wave forecasting for nonlinear sea–states

M. Galvagno  Centre for Mathematical Sciences, University of Plymouth, PL4 8AA
Plymouth, UK  D. Eeltink  Department of Mechanical Engineering, Massachusetts
Institute of Technology, Cambridge MA 02139-4307, USA  Department of
Engineering Science, University of Oxford, Oxford OX1 3PJ, UK  R. Stuhlmeier
<EMAIL_ADDRESS>Centre for Mathematical Sciences, University
of Plymouth, PL4 8AA Plymouth, UK

###### Abstract

We derive a simple algebraic form of the nonlinear wavenumber correction of
surface gravity waves in deep water, based on temporal measurements of the
water surface and the spatial Zakharov equation. This allows us to formulate
an improvement over linear deterministic wave forecasting with no additional
computational cost. Our new formulation is used to forecast both synthetically
generated as well as experimentally measured seas, and shows marked
improvements over the linear theory.

##  1 Introduction

The goal of deterministic wave forecasting is to determine what waves will
arrive at a distant point, or at some future time, based on spatial or
temporal measurements of the sea surface. To any beachgoer observing the
erratic nature of surface water waves this may seem an impossible task,
recalling Lord Rayleigh’s famous statement that “the basic law of the seaway
is the lack of any law.”

However, with advances in the remote sensing of the sea surface, and attendant
increases in computational power, the problem of deterministic forecasting of
water waves has become tractable, and has attracted significant recent
interest.

In the present work we set out to develop the theoretical basis of a
forecasting methodology that incorporates weakly nonlinear corrections to the
dispersion relation up to third order. While deep water waves undergo (nearly)
resonant interactions at third and higher orders, it is our goal to capture
only the corrected dispersion while neglecting the slow energy exchange
between wave modes. In order to compare with experiments, we consider spatial
evolution, and build our theory upon the spatial Zakharov equation developed
by Shemer et al [26] and employed in numerous subsequent studies [25, 14, 27].

Measurements of the free surface readily yield Fourier amplitudes, which form
the basis of the linear description of surface waves. Indeed, these Fourier
amplitudes can be used to construct simple and efficient linear forecasts,
which have been used in practical tests of deterministic forecasting systems
by Hilmer & Thornhill [10], Kusters et al [16], Al-Ani et al [2], and others.
The same Fourier amplitudes are the foundation of weakly nonlinear approaches,
where corrections to the linear description are sought as perturbations in the
(small) wave slope.

Weakly nonlinear approaches to wave forecasting include those based on PDEs
like the nonlinear Schrödinger equation and its modifications, derived under
assumptions of narrow bandwidth and used by numerous authors, including
Trulsen [32], Simanesew et al [28], Klein et al [15] and others.
Alternatively, the well developed higher-order spectral method (HOS) [7, 34]
presents an attractive computational technique which has gained much recent
attention by the deterministic forecasting community, see e.g. [35, 24, 4, 9,
18].

While they better capture the evolution of real waves, the principal drawback
of these approaches lies in an increased computational cost compared to linear
forecasting. For practical applications, forecasts are needed on scales of
seconds or minutes and tens or hundreds of meters, so speed is of the essence.
Our approach is to extract the correct third-order nonlinear dispersion and
include it as an essentially algebraic correction in the linear forecast. This
is computationally trivial, but we will show that it yields significant
advantages over the purely linear approach.

In what follows, we first review fundamental theory, including linear
forecasting, in section 2.1. Section 2.2 introduces the spatial Zakharov
equation, and contains the main theoretical results. Section 3 applies linear
and nonlinear forecasting methods to synthetically generated seas, simulated
in a numerical wave flume using HOS. Section 4 presents comparisons with
experimental measurements. Finally, section 5 presents a discussion of the
results and some concluding remarks.

##  2 Fundamental theory

Assuming unidirectional propagation of long-crested waves, we may write the
free surface elevation as $\eta(x,t)$ where $x$ is space and $t$ time. In
order to prepare a spatial forecast the sea surface must be measured at a
fixed location $x=x_{0}$ at $N$ times $t_{0},\,t_{1},\ldots t_{N-1}.$ It is
simplest to assume time intervals $\Delta t=T/N$ so that $t_{n}=nT/N,$ but
non-uniformly sampled data can be resampled using interpolation. This leads to
a record

$y_{0}=\eta(x_{0},t_{0}),\,y_{1}=\eta(x_{0},t_{1}),\,\ldots,y_{N-1}=\eta(x_{0},t_{N-1}).$

Taking the discrete-time Fourier transform of the sequence
$y_{0},\,y_{1},\ldots,y_{N-1}$ we find

$Y_{j}=\sum_{n=0}^{N-1}y_{n}\exp(-i2\pi
jn/N)=\sum_{n=0}^{N-1}y_{n}\exp(-it_{n}\omega_{j}),$ (1)

where $\omega_{j}=2\pi j/T.$ Note that $Y(j)=Y(j+N)$ so the Fourier
coefficients are $N$–periodic. The inverse transform is given by

$y_{n}=\frac{1}{N}\sum_{m=0}^{N-1}Y_{m}\exp(2\pi imn/N),$ (2)

which can be transformed into a continuous description using
$nT/N\longrightarrow t$:

$y(t)=\frac{1}{N}\sum_{m=0}^{N-1}Y_{m}\exp(i\omega_{m}t).$ (3)

Finally, this can be reformulated for $N$ even as

$y(t)=\frac{Y_{0}}{N}+\frac{1}{N}\sum_{m=1}^{N/2-1}\left[Y_{m}\exp(i\omega_{m}t)+Y^{*}_{m}\exp(-i\omega_{m}t)\right].$
(4)

Note that $y(t)=y(t+T).$ The term $Y_{0}/N=\frac{1}{N}\sum_{n=0}^{N-1}y_{n}$
is the mean elevation of the sampled points.

###  2.1 Linear forecasting

In the linear theory of water waves there is a one-to-one correspondence
between positive wavenumbers $k\in\mathbb{R}^{+}$ and positive frequency
$\omega\in\mathbb{R}^{+}$, given by the dispersion relation

$\omega^{2}=gk\tanh(kd),$

where $g$ is the acceleration of gravity and $d$ is the (constant) water
depth. For deep water $d\rightarrow\infty$ this dispersion relation reduces to
the simpler expression $\omega^{2}=gk.$ This correspondence allows for a
linear forecast to be constructed from the $N$ samples captured in (4). By
stipulating that a wave with measured frequency $\omega_{m}$ has wavenumber
$k_{m}=\omega_{m}^{2}/g$, it is immediately possible to write:

$\zeta_{L}(x,t)=\frac{Y_{0}}{N}+\frac{1}{N}\sum_{m=1}^{N/2-1}\left[Y_{m}\exp(i(k_{m}x-\omega_{m}t))+Y^{*}_{m}\exp(-i(k_{m}x-\omega_{m}t))\right].$
(5)

The waves accounted for in the forecast then have frequencies between
$\omega_{1}=\frac{2\pi}{T}$ and $\omega_{N/2-1}=\frac{2\pi(N/2-1)}{T}.$ The
energy associated with a given frequency moves at the group velocity, defined
as

$c_{g}:=\frac{d\omega(k)}{dk},$

with the simple form in deep water $c_{g}=0.5{g}/{\omega}.$

For a given measurement, the longest waves of interest $\omega_{L}$ will
travel fastest, and the shortest waves $\omega_{S}$ slowest (note that
practically $\omega_{L}$ may not be $\omega_{1},$ nor $\omega_{S}$ be
$\omega_{N/2-1},$ as there may be negligible energy associated with the
longest or shortest waves that can be theoretically resolved). This leads to
the concept of a predictable region in $(x,t)$ as shown in figure 1.

Figure 1: Predictable region (grey shaded area) based on measurements at $x=0$
in $[t_{0},t_{1}].$

The thick lines in figure 1 show the group velocities $c_{g,L}$ and $c_{g,S}$
of the longest and shortest waves $\omega_{L}$ and $\omega_{S},$ respectively.
Thinner lines in between these indicate the group velocities of waves of
length intermediate between $\omega_{L}$ and $\omega_{S}.$ For a measurement
at $x=0$ over time $t=[t_{0},t_{1}],$ all the waves in the shaded region
originate in the measurement domain, and are therefore predictable. The only
exceptions are waves longer than $\omega_{L}$ or shorter than $\omega_{S}$
that may encroach from $t<t_{0}$ or $t>t_{1},$ and are not accounted for in
the forecast.

###  2.2 Nonlinear forecasting

####  2.2.1 The spatial Zakharov equation

The discussion in the preceding section 2.1 is relevant only for waves of
small steepness, such that linear wave theory may be employed. Weakly
nonlinear theory (to third order) makes for dramatic changes to the dispersion
relation of waves in deep water, and complicates the forecast problem
considerably.

The theoretical basis for our nonlinear forecast will be the spatial Zakharov
equation developed in the early 2000s by Shemer et al [26]. This takes the
form

$\displaystyle ic_{g}\frac{\partial B(x,\omega)}{\partial x}=$
$\displaystyle\iiint
T(k,k_{1},k_{2},k_{3})B^{*}(x,\omega_{1})B(x,\omega_{2})B(x,\omega_{3})$
$\displaystyle\cdot\exp(-i(k+k_{1}-k_{2}-k_{3})x)\delta(\omega_{1}+\omega_{2}-\omega_{3}-\omega_{4})d\omega_{1}d\omega_{2}d\omega_{3}.$
(6)

where $c_{g}$ denotes the deep-water, linear group velocity. This equation can
be discretised as follows:

$ic_{g,j}\frac{dB_{j}(x)}{dx}=\sum_{l,m,n}T_{jlmn}B_{l}^{*}B_{m}B_{n}\exp(-i(k_{i}+k_{j}-k_{m}-k_{n})x)\delta(\omega_{j}+\omega_{l}-\omega_{m}-\omega_{n}),$
(7)

where $B_{i}=B(\omega_{i},x),$ and we abbreviate by $T_{jlmn}$ the kernel
$T(k_{j},k_{l},k_{m},k_{n})$ of the Zakharov equation. In (7) the function
$\delta$ is the ordinary Kronecker delta function. To extract the effect of
nonlinear dispersion we follow the procedure outlined by Stuhlmeier &
Stiassnie [30] for the conventional temporal Zakharov equation (see e.g.
Chapter 14 of Mei et al [21] for background).

We write the complex amplitude $B_{j}(x)$ as $|B_{j}|\exp(i\arg B_{j}),$ where
both magnitude and argument may depend on $x.$ Separating into real and
imaginary parts leads to:

$\displaystyle c_{g,j}\frac{d|B_{j}|}{dx}=-\sum
T_{jlmn}\delta_{jl}^{mn}|B_{l}||B_{m}||B_{n}|\sin(\theta_{jlmn}),$ (8)
$\displaystyle-c_{g,j}|B_{j}|\frac{d\arg(B_{j})}{dx}=\sum
T_{jlmn}\delta_{jl}^{mn}|B_{l}||B_{m}||B_{n}|\cos(\theta_{jlmn}),$ (9)

with

$\displaystyle\theta_{jlmn}=\Delta x+\arg B_{j}+\arg B_{l}-\arg B_{m}-\arg
B_{n},$
$\displaystyle\delta_{jl}^{mn}=\delta(\omega_{j}+\omega_{l}-\omega_{m}-\omega_{n}),\text{
and }$ $\displaystyle\Delta=k_{j}+k_{l}-k_{m}-k_{n}.$

Assuming that there is negligible evolution of the amplitudes, so that the
$|B_{j}|$’s may be replaced by their initial values $|B_{j}(0)|,$ we rewrite

$-c_{g,j}\frac{d}{dx}(\arg
B_{j})=\frac{1}{|B_{j}|}\left(\sum_{l}e_{lj}|B_{l}|^{2}|B_{j}|T_{jljl}+\sum_{l}\sum_{m\neq
j}\sum_{n\neq
j}T_{jlmn}\delta_{jl}^{mn}|B_{l}||B_{m}||B_{n}|\cos(\theta_{jlmn})\right)$
(10)

and, neglecting the second term on the right-hand side (which captures only
exactly resonant quartets) integrate:

$-c_{g,j}(\arg B_{j})=\sum_{l}e_{lj}|B_{l}|^{2}T_{jljl}x+\arg B_{j}(0),$ (11)

where $e_{np}=1$ for $n=p$ and $e_{np}=2$ for $n\neq p.$ The kernels of the
Zakharov equation reduce for two unidirectional waves to

$T(k,k_{1},k,k_{1})=\begin{cases}\frac{k_{x}k_{1x}^{2}}{4\pi^{2}}&\text{ for
}k_{1x}<k_{x},\\\ \frac{k_{x}^{2}k_{1x}}{4\pi^{2}}&\text{ for }k_{1x}\geq
k_{x}.\end{cases}$ (12)

This leads to a correction for the wavenumber:

$K_{n}=k_{n}-\frac{1}{c_{g,n}}\sum_{l}e_{ln}|B_{l}(0)|^{2}T_{lnln},$ (13)

which is the counterpart to the well-known Stokes’ correction to the
frequency. The effects of (weak) nonlinearity are thus to decrease the
wavenumber by an amount of $O(\epsilon^{2})$ compared to the linear theory. We
will explore the effect of this wavenumber correction on two explicit
solutions below, and see that it also impacts the predictable region discussed
in section 2.1 above.

####  2.2.2 Explicit solutions to the spatial Zakharov equation

The spatial Zakharov equation (6) can be easily solved in two special cases: a
single mode, or two modes. The former corresponds to the spatial evolution of
the well-known Stokes’ wave [29], and the latter to the spatial evolution of
the third-order two-wave system first considered by Longuet-Higgins & Phillips
[20]. Either of these cases trivially fulfil the resonance condition, since
$\delta(\omega_{a}+\omega_{a}-\omega_{a}-\omega_{a})=1$ and
$\delta(\omega_{a}+\omega_{b}-\omega_{a}-\omega_{b})=1.$ Because the viewpoint
of wavenumber correction (rather than frequency correction) is somewhat
unusual in water waves, it is instructive to consider these solutions.

In case of a single wave $\omega_{j},$ the spatial Zakharov equation (with
$T_{jjjj}$ abbreviated by $T_{j}$) becomes

$ic_{g,j}\frac{dB_{j}(x)}{dx}=T_{j}|B_{j}(x)|^{2}B_{j}(x),$ (14)

which admits the constant amplitude solution

$B_{j}(x)=A_{j}e^{-iA_{j}^{2}T_{j}x/c_{g,j}}.$ (15)

If two waves are present, say $\omega_{a}$ and $\omega_{b},$ the spatial ZE
becomes the coupled system

$\displaystyle ic_{g,a}\frac{dB_{a}}{dx}$
$\displaystyle=T_{a}|B_{a}|^{2}B_{a}+2T_{ab}|B_{b}|^{2}B_{a},$ (16)
$\displaystyle ic_{g,b}\frac{dB_{b}}{dx}$
$\displaystyle=T_{b}|B_{b}|^{2}B_{b}+2T_{ab}|B_{a}|^{2}B_{b},$ (17)

where the symmetry of the kernel $T_{abab}=T_{abba}$ has been used to simplify
the expressions and $T_{abab}$ has been abbreviated by $T_{ab}.$ Again, this
system admits a solution with constant amplitudes $A_{a}$ and $A_{b},$

$\displaystyle B_{a}(x)$
$\displaystyle=A_{a}\exp(-i(T_{a}A_{a}^{2}+2T_{ab}A_{b}^{2})x/c_{g,a}),$ (18)
$\displaystyle B_{b}(x)$
$\displaystyle=A_{b}\exp(-i(T_{b}A_{b}^{2}+2T_{ab}A_{a}^{2})x/c_{g,b}).$ (19)

The relationship between the complex amplitudes and the leading order free
surface elevation is given by

$\eta(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left(\frac{\omega}{2g}\right)^{1/2}\left[B(x,\omega)\exp(i(k(\omega)x-\omega
t))+\text{c.c.}\right]d\omega.$ (20)

Here “c.c.” stands for the complex conjugate of the preceding expression. For
a single mode $B(x,\omega)=B_{0}(x)\delta(\omega-\omega_{0})$ we find, using
(15),

$\eta(x,t)=\frac{1}{\pi}\left(\frac{\omega_{0}}{2g}\right)^{1/2}A_{0}\cos(([k(\omega_{0})-A_{0}^{2}T_{0000}x/c_{g,0}]x-\omega_{0}t)).$
(21)

This is a simple sinusoidal wave with a wavenumber altered due to the effect
of the weakly nonlinear dispersion relation. Normalising the constant
amplitude via

$A_{0}=\pi a_{0}\left(\frac{2g}{\omega_{0}}\right)^{1/2},$

this reduces to

$\eta(x,t)=a_{0}\cos(k_{0}[1-a_{0}^{2}k_{0}^{2}]x-\omega_{0}t).$

An exactly analogous procedure for two waves ($k_{a}<k_{b}$) yields the free
surface

$\eta(x,t)=a_{a}\cos\left(k_{a}\left[1-a_{a}^{2}k_{a}^{2}-2a_{b}^{2}k_{a}^{3/2}k_{b}^{1/2}\right]x-\omega_{a}t\right)+a_{b}\cos\left(k_{b}\left[1-a_{b}^{2}k_{b}^{2}-2a_{a}^{2}k_{a}^{3/2}k_{b}^{1/2}\right]x-\omega_{b}t\right).$

This shows clearly that the dispersion of one mode $\omega_{a}$ is influenced
both by its own steepness $a_{a}$ and also by that of the second mode
$\omega_{b}.$ It is easy to verify that these wavenumber corrections can be
obtained from equation (13), making liberal use of (12) to simplify the
kernels.

####  2.2.3 Nonlinear forecasts with wavenumber correction

The wavenumber correction (13) gives rise to a simple improved forecast

$\zeta_{N}(x,t)=\frac{Y_{0}}{N}+\frac{1}{N}\sum_{m=1}^{N/2-1}\left[Y_{m}\exp(i(K_{m}x-\omega_{m}t))+Y^{*}_{m}\exp(-i(K_{m}x-\omega_{m}t))\right],$
(22)

where $K_{n}$ denotes the corrected wavenumber. This is otherwise cosmetically
identical to the linear forecast (5), a kinship which demonstrates the
advantages of the formulation. Indeed, the only additional computational cost
consists of calculating the $K_{n}$’s from the Fourier amplitudes via simple
algebra.

In addition, the change in wavenumber leads to a change in phase and group
velocities, which become $c_{p,j}=\omega_{j}/K_{j}$ and
$c_{g,j}=d\omega_{j}/dK_{j},$ respectively. The latter of these must be
evaluated numerically in practice. The nonlinear corrected velocities are
somewhat larger, effectively because of a smaller denominator ($K_{j}$ is
smaller than the linear wavenumber $k_{j},$ and this effect is more pronounced
for shorter waves than for longer waves), which means that the spatial
predictable zone increases in size slightly when nonlinear corrections are
taken into account. The extent of this increase depends on both the
frequencies and the amplitudes present in the measured sea.

##  3 Deterministic forecasting of synthetic seas

To test our two deterministic forecasting methods it is necessary to generate
a wave field, measure it at a point, and compare the forecast with other
measurements. This procedure can be undertaken either in a wave flume or
computationally, and we devote this section to the latter. The advantage of
synthetic seas lies in the ease of tuning inputs, and the ability to easily
produce many realisations of different cases. To this end, we employ the open
source HOS-NWT code [8] which implements the high-order spectral method (HOS).

Figure 2: Diagram of the HOS-NWT numerical wave flume configuration.
$P_{1}$–$P_{4}$ denote the probes at which measurements are taken. Right-going
waves are generated in water of 3 m depth, and fully absorbed by a numerical
beach at the right-hand side of the flume.

Throughout we use a numerical wave flume to generate purely unidirectional,
random wave fields initialised by JONSWAP spectra with peak frequency
$f_{p}=1.3$ Hz, values of significant wave–height $H_{s}$ between 0.02 – 0.04
m, and peak-sharpening parameters $\gamma=1,\,3.3$ and 7. The numerical wave
flume is 15 m long and 3 m deep, with probes located at 3, 6, 9 and 12 m, and
has a fully absorbing beach, as depicted in figure 2. HOS is used to model the
nonlinear propagation of waves along the flume, and includes nonlinear
dispersion, (near) resonant energy exchange, and the effects of bound modes –
considerably more physics than captured by either of our simple forecasts.

Figure 3: The nonlinear predictable zone for a JONSWAP spectrum with
$H_{s}=0.03$ m, $\gamma=3.3$ generated in the 15 m long HOS-NWT numerical wave
flume. Waves are generated at the wavemaker ($x=0$ m) between $t=0$ and
$t=T=72$ s.

To produce a forecast, measurements from probe $P_{1}$ are sampled, the
Fourier amplitudes extracted via FFT, and inserted into either (5) or (22).
The resulting forecasts can then be compared to the measured time series at
probes $P_{2}$, $P_{3}$, and $P_{4}$. Figure 3 depicts the nonlinear
predictable zone based on an example for $H_{s}=0.03$ and $\gamma=3.3.$ Waves
are generated by the wavemaker at $x=0$ m and propagate along the 15 m long
numerical flume. The longest waves resolved are $\omega_{L}=3.1\,\text{rad}/s$
and the shortest $\omega_{S}=25.1\,\text{rad}/s,$ and the associated nonlinear
group velocities determine the edges of the predictable zone. The data at
$P_{1}$ is then used to generate forecasts as shown in figure 4. The
horizontal (time) axis has a different starting point for each panel of figure
4, reflecting the narrowing of the predictable region seen in figure 3 (for
convenience we have here set $t=0$ s as the beginning of the record at
$P_{1},$ rather than the start of the wavemaker).

Figure 4: Comparison of linear ($\zeta_{L}$, red curve) and nonlinear
forecasts ($\zeta_{N}$, yellow curve) with measurements ($P_{i},$ blue curve)
for a numerically generated JONSWAP spectrum with $H_{s}=0.03$ and
$\gamma=3.3.$ The horizontal axis shows time (a section through the
nonlinearly corrected predictable region), while the vertical shows free
surface elevation. (Top panel) comparison with measurements at probe $P_{2}.$
(Middle panel) comparison with measurements at probe $P_{3}.$ (Bottom panel)
comparison with measurements at probe $P_{4}.$

Figure 4 shows excellent agreement between measurement and both forecasts for
the closest probe $P_{2},$ while the forecast grows progressively less
accurate as distance from the measurement point increases. However, the
nonlinear forecast $\zeta_{N}$ remains considerably closer to the measured
data at $P_{4},$ some ten peak wavelengths from the initial measurement.

In order to obtain a measure of the aggregate quality of a forecast it is
useful to compare the fit over several realisations. Table 1 considers the
mean correlation (denoted $\bar{\rho}$) between a linear or nonlinear forecast
produced from probe 1 and the predictable portion of the measured time series
at probe $n,$ for $n=2,\,3,$ or 4. Here and elsewhere the predictable region
is calculated based on the nonlinear group velocities; it is therefore
strictly slightly larger than the comparable predictable zone based on linear
group velocities, see section 2.2.3. In each case 20 realisations of a given
JONSWAP spectrum with random, uniformly distributed phases are considered, and
the forecasting procedure applied to each realisation in turn, before taking
the arithmetic mean of the linear correlations.

| $\gamma=1$ | $\gamma=3.3$ | $\gamma=7$
---|---|---|---
$H_{s}$ | 0.02 | 0.03 | 0.04 | 0.02 | 0.03 | 0.04 | 0.02 | 0.03 | 0.04
$\bar{\rho}_{1,2}^{L}$ | 0.9663 | 0.8498 | 0.7204 | 0.9802 | 0.9015 | 0.8012 | 0.9823 | 0.9274 | 0.8382
$\bar{\rho}_{1,2}^{N}$ | 0.9883 | 0.9432 | 0.8914 | 0.9933 | 0.9637 | 0.9297 | 0.9928 | 0.9704 | 0.9368
$\bar{\rho}_{1,3}^{L}$ | 0.8895 | 0.6375 | 0.4327 | 0.9373 | 0.7419 | 0.5323 | 0.9449 | 0.7930 | 0.6014
$\bar{\rho}_{1,3}^{N}$ | 0.9665 | 0.8736 | 0.7898 | 0.9838 | 0.9227 | 0.8391 | 0.9834 | 0.9311 | 0.8630
$\bar{\rho}_{1,4}^{L}$ | 0.8096 | 0.4499 | 0.2574 | 0.8656 | 0.5604 | 0.3176 | 0.9036 | 0.6314 | 0.3667
$\bar{\rho}_{1,4}^{N}$ | 0.9530 | 0.8071 | 0.7057 | 0.9660 | 0.8708 | 0.7392 | 0.9762 | 0.8700 | 0.7693
Table 1: Mean correlations $\bar{\rho}$ over 20 realisations between synthetic
simulated seas and forecasts based on linear (superscript $L$) and nonlinear
(superscript $N$) forecasts. Subscripts $i,j$ denote that measurements from
probe $i$ are used to forecast probe $j.$ Significant wave height $H_{s}$
ranges from 0.02 to 0.04, and peak-sharpening parameters $\gamma=1,\,3.3$ and
7.

As anticipated from figure 4, the average correlation between measurement and
forecast decreases with distance along the flume. For higher wave steepness,
nonlinear dispersion is well captured by the corrected forecast $\zeta_{N}$
(see rows $\rho^{N}_{i,j}$), as evidenced by the good agreement for moderate
distances. For the steepest waves and longest propagation distances the
forecast quality degrades markedly, nevertheless the nonlinear forecast
retains a clear advantage. Indeed, for $H_{s}=0.3$ and 0.04 and over all
values of $\gamma,$ the nonlinear forecast can predict twice as far (6 m vs 3
m) as the linear forecast with the same average accuracy. It is also
interesting to note that prediction is consistently easier for narrower
spectra ($\gamma=3.3,\,7),$ with accuracy of both linear and nonlinear
forecasts increasing at a given distance as $\gamma$ increases.

##  4 Deterministic forecasting in a wave flume

Figure 5: Diagram of the experimental configuration. $P_{1}$–$P_{4}$ denote
the probes at which measurements are taken. Right-going waves are generated in
water of 0.8 m depth.

To assess the accuracy of our new nonlinear forecasting approach, we also
compare with experimental data from the 40 m long, 2.7 m wide flume at
IRPHE/Pytheas Aix Marseille University. Data is taken from four probes placed
at distances of 3.79, 6.64, 11.63, and 16.11 m from a piston wave maker in
water of depth $d=0.8$ m, as shown in figure 5. As above, measurements from
probe 1 will supply the data necessary to produce a forecast, which is then
compared with data from probes 2–4. We will consider three cases: (J1) a
JONSWAP spectrum with $f_{p}=1.10$ Hz, $\gamma=3.3$ and $H_{s}=0.01$ m, (J2) a
steeper JONSWAP spectrum with $f_{p}=1.11$ Hz, $\gamma=3.3$ and $H_{s}=0.04$
m, and (M) a modulated plane wave, consisting of a plane wave with
$f_{p}=1.42$ and slope $ak=0.16,$ and two side bands.

(a) JONSWAP Case J1.

(b) JONSWAP Case J2.

Figure 6: Comparison of linear ($\zeta_{L}$, red curve) and nonlinear
($\zeta_{N}$, yellow curve) forecasts with measured probe data ($P_{i},$ blue
curve) for JONSWAP cases J1 and J2. Measurements taken at probe $P_{1}$ supply
the Fourier amplitudes for the forecasts at probe $P_{2}$ (top panel), probe
$P_{3}$ (middle panel) and probe $P_{4}$ (bottom panel).

Forecasts for the two JONSWAP cases J1 and J2 are depicted in figures 6(a) and
6(b). For mild seas with $H_{s}=0.01$ either linear or nonlinear forecasting
produces excellent agreement with measurements up to probe $P_{3},$ nearly 8 m
(ca. 7 peak wavelengths) away, as seen in figure 6(a). Due to the low
steepness, the nonlinear correction is essentially negligible, and $\zeta_{L}$
is barely distinguishable from $\zeta_{N}.$ Akin to what was observed for
synthetic data generated by HOS in section 3, as the wave steepness is
increased to $H_{s}=0.04$ m the forecasts begin to depart from the measured
data. The difference in linear and nonlinear forecasts is clearer here, with
the quality of the dispersion-corrected forecast $\zeta_{N}$ outstripping the
simple linear case $\zeta_{L}.$ This information is also captured by the
correlation, used in table 2 to provide a measure of forecast quality (here
only over a single experimental realisation).

Figure 7: Comparison of linear ($\zeta_{L}$, red curve) and nonlinear
($\zeta_{N}$, yellow curve) forecasts with measured probe data ($P_{i},$ blue
curve) for modulated plane wave case M. Measurements taken at probe $P_{1}$
supply the Fourier amplitudes for the forecasts at probe $P_{2}$ (top panel,
2.85 m propagation distance), probe $P_{3}$ (second panel, 7.84 m propagation
distance) and probe $P_{4}$ (third panel, 12.32 m propagation distance). The
bottom panel is an enlargement of the region outlined in black in panel 2. The
Fourier amplitude spectrum at the three measurement gauges $P_{2},\,P_{3}$ and
$P_{4}$ is shown adjacent to panels 1–3.

Forecasts for a modulated plane wave are shown in figure 7. This case exhibits
the well-known modulational instability of a degenerate quartet consisting of
a carrier and two side bands, and the side-band growth with propagation
distance can be clearly seen in the insets depicting the Fourier amplitude
spectrum at probes $P_{2},\,P_{3}$ and $P_{4}.$ As the wave-field propagates
along the flume, the side-band amplitudes to grow at the expense of the
carrier, while also influencing the modes’ dispersion.

The nonlinear forecast $\zeta_{N}$ employs only the initial mode amplitudes to
calculate the corrections to the dispersion relation (the procedure described
in section 2.2.1 effectively neglects the energy exchange between modes,
employing the Fourier amplitudes of $P_{1}$ throughout). The effect of this is
clearly visible in the enlarged forecast in figure 7 (bottom panel) and the
data in table 2: the nonlinear forecast $\zeta_{N}$ shows excellent agreement
with the measured phases, but fails to capture the evolving amplitudes at
larger distances. Because the correlation is insensitive to amplitudes (i.e.
the signals $\sin(x)$ and $2\sin(x)$ have correlation $\rho=1$), good
agreement is found for all $\rho_{i,j}^{N}$ in table 2. For yet longer
propagation distances, the coupling between the phase and amplitude evolution
eventually degrades the otherwise good match between the phases in $\zeta_{N}$
and the experimental record.

| Case M | Case J1 | Case J2
---|---|---|---
$\rho_{1,2}^{L}$ | 0.799 | 0.981 | 0.914
$\rho_{1,2}^{N}$ | 0.967 | 0.982 | 0.955
$\rho_{1,3}^{L}$ | 0.307 | 0.957 | 0.675
$\rho_{1,3}^{N}$ | 0.837 | 0.950 | 0.871
$\rho_{1,4}^{L}$ | -0.178 | 0.929 | 0.551
$\rho_{1,4}^{N}$ | 0.663 | 0.925 | 0.854
Table 2: Correlation between linear (superscript $L$) and nonlinear
(superscript $N$) forecast for the three experimental cases M, J1, and J2.
Subscripts $i,j$ denote that measurements from probe $i$ are used to forecast
probe $j.$

##  5 Discussion

We have derived a compact and theoretically simple wavenumber correction from
the spatial Zakharov equation, and demonstrated its utility in simple cases of
wave forecasting from synthetic and experimentally generated waves. We have
seen that this method accurately captures the most important aspects of
nonlinear dispersion, and is an spatial analogue of the techniques developed
by Stuhlmeier & Stiassnie [31] for the temporal forecasting problem. These
corrections have applicability beyond the immediate context of deterministic
forecasting, for example to the so-called Molin lensing effect [22] which has
recently been studied in the context of wave run-up [36].

Deterministic wave forecasting has numerous uses, ranging from the control of
wave energy converters to improve power capture [19, 12], to ship motion
forecasting for maritime operations [2, 17]. Wave prediction theories, such as
that presented in this manuscript, are only one part of this story: ocean
waves must first be properly measured, and considerable work is being
undertaken on X-band marine radar [37, 23]. Measurement errors and sources of
noise must be dealt with efficiently, as discussed recently by Desmars et al
[6], for example by continuous data assimilation and ensemble Kalman filtering
[33]. Subsequently prediction can be accomplished with a wide variety of
propagation techniques.

We have used the fast Fourier transform throughout, and tacitly assumed that
it introduces no appreciable errors into the forecasting methodology. This is
not quite the case, as the sampled water surface is not a strictly periodic
signal, and we have only a finite-length snapshot at each probe location. This
induces a rectangular windowing and results in spectral leakage, recently
addressed in the context of forecasting [11, 1].

It is interesting to observe that in our synthetic forecasts a decrease in
spectral width (by increasing $\gamma$) increases the average accuracy of our
forecasts, as measured by the linear correlation and presented in table 1. Of
the two phenomena associated with cubically nonlinear wave propagation in deep
water, energy exchange is expected to be more significant for a narrow
spectrum, while frequency correction is expected to be less significant. The
former phenomenon is connected to the Benjamin-Feir index introduced by
Janssen (see [13]), connecting spectral width and scale of nonlinearity to the
appearance of modulational instability. The latter is a consequence of the
asymmetry of (13): energy in long waves has a large effect on the dispersion
of short waves, but not vice versa. For a broader spectrum, with energy
distributed among modes further from the spectral peak (especially in higher
frequencies), these dispersion corrections will therefore be more significant
(see [30]).

The extremely narrow and discrete spectrum of the modulated plane wave in
figure 7, allows to track energy transfer clearly, and makes it a popular
laboratory wave, although it is unlikely to be found on the ocean. The energy
exchange associated with the modulational instability [3, 5] drives
significant changes in the spectral amplitudes between one probe and the next,
while both linear and nonlinear forecasts implicitly assume the spectrum
remains unchanged throughout the propagation. Therefore, it is interesting to
see that as expected, the amplitudes of the waves are not well predicted, yet,
the phases are matched very well for the nonlinear prediction.

## Acknowledgements

RS gratefully acknowledges support by EPSRC grant EP/V012770/1 and a QJMAM
Fund grant from the IMA. Research visits to the University of Plymouth by MG
were supported by a QJMAM Fund grant from the IMA. DE acknowledges financial
support from the Swiss National Science Foundation (Fellowship P2GEP2-191480).
The authors thank Christopher Luneau and Hubert Branger from the IRPHE/Pytheas
Aix Marseille University for their help during the experiments.

## References

  * [1] L. Abusedra and M. R. Belmont.  Prediction diagrams for deterministic sea wave prediction and the introduction of the data extension prediction method.  Int. Shipbuild. Prog., 58(1):59–81, 2011.
  * [2] M. Al-Ani, M. Belmont, and J. Christmas.  Sea trial on deterministic sea waves prediction using wave-profiling radar.  Ocean Eng., 207:107297, 2020.
  * [3] T. B. Benjamin and J. E. Feir.  The disintegration of wave trains on deep water Part 1. Theory.  J. Fluid Mech., 27(03):417–430, 1967.
  * [4] E. Blondel, F. Bonnefoy, and P. Ferrant.  Deterministic non-linear wave prediction using probe data.  Ocean Eng., 37(10):913–926, 2010.
  * [5] A. Chabchoub, B. Kibler, C. Finot, G. Millot, M. Onorato, J. Dudley, and A. Babanin.  The nonlinear Schrödinger equation and the propagation of weakly nonlinear waves in optical fibers and on the water surface.  Ann. Phys., 361:490–500, 2015.
  * [6] N. Desmars, F. Bonnefoy, S. T. Grilli, G. Ducrozet, Y. Perignon, C. A. Guérin, and P. Ferrant.  Experimental and numerical assessment of deterministic nonlinear ocean waves prediction algorithms using non-uniformly sampled wave gauges.  Ocean Eng., 212(January):107659, 2020.
  * [7] D. G. Dommermuth and D. K. Yue.  A high-order spectral method for the study of nonlinear gravity waves.  J. Fluid Mech., 184(1987):267–288, 1987.
  * [8] G. Ducrozet, F. Bonnefoy, D. Le Touzé, and P. Ferrant.  A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank.  Eur. J. Mech. B/Fluids, 34:19–34, 2012.
  * [9] C. A. Guérin, N. Desmars, S. T. Grilli, G. Ducrozet, Y. Perignon, and P. Ferrant.  An improved Lagrangian model for the time evolution of nonlinear surface waves.  J. Fluid Mech., 876:527–552, 2019.
  * [10] T. Hilmer and E. Thornhill.  Observations of predictive skill for real-time Deterministic Sea Waves from the WaMoS II.  Ocean. 2015 - MTS/IEEE Washingt., pages 1–7, 2016.
  * [11] T. Hlophe, H. Wolgamot, A. Kurniawan, P. H. Taylor, J. Orszaghova, and S. Draper.  Fast wave-by-wave prediction of weakly nonlinear unidirectional wave fields.  Appl. Ocean Res., 112(March):102695, 2021.
  * [12] M. Huchet, A. Babarit, G. Ducrozet, J. C. Gilloteaux, and P. Ferrant.  Nonlinear deterministic sea wave prediction using instantaneous velocity profiles.  Ocean Eng., 220(November 2020), 2021.
  * [13] P. A. Janssen.  The Interaction of Ocean Waves and Wind.  Cambridge University Press, 2004.
  * [14] E. Kit and L. Shemer.  Spatial versions of the Zakharov and Dysthe evolution equations for deep-water gravity waves.  J. Fluid Mech., 450:201–205, 2002.
  * [15] M. Klein, M. Dudek, G. F. Clauss, S. Ehlers, J. Behrendt, N. Hoffmann, and M. Onorato.  On the deterministic prediction of water waves.  Fluids, 5(1):1–19, 2020.
  * [16] J. G. Kusters, K. L. Cockrell, B. S. Connell, J. P. Rudzinsky, and V. J. Vinciullo.  FutureWaves™: A real-time Ship Motion Forecasting system employing advanced wave-sensing radar.  Ocean. 2016 MTS/IEEE Monterey, OCE 2016, 2016.
  * [17] J. G. Kusters, B. S. Connell, W. M. Milewski, V. J. Vinciullo, and R. Van Dijk.  Wave characterization and timing using doppler radar - Update on the futurewaves™ wave and vessel motion forecasting system.  Ocean. 2019 MTS/IEEE Seattle., (October), 2019.
  * [18] Y. Law, H. Santo, K. Lim, and E. Chan.  Deterministic wave prediction for unidirectional sea-states in real-time using Artificial Neural Network.  Ocean Eng., 195:106722, 2020.
  * [19] G. Li, G. Weiss, M. Mueller, S. Townley, and M. R. Belmont.  Wave energy converter control by wave prediction and dynamic programming.  Renew. Energy, 48:392–403, 2012.
  * [20] M. S. Longuet-Higgins and O. M. Phillips.  Phase velocity effects in tertiary wave interactions.  J. Fluid Mech., 12(03):333–336, 1962.
  * [21] C. C. Mei, M. A. Stiassnie, and D. K.-P. Yue.  Theory and applications of ocean surface waves.  World Scientific Publishing Co., 3rd edition, 2018.
  * [22] B. Molin, F. Remy, O. Kimmoun, and E. Jamois.  The role of tertiary wave interactions in wave-body problems.  J. Fluid Mech., 528:323–354, 2005.
  * [23] M. Previsic, A. Karthikeyan, and D. Lyzenga.  In-Ocean Validation of a Deterministic Sea Wave Prediction ( DSWP ) System leveraging X-Band Radar to Enable Optimal Control in Wave Energy Conversion Systems.  Appl. Ocean Res., 114(July):102784, 2021.
  * [24] Y. Qi, G. Wu, Y. Liu, M.-H. Kim, and D. K. P. Yue.  Nonlinear phase-resolved reconstruction of irregular water waves.  J. Fluid Mech, 838:544–572, 2018.
  * [25] L. Shemer and A. Chernyshova.  Spatial evolution of an initially narrow-banded wave train.  J. Ocean Eng. Mar. Energy, 3(4):333–351, 2017.
  * [26] L. Shemer, H. Jiao, E. Kit, and Y. Agnon.  Evolution of a nonlinear wave field along a tank: experiments and numerical simulations based on the spatial Zakharov equation.  J. Fluid Mech., 427:107–129, 2001.
  * [27] L. Shemer, E. Kit, and H. Jiao.  An experimental and numerical study of the spatial evolution of unidirectional nonlinear water-wave groups.  Phys. Fluids, 14(10):3380–3390, 2002.
  * [28] A. Simanesew, K. Trulsen, H. E. Krogstad, and J. C. Nieto Borge.  Surface wave predictions in weakly nonlinear directional seas.  Appl. Ocean Res., 65:79–89, 2017.
  * [29] G. G. Stokes.  On the theory of oscillatory waves.  Trans. Camb. Phil. Soc, 8:441–455, 1847.
  * [30] R. Stuhlmeier and M. Stiassnie.  Nonlinear dispersion for ocean surface waves.  J. Fluid Mech., 859:49–58, 2019.
  * [31] R. Stuhlmeier and M. Stiassnie.  Deterministic wave forecasting with the Zakharov equation.  J. Fluid Mech., 913:1–22, 2021.
  * [32] K. Trulsen.  Spatial Evolution Of Water Surface Waves.  In Fifth Int. Symp. WAVES 2005, number 127, pages 1–10, 2005.
  * [33] G. Wang and Y. Pan.  Phase-resolved ocean wave forecast with ensemble-based data assimilation.  J. Fluid Mech., 918:A19, 2021.
  * [34] B. J. West, K. A. Brueckner, R. S. Janda, D. M. Milder, and R. L. Milton.  A new numerical method for surface hydrodynamics.  J. Geophys. Res., 92(C11):11803, 1987.
  * [35] G. Wu.  Direct Simulation and Deterministic Prediction of Large-scale Nonlinear Ocean Wave-field.  MIT Ph.D. thesis, (1994):258, 2004.
  * [36] W. Zhao, P. H. Taylor, H. A. Wolgamot, and R. Eatock Taylor.  Amplification of random wave run-up on the front face of a box driven by tertiary wave interactions.  J. Fluid Mech., 869:706–725, 2019.
  * [37] V. Zinchenko, L. Vasilyev, S. O. Halstensen, and Y. Liu.  An improved algorithm for phase-resolved sea surface reconstruction from X-band marine radar images.  J. Ocean Eng. Mar. Energy, 7(1):97–114, 2021.