# Space-Time Quantum Metasurfaces

Wilton J. M. Kort-Kamp  Los Alamos National Laboratory, Los Alamos, NM 87545,
USA  Abul K. Azad  Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Diego A. R. Dalvit∗ Los Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

Metasurfaces are a key photonic platform to manipulate classical light using
sub-wavelength structures with designer optical response. Static metasurfaces
have recently entered the realm of quantum photonics, showing their ability to
tailor nonclassical states of light. We introduce the concept of space-time
quantum metasurfaces for dynamical control of quantum light. We provide
illustrative examples of the impact of spatio-temporally modulated
metasurfaces in quantum photonics, including the creation of frequency-spin-
path hyperentanglement on a single photon and the realization of space-time
asymmetry at the deepest level of the quantum vacuum. Photonic platforms based
on the space-time quantum metasurface concept have the potential to enable
novel functionalities, such as on-demand entanglement generation for quantum
communications, nonreciprocal photon propagation for free-space quantum
isolation, and reconfigurable quantum imaging and sensing.

The generation, manipulation, and detection of nonclassical states of light is
at the heart of quantum photonics. As quantum information can be encoded into
the different degrees of freedom of a single photon, it is highly desirable to
develop photonic platforms that allow to control them while maintaining
quantum coherence. Metasurfaces Kildishev2013 ; Chen2016  have recently
transitioned from the classical to the quantum domain Solntsev2020  and
enabled enhanced light-matter interactions facilitated by the ultrathin
subwavelength nature of their constituent scatterers. Spin-orbital angular
momentum entanglement of a single photon has been demonstrated Stav2018  using
geometrically tailored metasurfaces that induce spin-orbit coupling of light
via the Pancharatnam-Berry phase Bomzon2002 , and multiphoton interferences
and polarization-state quantum reconstruction has also been achieved in a
single geometric phase metasurface Wang2018 . A metasurface-based
interferometer has been demonstrated for generating and probing entangled
photon states Georgi2019 , opening opportunities for implementing quantum
sensing and metrology protocols using metasurface platforms. Embedded quantum
building blocks into arrays of meta-atoms, including quantum dots,
semiconductor emitters, and nitrogen-vacancy centers, result in enhanced
Purcell factors Vaskin2019 , directional lasing Xie2020 , and circularly-
polarized single-photon emission Kan2020 . Recently, quantum metasurfaces
based on atomic arrays have been proposed Bekenstein2020 .

Most demonstrations of metasurfaces in quantum photonics are based on static
meta-atoms whose optical properties are determined by their material
composition and geometrical design that cannot be changed on demand. A few
realizations of active quantum metasurfaces have been reported, for example,
for tuning spontaneous emission from a Mie-resonant dielectric metasurface
using liquid crystals Bohn2018 . However, a fully tailored response requires
quantum metasurfaces that can continuously alter their scattering properties
simultaneously in space and time. At the classical level, spatio-temporally
modulated metasurfaces Shaltout2019  have been shown to provide that higher
degree of control, both by reconfigurable and fully-dynamic tailoring of the
optical response of meta-atoms using analog and digital modulation schemes
Cardin2020 ; Zhang2018 .

Here, we put forward the concept of space-time quantum metasurfaces (STQMs)
for spatio-temporal control of quantum light. In order to highlight the broad
implications of this concept in different areas of quantum science and
technology, we discuss two instances of how STQMs operate both at the single-
photon and virtual-photon levels. We describe STQM-enabled hyperentanglement
Barreiro2005  manipulation of nonclassical states of light and STQM-induced
photon pair generation (Fig. 1) in a process analogue to the dynamical Casimir
effect Moore1970 .

Results:

We introduce the STQM concept by considering the transmission of a single
photon through a metasurface whose meta-atoms are modulated in space and time.
The metasurface is composed of identical anisotropic scatterers (Fig. 2a)
suitably rotated with respect to each other. The combination of anisotropy and
rotation results in circular cross-polarization conversion and a spin-
dependent geometric phase distribution $\Psi({\bf r})$ akin to spin-orbit
coupling. To minimize photon absorption the metasurface is assumed to be
comprised of low-loss dielectric meta-atoms. The spatio-temporal modulation is
modeled as a perturbation of the electric permittivity, $\epsilon({\bf
r},t)=\epsilon_{um}+\Delta\epsilon\cos(\Omega t-\Phi({\bf r}))$, where
$\epsilon_{um}$ is the unmodulated permittivity, $\Delta\epsilon$ the
modulation amplitude, $\Omega$ the modulation frequency, and $\Phi({\bf r})$
is a “synthetic” phase. Such type of modulation has been recently demonstrated
using a heterodyne laser-induced dynamical grating on an amorphous Si
metasurface via the nonlinear Kerr effect Guo2019 , setting a traveling-wave
permittivity perturbation with $\Phi({\bf r})=\bm{\beta}\cdot{\bf r}$
($\bm{\beta}$ is an in-plane momentum “kick”). Note that the geometric phase
is fixed by the design of the metasurface while the synthetic phase is
reconfigurable on-demand.

Figure 1: Conceptual representation of a space-time quantum metasurface. A
spatio-temporal spinning modulation of graphene nanostructures generates
entangled vortex photon pairs out of the quantum vacuum.

STQMs for on-demand entanglement manipulation: The geometry of the dielectric
nanoresonator can be tailored so that it has maximal cross-polarized
transmission (Fig. 2b) and at the same time so that its Mie electric and
magnetic dipolar resonances dominate the optical response of the metasurface.
One can then describe the interaction of each Mie resonator with light using
the effective Hamiltonian $H_{int}=-{\bf p}\cdot{\bf E}-{\bf m}\cdot{\bf B}$
Kuznetsov2016 ; Novotny2007 , where ${\bf p}$ and ${\bf m}$ are the electric
and magnetic dipole operators and ${\bf E}$ and ${\bf B}$ are the local
quantized electromagnetic fields. Higher-order Mie resonances can be neglected
because the transmissivity and reflectivity of the metasurface are well-
described by that of an array of electric and magnetic dipoles corresponding
to the two lowest Mie multipoles. It is convenient to trace over the
nanostructure’s degrees of freedom to express the Hamiltonian only in terms of
photonic modes by relating dipoles and fields via effective electric
$\bm{\alpha}_{E}$ and magnetic $\bm{\alpha}_{M}$ polarizability tensors (see
Supplementary Information for the derivation of the polarizabilities). We show
in Fig. 2c the relevant polarizability components for describing the coupling
with the normally-incident photon and in Fig. 2d the electric field
distributions at the Mie resonance frequencies. The unmodulated Hamiltonian
describing cross-polarized transmission has an effective coupling strength
$\alpha^{(cr)}_{um}(\omega)$ that is a simple combination of the electric and
magnetic polarizabilities (see Methods).

Upon spatio-temporal modulation, the effective polarizabilities adiabatically
follow the harmonic driving because the response times of semiconductors
($<100$ fs for the nonlinear Kerr response time in amorphous Si) are much
faster than THz modulations achievable with all-optical schemes. Hence,

$\alpha^{(cr)}(\omega;{\bf
r},t)=\alpha^{(cr)}_{um}(\omega)+\Delta\alpha^{(cr)}(\omega)\cos(\Omega
t-\Phi({\bf r})).$ (1)

We calculate the polarizability modulation amplitude
$\Delta\alpha^{(cr)}(\omega)$ from the dependency of transmissivity on
permittivity modulation (Fig. 2e). For a $1\%$ permittivity modulation depth
the resulting polarizability change is approximately $20\%$, the increase
originating from the strong dispersion of the unmodulated polarizability close
to the input frequency. The STQM Hamiltonian $H_{1}(t)$ is the sum of the
unmodulated part plus a modulation contribution that annihilates the input
photon and creates a new one with Doppler-shifted frequency and synthetic
phase, in addition to flipping its spin components and adding geometric phases
in the same way as the unmodulated part (see Methods). In this work we
restrict to unitary evolution as photons do not suffer from severe decoherence
problems and absorption is negligible in high-index dielectrics Stav2018 .

Figure 2: Effective polarizabilities of all-dielectric space-time quantum
metasurfaces. (a) Anisotropic amorphous Si Mie nanocross meta-atom with
optimized geometrical parameters for maximal cross-polarization transmission
for a normally-incident $\lambda_{in}=1550$ nm input photon. Parameters are
$L_{1}=950$ nm, $L_{2}=435$ nm, $h=300$ nm, $w=200$ nm, and square unit cell
with period $P=1200$ nm. (b) Co- and cross-polarized reflectivity and
transmissivity for the full metasurface (solid) and electric/magnetic dipole
array (dashed). (c) Real parts of the electric and magnetic polarizabilities
normalized by the meta-atom volume. Solid line is the effective unmodulated
coupling strength for cross-polarized transmission: $\alpha^{(cr)}_{um}\approx
0.6\mu{\rm m}^{3}$ at the input frequency. (d) Electric field distribution for
the two electric and the two magnetic Mie resonances. (e) Polarizability
modulation amplitudes for permittivity modulation depth
$\Delta\epsilon/\epsilon_{um}=1\%$. Solid line is the polarizability
modulation amplitude for cross-polarized transmission:
$\Delta\alpha^{(cr)}/\alpha^{(cr)}_{um}\approx 0.2$ at the input frequency.

When the geometric phase is a linear function of the meta-atoms’ positions it
generates spin-momentum correlations, while a linear synthetic phase creates
momentum-frequency correlations. The two correlations are intertwined through
momentum and the photon evolves into a state that is hyperentangled in spin,
path, and frequency

$\\!\\!|\psi(t)\rangle\\!=\\!\sum_{p,q}[c^{(R)}_{p,q}(t)|\omega_{p};{\bf
k}_{p,q};\\!R\rangle+c^{(L)}_{p,q}(t)|\omega_{p};{\bf k}_{p,-q};\\!L\rangle],$
(2)

where $p$ are integers, $q=0,1$, $R(L)$ denotes right (left) circular
polarization, $\omega_{p}=\omega_{in}+p\Omega$ are harmonics of the input
frequency $\omega_{in}$, ${\bf k}_{p,q}={\bf
k}_{in}+p\bm{\beta}+q\bm{\beta}_{g}$ are in-plane momentum harmonics of the
in-plane input wave-vector ${\bf k}_{in}$, and $\bm{\beta}_{g}$ is the
momentum kick induced by the linear geometric phase. We will denote states in
the first term as $(p,q,R)$ and in the second term as $(p,-q,L)$, highlighting
that the geometric-phase-induced momentum kicks for right- and left-polarized
photons have opposite directions. To calculate the probability amplitudes we
consider a normally-incident single-photon pulse and assume modulation
frequencies and in-plane momentum kicks much smaller than the input frequency
and input wave-vector. Since the dielectric metasurface enables large
polarizability modulation amplitudes for modest permittivity variations, it is
possible for the input photon to transition to multiple frequency/momentum
harmonics during its transit within the metasurface. For input linear
polarization, the transition probabilities to states $(p,q,R)$ and $(p,-q,L)$
are identical and are given by

$|c^{(R/L)}_{p,q}(t)|^{2}=\frac{1}{2}\cos^{2}\Big{(}\frac{\omega_{in}t\alpha^{(cr)}_{um}}{2hP^{2}}\Big{)}J^{2}_{p}\Big{(}\frac{\omega_{in}t\Delta\alpha^{(cr)}}{2hP^{2}}\Big{)}$
(3)

when $p$ and $q$ have the same parity; for opposite parity the cosine is
replaced by a sine. $J_{p}(x)$ is the Bessel function and probabilities for
$\pm p$ are the same (see Supplementary Information for details on the state
evolution).

Figure 3:  Entanglement manipulation with space-time quantum metasurfaces. (a)
Conversion probability into an output photon in frequency harmonic
$\omega_{in}+p\Omega$ and momentum harmonic $p\bm{\beta}+q\bm{\beta}_{g}$
versus polarizability modulation depth for the STQM of Fig. 2. (b) Density
matrices of input (i) and output photon featuring (ii) spin-path entanglement,
(iii-iv) frequency-path entanglement, and (v-vi) frequency-spin-path
hyperentanglement. Larger modulation depths result in more harmonics involved
in the entangled states (iv) and (vi). (c) Population dynamics of geometric-
phase kicked $q=1$ and unkicked $q=0$ states for in-transit photon at $0.1$
modulation depth, showing Rabi oscillations for the fundamental frequency
harmonic. Inset: Rabi dynamics in higher harmonics. Envelopes (dashed black)
are the populations of $p$-harmonics in the absence of geometric phase.

The probability that the output photon is in a given frequency/momentum
harmonic as a function of the modulation depth is shown in Fig. 3a. At zero
modulation, the output has the same frequency as the input and is
approximately an equal superposition of right- and left-polarized geometric-
phase-kicked states, with a small overlap with unkicked states due residual
co-polarized transmission. As the modulation increases, transitions to only
the first few frequency/momentum harmonics occur and a larger amount of the
available Hilbert space is explored at large modulation depths. Figure 3b
depicts the density matrices of the input (panel (i)) and output photons for
different configurations of the STQM, resulting in distinct kinds of quantum
correlations: (ii) Geometric phase with spatio-temporal modulation off, giving
a spin-path entangled output of same frequency as input; (iii-iv) No geometric
phase and spatio-temporal modulation on, resulting in frequency-path entangled
cross-polarized output; (v-vi) Geometric phase with spatio-temporal modulation
on, delivering a frequency-spin-path hyperentangled output. It is possible to
tailor the modulation depth to completely suppress the contribution of a given
harmonic to the output state, as shown in (iv, vi) for the fundamental
frequency. Under temporal modulation only, i.e., null synthetic phase (not
shown), the output photon is unentangled (hyperentangled) in the absence
(presence) of geometric phase. Figure 3c shows the population dynamics of
different harmonics while the photon is in-transit inside the STQM.
Interestingly, the evolution of populations with and without geometric phase
are fundamentally different. Due to spin-orbit coupling the photon undergoes
Rabi-flopping between state pairs $(p,0,R)\leftrightarrow(p,-1,L)$ and
$(p,0,L)\leftrightarrow(p,1,R)$, and this population exchange cannot take
place at zero geometric phase. Unmodulated and modulated polarizabilities
control the time-scales of Rabi and synthetic-phase dynamics, respectively.

When both phase distributions are azimuthally varying, i.e., $\Psi({\bf
r})=\ell_{g}\varphi$, $\Phi({\bf r})=\ell\varphi$ ($\ell_{g}$ and $\ell$
integers), the input photon becomes hyperentangled in frequency, spin, and
orbital angular momentum (OAM) Calvo2006 . The state of the photon can be
written as in Eq. (2) replacing linear momentum harmonics ${\bf k}_{p,q}$ by
OAM harmonics $\ell_{p,q}=p\ell+q\ell_{g}$. Such a rotating synthetic phase
could be implemented, e.g., via a heterodyne laser-induced dynamical grating
with Laguerre-Gauss petal modes Eichler1986 ; Naidoo2012  to generate an all-
optical spinning perturbation of the meta-atoms’ refractive index. As STQMs
offer the possibility to reconfigure the synthetic phase on-demand, the
question naturally arises as to what happens when the synthetic and geometric
phase distributions have utterly different symmetry, for instance one is
linear and the other spinning. It is then necessary to expand one phase in
terms of a mode basis with symmetries appropriate for the other phase, e.g.,
plane waves into cylindrical waves (see Supplementary Information for details
of mixed-phase STQMs). As the synthetic phase creates frequency-path
correlations and the geometric phase spin-OAM correlations, the two
correlations are not intertwined and the STQM does not produce
hyperentanglement but bipartite entanglement between pairs of degrees of
freedom of the single photon. Finally, we mention that all the analysis
presented in this section can be extended to other nonclassical inputs, such
as two-photon Fock states.

STQMs for tailored photon-generation out of quantum vacuum: Space-time quantum
metasurfaces can produce other nonclassical states of light and even induce
nonreciprocity Sounas2017  on quantum vacuum fluctuations. In addition to the
photon-number-conserving Hamiltonian $H_{1}(t)$ discussed above, STQMs couple
to the quantum electromagnetic field via a photon-number-non-conserving
Hamiltonian $H_{2}(t)$ that creates photon pairs out of the quantum vacuum
(see Methods). Their frequencies add up to the modulation frequency,
$\omega+\omega^{\prime}=\Omega$, thereby conserving energy, and this process
is essentially an analogue of the dynamical Casimir effect (DCE) in which an
oscillating boundary parametrically excites virtual into real photons
Dalvit2006 ; Dodonov2010 . Although the mechanical DCE effect has not been
detected because it requires unfeasibly large mechanical oscillation
frequencies, various analogue DCE systems have been demonstrated Wilson2011 ;
Jaskula2012 ; Lahteenmaki2013 ; Vezzoli2019 . STQMs allow for a novel degree
of dynamical control over the quantum vacuum through the synthetic phase: The
scattering matrix for the DCE process Maghrebi2013  becomes asymmetric via the
spatio-temporal modulation, reflecting that Lorentz reciprocity is broken at
the level of quantum vacuum fluctuations.

Figure 4: Steered quantum vacuum. (a) A linear synthetic phase is imprinted on
a metasurface through a traveling-wave modulation and is tuned on demand
(colored arrows) to steer the emitted dynamical Casimir photons. (b) Emission
lobes of one photon for varying momentum kick and fixed (vertical) emission
direction of its twin. (c) Density polar plots of angular emission spectrum
for various $\beta=(0,0.2,0.3,0.38,0.5)\Omega/c$ from left to right. The areas
to the right (left) of the vertical solid line correspond to the angular
emission spectrum of the high- (low-) frequency photon in a pair. Frequencies
are $\omega/\Omega=0.7$ and $\omega^{\prime}/\Omega=0.3$. Shaded zones
correspond to forbidden photon emission directions. Between the two rightmost
panels two special events simultaneously happen: the merge of the emission
“island” of the high-frequency photon with the grazing line and the birth of
forbidden regions for the low-frequency photon. (d) Spherical polar plots for
the same panels as in (c).  Figure 5: Photo-emission rates for various
synthetic phases. (a) Spectral weight function for linear synthetic phase.
Sharp edges of each plateau correspond to the special events of Fig. 4(c).
Inset: Crossings responsible for non-monotonicity in (f). (b) Spectral weight
function for rotating synthetic phase. Solid lines correspond to a finite
radius metasurface ($\Omega R/c=30$) and dashed line is the $\ell=0$ case for
an infinite metasurface. (c) Angular-momentum spectra for finite radius
metasurface for the high-frequency photon. (d) Same as (c) for the low-
frequency photon. (e) Spectral photo-production rate for null synthetic phase
for a graphene-disk STQM. Inset: unmodulated electric polarizability
$\alpha_{um}(\omega)$ (solid) and modulation amplitude $\Delta\alpha(\omega)$
(dashed). (f) Emission rate for linear synthetic phase. The profile on the
left shows the rate at null synthetic phase. The black thick curve joins peaks
of maximal emission and the thin black curve
$c\beta=2\omega_{res}(E_{F})-\Omega$ is its projection on the $\beta-E_{F}$
plane. The rate decreases non-monotonically to zero at $\beta_{max}$. In (e-f)
parameters are: $\Omega/2\pi=10$ THz, $\Delta E_{F}/E_{F}=1\%$,
$n_{MS}=10^{3}\;{\rm mm}^{-2}$, $D=5\;\mu$m, and graphene mobility
$\mu=10^{4}\;{\rm cm}^{2}\;{\rm V}^{-1}\;{\rm s}^{-1}$.

We first consider the case of the linear synthetic phase (Fig. 4a) and set the
geometric phase to zero. Momentum conservation dictates that the emitted
photons must have in-plane momenta that add up to the imprinted kick,
$\bf{k}+\bf{k}^{\prime}=\bm{\beta}$, and the emitted photons are frequency-
path entangled. In Fig. 4b we show the one-photon angular emission
distribution for a fixed propagation direction of its twin, indicating how the
externally imprinted momentum controls the directivity of the emission
process. Figure 4c contains polar plots of the emission distributions for a
given circularly-polarized photon pair (see Methods). In the absence of kick,
the high-frequency photon can be emitted in any azimuthal direction but it has
a maximum polar angle of emission, while no such a constraint exists for the
low-frequency photon. As the magnitude of the momentum kick $\beta$ increases,
the distributions undergo intricate changes. The region of allowed emission
for the first photon gets deformed when the kick is non-zero and at a critical
value of the kick an “island” of emission appears surrounded by a sea of
forbidden emission directions (shaded areas). The island drifts to higher
polar angles until it touches the grazing emission line, starts to shrink in
size, and finally at $\beta_{max}=\Omega/c$ it collapses to a point and the
photon is only emitted parallel to the kick. Far-field emission above that
value of the kick is not possible. Regarding the second photon, its emission
distribution remains mostly unperturbed until two areas of forbidden emission
appear at large polar angles and opposite to the kick direction. The forbidden
region grows until it engulfs its allowed emission region and a second island
forms (not shown). Finally, it ends up being emitted at a grazing angle but in
a direction anti-parallel to the kick. The corresponding spherical plots are
shown in Fig. 4d, with emission profiles resembling cone- (dome-) like shapes
for the high- (low-) frequency photon and become increasingly distorted as the
kick grows. When both photons are emitted with the same frequency, i.e., twin
photons, the emission distribution is disk-shaped and gets elongated in a
direction parallel to the kick as this increases in magnitude (not shown). The
modulation also excites hybrid entangled pairs composed of one photon and one
evanescent surface wave (shaded areas in Fig. 4c), and when
$\beta>\beta_{max}$ only evanescent modes are created and subsequently decay
via non-radiative loss mechanisms.

The two-photon emission rate from an STQM of area $A$ with arbitrary synthetic
phase $\Phi({\bf r})$ is

${\Gamma}_{\Phi}=\frac{An_{MS}^{2}\Omega^{4}}{512\pi^{3}c^{4}}\\!\int_{0}^{\Omega}\\!\\!\\!d\omega|\Delta\alpha(\omega)+\Delta\alpha(\Omega-\omega)|^{2}f_{\Phi}(\omega).$
(4)

The rate scales as the square of the meta-atoms number surface density
$n_{MS}$ indicating coherent emission of photon pairs. Electro-optical
properties of the meta-atoms are contained in the modulated electric
polarizability amplitude $\Delta\alpha(\omega)$. The spectral weight function
$f_{\Phi}$ results from the angular integration of all emission events and is
plotted in Fig. 5a for the case of the linear synthetic phase.
$f_{\beta}(\omega)$ has a central plateau-like form with sharp edges that at
zero kick coalesce into a single logarithmic integrable divergency at the
center of the spectrum and corresponds to the emission of twin photons
MaiaNeto1996 . As the kick grows, the plateau becomes lower and at the maximum
allowed kick the spectral weight function vanishes.

STQMs can affect the quantum vacuum in more exotic ways, e.g., a modulation
with a spinning synthetic phase “stirs” the vacuum (Fig. 1) and induces
angular momentum nonreciprocity Sounas2013  at the level of quantum
fluctuations. The rotating modulation generates vortex photon pairs that carry
angular momenta satisfying $m+m^{\prime}=\ell$. For $\ell\neq 0$ the average
of the Poynting vector over all possible emission events results in a single
vortex line along the synthetic spinning axis. Photon pairs are frequency-
angular momentum entangled and their quantum correlations could be probed
using photo-coincidence detection and techniques based on angular momentum
sorting of light Berkhout2010 ; Mirhosseini2013 . The spectral weight function
$f_{\ell}(\omega)$ is reported in Fig. 5b for a finite-radius metasurface,
showing plateau-like structures as in Fig. 5a but without the sharp features
on the edges, and with decreasing height as the spinning grows. There is a
drastic but subtle difference between the two spectral weight functions
$f_{\beta}(\omega)$ and $f_{\ell}(\omega)$ that is not apparent in the plots:
The former vanishes beyond the finite kick threshold $\beta_{max}$, while
there is no finite spinning threshold for the latter. Figures 5c-d show the
angular momentum spectra of high- and low-frequency photons in an emitted pair
(see Methods). When the STQM does not imprint any spinning, the spectra are
symmetric around the peak at $m=0$, with oppositely twisted photons in each
emitted pair. When spinning is imprinted, the two spectra are related as
$f_{\ell}(m,\omega)=f_{\ell}(\ell-m,\Omega-\omega)$ and the spectrum for the
high- (low-) frequency photon is centered around $m=\ell$ ($m=0$). This is the
angular-momentum equivalent of asymmetric linear momentum emission in Fig. 4d.

Photo-emission rates can be boosted with suitable functional meta-atoms, such
as atomically-thin nanostructures made of plasmonic materials that can support
highly localized plasmons Yu2017 ; Abajo2015 ; Muniz2020  and enable large
electric polarizabilities conducive to enhanced coupling of the STQM with the
quantum vacuum. As an example, we consider a STQM based on graphene disks
whose Fermi energy $E_{F}$ is spatio-temporally modulated (Fig. 1). Changing
the Fermi energy it is possible to tune the plasmonic resonances into the DCE
spectral range and to modify the electric polarizability modulation amplitude
(inset of Fig. 5e). Furthermore, the use of ultra-high mobility graphene
samples minimizes photon absorption and substantially enhances photo-
production rates. Figure 5e depicts the spectral rate for a STQM for null
synthetic phase at selected Fermi energies, featuring Lorentzian peaks at
complementary frequencies. For high-Q resonances the emission rate for
arbitrary synthetic phase can be approximated as

${\Gamma}_{\Phi}\approx
g\Omega\;(An_{MS}^{2}D^{6}\omega^{4}_{res}/c^{4})\;f_{\Phi;res}\Big{(}\frac{\Delta
E_{F}}{E_{F}}\Big{)}^{2}\Big{(}\frac{\Omega}{\gamma}\Big{)}^{3}$ (5)

Here, $g$ is a numerical factor determined by the plasmon eigenmode, $D$ the
disk diameter, $\omega_{res}$ is the plasmonic resonance frequency,
$f_{\Phi;res}$ the spectral weight on resonance, $\Delta E_{F}/E_{F}$ the
Fermi energy modulation depth, and $\Omega/\gamma\gg 1$ with $\gamma$ the
scattering rate of graphene. (see Supplementary Information for the derivation
of the polarizability and emission rate). Figure 5f shows the emission rate
for the linear synthetic phase as a function of momentum kick and Fermi
energy. Giant photon-pair production rates on the order of $10^{12}$
photons$/{\rm cm}^{2}{\rm s}$ are obtained at low-THz driving frequencies and
modest modulation depths. Conventional electrical doping may not allow to
reach the large Fermi energies where the rate is maximized, but it suffices
for exploring the lower Fermi energy region where photon-pairs are already
produced in giant numbers. A heterodyne dynamical grating based on ultrafast
all-optical THz modulation of graphene conductivity Tasolamprou2019  could
enable giant and steered photo-pair emission out of the quantum vacuum.
Finally, we note that electro-optical ultrafast on-chip graphene modulators
Li2014 ; Phare2015 ; Kovacevic2018  could potentially be employed to
independently bias different STQM graphene pixels with designer temporal
delays to implement complex synthetic phases.

Discussion: Metasurfaces are crossing the classical-quantum divide to offer
novel possibilities for flat quantum optics and photonics. On the quantum
side, they can become an enabler platform for generating and manipulating
nonclassical states of light in real time. We uncovered a key property of
space-time quantum metasurfaces relevant for potential applications: On-demand
reconfiguration of the synthetic phase allows dynamically tunable quantum
correlations, enabling to tailor the nature of entanglement depending on the
symmetry properties of both geometric and synthetic phases. We also
illustrated a second key property of space-time quantum metasurfaces with
fundamental relevance: Lorentz nonreciprocity at the deepest level of vacuum
fluctuations is attained through joint space and time modulations of optical
properties and can be interpreted as an asymmetric quantum vacuum.

Spatio-temporally modulated quantum metasurfaces have the potential to become
a flexible photonic platform for generating nonclassical light with designer
spatial and spectral shapes, for on-demand manipulation of entanglement for
free-space communications, and for reconfigurable sensing and imaging systems.
Conversion efficiencies into specific frequency, linear momentum, or orbital
angular momentum harmonics for selective quantum information encoding could be
enhanced through advanced modulation protocols. Incorporation of quantum
matter building blocks into space-time metasurfaces may further expand the
possibilities afforded by the proposed platform. As such, space-time quantum
metasurfaces can provide breakthrough advances in quantum photonics.

Methods: The effective polarizability tensors $\bm{\alpha}_{E}(\omega)$ and
$\bm{\alpha}_{M}(\omega)$ of the dielectric nanostructure are obtained using a
Cartesian multipole expansion Evlyukhin2013  of the full-wave simulated
electromagnetic field under a plane wave excitation, and computing ratios of
the resulting Mie electric and magnetic dipoles to the incident field at the
nanostructure’s center. The Hamiltonian for the all-dielectric STQM in cross-
polarized transmission is

$\displaystyle
H_{1}(t)=-\sum_{j,\gamma,\gamma^{\prime}}\;[\alpha^{(cr)}_{um}(\omega)+\Delta\alpha^{(cr)}(\omega)\cos(\Omega
t-\Phi_{\\!j})]$ (6) $\displaystyle\times
A^{*}_{\gamma;j}A_{\gamma^{\prime};j}e^{i(\omega-\omega^{\prime})t}[e^{i\Psi_{\\!j}}a^{\dagger}_{\gamma,R}a_{\gamma^{\prime},L}\\!+\\!e^{-i\Psi_{\\!j}}a^{\dagger}_{\gamma,L}a_{\gamma^{\prime},R}]\\!+\\!h.c.$

The sums are over all meta-atoms and field modes, the geometric $\Psi({\bf
r})$ and synthetic $\Phi({\bf r})$ phase distributions are evaluated at the
position of the meta-atoms, $\Omega$ is the modulation frequency,
$A_{\gamma}$, $A_{\gamma^{\prime}}$ are spatial modes, and
$a_{\gamma^{\prime},L/R}$ and $a^{\dagger}_{\gamma,R/L}$ are annihilation and
creation operators of circularly polarized photons. The unmodulated coupling
strength is

$\alpha^{(cr)}_{um}(\omega)\\!=\\!{\rm
Re}[\alpha_{E,xx}(\omega)+\alpha_{M,yy}(\omega)-\alpha_{E,yy}(\omega)-\alpha_{M,xx}(\omega)]$

and $\Delta\alpha^{(cr)}(\omega)$ is the modulation coupling strength obtained
by replacing in the above equation each effective polarizability by its
respective modulation amplitude.

The Hamiltonian for the all-plasmonic STQM for two-photon emission is

$\displaystyle H_{2}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sum_{j,\gamma,\gamma^{\prime}}\sum_{\lambda,\lambda^{\prime}}\;[\Delta\alpha(\omega-\Omega)+\Delta\alpha(\omega^{\prime}-\Omega)]$
$\displaystyle\times$ $\displaystyle
A^{*}_{\gamma;j}A^{*}_{\gamma^{\prime};j}\;e^{i\Phi_{j}}e^{i(\omega+\omega^{\prime}-\Omega)t}\;a^{\dagger}_{\gamma,\lambda}a^{\dagger}_{\gamma^{\prime},\lambda^{\prime}}+h.c.$

where we neglected multiscattering between meta-atoms Holloway2005 .
$\lambda,\lambda^{\prime}$ are polarization states of the two photons and
$\Delta\alpha(\omega)$ is the modulated electric polarizability amplitude of
the meta-atom computed with the plasmon wavefunction formalism Yu2017 . For a
graphene disk of diameter $D$ with a high-Q localized bright-mode plasmonic
resonance

$\Delta\alpha(\omega)\\!\approx\\!\frac{\pi^{2}a_{1}^{2}\alpha_{fs}cD^{2}\Delta
E_{F}}{512\hbar\Omega^{2}}\frac{(\gamma/2\Omega)^{2}}{[((\omega\\!-\\!\omega_{res})/\Omega)^{2}\\!+\\!(\gamma/2\Omega)^{2}]^{2}}.$

$\omega_{res}(E_{F})=(\alpha_{fs}cE_{F}/\pi|\xi_{1}|\hbar D)^{1/2}$ is the
resonance frequency of the lowest bright-mode, $E_{F}$ and $\Delta E_{F}$ are
the Fermi energy and its modulation amplitude, $\alpha_{fs}$ is the fine
structure constant, $\gamma=ev^{2}_{F}/E_{F}\mu$ the scattering rate of
graphene, $v_{F}$ the Fermi velocity, and $\mu$ the mobility. The eigenmode
coefficients $a_{1}=6.1$ and $\xi_{1}=-0.072$ determine
$g=5\pi^{4}a_{1}^{4}\xi_{1}^{2}/2(512)^{3}$ in Eq. (5). Emission rates are
computed with time-dependent perturbation theory. Non-paraxial quantization of
the electromagnetic field with angular momentum is employed for the STQM with
rotating synthetic phase Enk1994 . The spectral weight functions for the
linear and spinning synthetic phases are respectively decomposed into in-plane
linear momentum $f_{\bm{\beta}}({\bf k},\omega)$ and angular momentum
${f}_{\ell}(m,\omega)$ spectra

$\displaystyle f_{\bm{\beta}}(\omega)$ $\displaystyle=$ $\displaystyle\int
d{\bf k}\;(c/\Omega)^{2}\;(1-|c{\bf k}/\omega|^{2})^{-1/2}f_{\bm{\beta}}({\bf
k},\omega),$ (8) $\displaystyle f_{\ell}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{m}{f}_{\ell}(m,\omega).$

Explicit expressions for these spectra can be found in the Supplementary
Information.

Acknowledgements: This work was supported by the DARPA QUEST program. We are
grateful to A. Efimov, M. Julian, C. Lewis, M. Lucero, and A. Manjavacas for
discussions.

Author Contributions: D.A.R.D. and W. K.-K. conducted the theory work and A.
K. A. analyzed experimental feasibility. All authors discussed the findings
and contributed to writing the paper.

Competing Interests: The authors declare no competing interests.

∗Correspondence<EMAIL_ADDRESS>
## References

  * (1) A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, Planar photonics with metasurfaces, Science 339, 1232009 (2013).
  * (2) H.-T. Chen, A. J. Taylor, and N. Yu, A review of metasurfaces: physics and applications, Rep. Prog. Phys. 79, 076401 (2016).
  * (3) A. S. Solntsev, G. S. Agarwal, and Y. S. Kivshar, Metasurfaces for quantum photonics, arXiv:2007.14722.
  * (4) T. Stav, A. Faerman, E. Maguid, D. Oren, V. Kleiner, E. Hasman, and M. Segev, Quantum entanglement of the spin and orbital angular momentum of photons using metamaterials, Science 361, 1101-1104 (2018).
  * (5) Z. Bomzon, G. Biener, V. Kleiner, and E. Hasman, Space-variant Pancharatnam–Berry phase optical elements with computer-generated subwavelength gratings, Opt. Lett. 27, 1141-1143 (2002).
  * (6) K. Wang, J. G. Titchener, S. S. Kruk, L. Xu, H.-P. Chung, M. Parry, I. I. Kravchenko, Y.-H. Chen, A. S. Solntsev, Y. S. Kivshar, D. N. Neshev, and A. A. Sukhorukov, Quantum metasurface for multiphoton interference and state reconstruction, Science 361, 1104–1108 (2018).
  * (7) P. Georgi, M. Massaro, K.-H. Luo, B. Sain, N, Montaut, H. Herrmann, T. Weiss, G. Li, C. Silberhorn, and T. Zentgraf, Metasurface interferometry toward quantum sensors, Light: Science & Applications 8, 70 (2019).
  * (8) A. Vaskin, R. Kolkowskia, A. F. Koenderink, and I. Staude, Light-emitting metasurfaces, Nanophotonics 8, 1151–1198 (2019).
  * (9) Y.-Y. Xie, P.-N. Ni, Q.-H. Wang, Q. Kan, G. Briere, P.-P. Chen, Z.-Z. Zhao, A. Delga, H.-R. Ren, H.-D. Chen, C. Xu, and P. Genevet, Metasurface-integrated vertical cavity surface-emitting lasers for programmable directional lasing emissions, Nat. Nanotech. 15, 125–130 (2020).
  * (10) Y. Kan, S. K. H. Andersen, F. Ding, S. Kumar, C. Zhao, and S. I. Bozhevolnyi, Metasurface-enabled generation of circularly polarized single photons, Adv. Mater. 32, 1907832 (2020).
  * (11) R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin, and M. D. Lukin, Quantum metasurfaces with atom arrays, Nat. Phys. 16, 676–681 (2020).
  * (12) J. Bohn, T. Bucher, K. E. Chong, A. Komar, D.-Y. Choi, D. N. Neshev, Y. S. Kivshar, T. Pertsch, and I. Staude, Active tuning of spontaneous emission by Mie-resonant dielectric metasurfaces, Nano Lett. 18, 3461–3465 (2018).
  * (13) A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, Spatiotemporal light control with active metasurfaces, Science 364, 648 (2019).
  * (14) A. E. Cardin, S. R. Silva, S. R. Vardeny, W. J. Padilla, A. Saxena, A. J. Taylor, W. J. M. Kort-Kamp, H.-T. Chen, D. A. R. Dalvit, and A. K. Azad, Surface-wave-assisted nonreciprocity in spatio-temporally modulated metasurfaces, Nat. Commun. 11, 1469 (2020).
  * (15) L. Zhang, X. Q. Chen, S. Liu, Q. Zhang, J. Zhao, J. Y. Dai, G. D. Bai, X. Wan, Q. Cheng, G. Castaldi, V. Galdi, and T. J. Cui, Space-time-coding digital metasurfaces, Nat. Commun. 9, 4334 (2018).
  * (16) J. T. Barreiro, N. K. Langford, N. A. Peters, and P. G. Kwiat, Generation of Hyperentangled Photons Pairs, Phys. Rev. Lett. 95, 260501 (2005).
  * (17) G. T. Moore, Quantum theory of the electromagnetic field in a variable‐length one‐dimensional cavity, Journal of Mathematical Physics 11, 2679 (1970).
  * (18) X. Guo, Y. Ding, Y. Duan, and X. Ni, Nonreciprocal metasurface with space–time phase modulation, Light: Science & Applications 8, 123 (2019).
  * (19) A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S Kivshar, and B. Luk’yanchuk, Optically resonant dielectric nanostructures, Science 354, aag2472 (2016).
  * (20) L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University Press, New York, 2007).
  * (21) G. F. Calvo, A. Picón, and E. Bagan, Quantum field theory of photons with orbital angular momentum, Phys. Rev. A 73, 013805 (2006).
  * (22) H. J. Eichler, P. Günter, and D. H. Pohl, Laser induced dynamical gratings (Springer-Verlag, Heidelberg, 1986).
  * (23) D. Naidoo, K. Ait-Ameur, M. Brunel, and A. Forbes, Intra-cavity generation of superpositions of Laguerre–Gaussian beams, Appl. Phys. B, 106, 683–690 (2012).
  * (24) D. L. Sounas and A. Alù, Non-reciprocal photonics based on time modulation, Nat. Photonics 1, 774–783 (2017).
  * (25) D. A. R. Dalvit, P. A. Maia Neto, and F. D. Mazzitelli, in Casimir Physics, Lecture Notes in Physics 834 (eds. D. Dalvit, P. Milonni, D. Roberts, and F. Da Rosa), (Springer, Heidelberg, 2006).
  * (26) V. V. Dodonov, Current status of the dynamical Casimir effect, Physica Scripta 82, 038105 (2010).
  * (27) C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, and P. Delsing, Observation of the dynamical Casimir effect in a superconducting circuit, Nature 479, 376–379 (2011).
  * (28) J.-C. Jaskula, G. B. Partridge, M. Bonneau, R. Lopes, J. Ruaudel, D. Boiron, and C. I. Westbrook, Acoustic Analog to the Dynamical Casimir Effect in a Bose-Einstein Condensate, Phys. Rev. Lett. 109, 220401 (2012).
  * (29) P. Lahteenmaki, G. S. Paraoanu, J. Hassel, and P. J. Hakonen, Dynamical Casimir effect in a Josephson metamaterial, Proc. Nat. Acad. Sci. USA 110, 4234-4238 (2013).
  * (30) S. Vezzoli, A. Mussot, N. Westerberg, A. Kudlinski, H. D. Saleh, A. Prain, F. Biancalana, E. Lantz, and D. Faccio, Optical analogue of the dynamical Casimir effect in a dispersion-oscillating fibre, Commun. Phys. 2, 84 (2019).
  * (31) M. F. Maghrebi, R. Golestanian, and M. Kardar, Scattering approach to the dynamical Casimir effect, Phys. Rev. D 87, 025016 (2013).
  * (32) P. A. Maia Neto and L. A. S. Machado, Quantum radiation generated by a moving mirror in free space, Phys. Rev. A 54, 3420 (1996).
  * (33) D. L. Sounas, C. Caloz, and A. Alù, Giant non-reciprocity at the subwavelength scale using angular momentum-biased metamaterials, Nat. Commun. 4, 2407 (2013).
  * (34) G. C. Berkhout, M. P. Lavery, J. Courtial, M. W. Beijersbergen, and M. J. Padgett, Efficient Sorting of Orbital Angular Momentum States of Light, Phys. Rev. Lett. 105, 153601 (2010).
  * (35) M. Mirhosseini, M. Malik, Z. Shi, and R. W. Boyd, Efficient separation of the orbital angular momentum eigenstates of light, Nat. Commun. 4, 2783 (2013).
  * (36) R. Yu, J. D. Cox, J. R. M. Saavedra, and F. J. García de Abajo, Analytical modeling of graphene plasmons, ACS Photonics 4, 3106 (2017).
  * (37) F. J. García de Abajo and A. Manjavacas, Plasmonics in atomically thin materials, Faraday Discuss. 178, 87-1073548 (2015).
  * (38) Y. Muniz, A. Manjavacas, C. Farina, D. A. R. Dalvit, and W. J. M. Kort-Kamp, Two-Photon Spontaneous Emission in Atomically Thin Plasmonic Nanostructures, Phys. Rev. Lett. 125, 033601 (2020).
  * (39) A. C. Tasolamprou, A. D. Koulouklidis, C. Daskalaki, C. P. Mavidis, G. Kenanakis, G. Deligeorgis, Z. Viskadourakis, P. Kuzhir, S. Tzortzakis, M. Kafesaki, E. N. Economou, and C. M. Soukoulis, Experimental demonstration of ultrafast THz modulation in a graphene-based thin film absorber through negative photoinduced conductivity, ACS Photonics 6, 720–727 (2019).
  * (40) W. Li, B. Chen, C. Meng, W. Fang, Y. Xiao, X. Li, Z. Hu, Y. Xu, L. Tong, H. Wang, W. Liu, J. Bao, and Y. R. Shen, Ultrafast all-optical graphene modulator, Nano Lett. 14, 955–959 (2014).
  * (41) C. T. Phare, Y.-H. D. Lee, J. Cardenas, and M. Lipson, Graphene electro-optic modulator with 30 GHz bandwidth, Nat. Photonics 9, 511–514 (2015).
  * (42) G. Kovacevic, C. Phare, S. Y. Set, M. Lipson, and S. Yamashita, Ultra-high-speed graphene optical modulator design based on tight field confinement in a slot waveguide, Appl. Phys. Express 11, 065102 (2018).
  * (43) A. B. Evlyukhin, C. Reinhardt, E. Evlyukhin, and B. N. Chichkov, Multipole analysis of light scattering by arbitrary-shaped nanoparticles on a plane surface, J. Opt. Soc. Am. B 30, 2589 (2013).
  * (44) S. J. Van Enk and G. Nienhuis, Commutation rules and eigenvalues of spin and orbital angular momentum of radiation fields, J. Mod. Opt. 41, 963-977 (1994).
  * (45) C. L. Holloway, M. A. Mohamed, E. F. Kuester, and A. Dienstfrey, Reflection and transmission properties of a metafilm: with an application to a controllable surface composed of resonant particles, IEEE Transactions on Electromagnetic Compatibility 47, 853-865 (2005).


solution to this quandry, although it is clearly related to the discussion in
Section 8. We do not elaborate further, since this discrepancy is of only
marginal statistical significance.

Perhaps our most significant result — strong evidence for source evolution
with redshift, i.e. $n>0$ — is most strongly disfavoured by $p_{n}$, such that
when this is removed, we prefer a strongly evolving population. This,
predictably, is related to the low number of bursts observed by Parkes/Mb
relative to ASKAP/FE compared to expectations discussed above, which alone
argues against source evolution. Removing any of the three FRB surveys, or
$p_{\rm snr}$, simply results in mildly less evidence for $n>0$ — there is no
single result strongly favouring $n>0$. We thus conclude that our conclusion
on $n$ robust, especially considering uncertainties discussed in Appendix A.1.

Our host galaxy parameters, $\mu_{\rm host}$ and $\sigma_{\rm host}$ — recall
these are in log DM space — are most affected by the highest- and lowest-
sensitivity surveys, ASKAP/FE and Parkes/Mb. Without these observations, very
low values of $\mu_{\rm host}$ (and hence very large values of $\sigma_{\rm
host}$) become possible. This validates the use of DM distributions only when
studying these parameters. Neither $p_{\rm snr}$ nor $p_{n}$ have a
significant effect.


# Improved calculation of the second-order ocean Doppler spectrum for sea
state inversion

Charles-Antoine Guérin

###### Abstract

We propose a simple method, based on an original change of variables, for the
fast and accurate calculation of the second-order ocean Doppler spectrum that
describes the sea echo of High-Frequency radars. A byproduct of the technique
is the derivation of an improved weighting function which can be used for the
inversion of the main sea state parameters. For this we revisit Barrick’s
method for the estimation of the significant wave height and the mean period
from the ocean Doppler spectrum. On the basis of numerical simulations we show
that a better estimation of these parameters can be reached but necessitates a
preliminary bias correction that depends only on the radar frequency. A second
consequence of our formulation is the derivation of a simple yet analytical
nonlinear approximation of the second-order ocean Doppler spectrum when the
Doppler frequency is larger than the Bragg frequency. This opens new
perspective for the inversion of directional wave spectra from High-Frequency
radar measurements.

###### Index Terms:

High-Frequency radar, Doppler spectrum, sea state inversion

##  I Introduction

High-Frequency radars have been used for half a century as an efficient tool
for measuring surface currents and waves in the coastal region (see e.g. [1]
for a recent review). In this range of radio frequencies (3-30 MHz), the
backscattered Doppler spectrum from the sea surface is accurately described by
the second-order electromagnetic and hydrodynamic perturbation theory, whose
complete equations were published by Barrick half a century ago [2].
Hasselmann [3] first suggested that the continuous component of the ocean
Doppler spectrum is essentially a replica of the ocean wave spectrum
translated by the Bragg frequency and thus could be used for sea state
inversion. As many followers have shown, this is only a coarse approximation
and estimating the wavenumber spectrum has both mathematical and practical
limitations (e.g. [4]). There is no unique relationship between the ocean
spectrum and the Doppler spectrum, so that the inversion can only be performed
with some additional constraints. The numerical inversion of the Doppler
spectrum, using either Barrick’s linearized equations [5, 6] or their fully
nonlinear expression [7], is computationally intensive. It requires solving a
large number of direct problems involving the tricky calculation of a double
integral with singular kernel and constrained variables. In this paper, we
propose an improved method for the fast and accurate calculation of the
second-order ocean Doppler spectrum based on an appropriate, original change
of variables (section II). For technical reasons we limit this analysis to the
backscattering configuration and to the deep-water gravity wave dispersion
relation for the hydrodynamic processes at play, a double assumption which
corresponds to the majority of practical situations. Equipped with this new
formulation, we address the classical problem of estimating the main sea state
parameters from the integral properties of the Doppler spectrum, following the
method first proposed by Barrick [8, 9]. To do this, we introduce a function,
referred to as the “Zeta” function, which has the same structure as the
second-order Doppler spectrum but is free of the coupling coefficient. We show
that the Zeta function can, in principle, be used for the estimation of the
significant wave height (SWH) and the mean period. However, numerical
simulations using a typical ocean wave spectrum model show that the estimation
of these 2 parameters must be corrected by a frequency-dependent
multiplicative and additive bias, respectively (section III). Transferring
these results to the actual Doppler spectrum requires the use of a weighting
function, as first proposed by Barrick, to account for the influence of the
coupling coefficient. In the light of our improved formulation we derive a new
weighting function (Section IV) to approximate the Zeta function, resulting in
a more accurate estimation of the sea state parameters. Our last result
concerns the direct estimation of the ocean wave spectrum. For Doppler
frequencies higher than the Bragg frequency we obtain, as a consequence of the
new integral formulation, a simple analytical approximation which is verified
to be extremely accurate (Section V).

##  II Evaluation of the second-order spectrum

We consider the classical idealized problem of a High-Frequency radar
measuring the backscattered ocean Doppler spectrum in vertical polarization at
grazing incidence in the limit of an infinite sea surface patch. We denote
$f_{0}$ the radar frequency, $k_{0}=2\pi f_{0}/c_{0}$ the associated
electromagnetic wavenumber, where $c_{0}=3.10^{8}$ m/s is the speed of light
in vacuum, and $\boldsymbol{k}_{0}$ the horizontal electromagnetic wave vector
of the incident electric field. As it is well known in ocean radar remote
sensing, the so-called Bragg wave vector
$\boldsymbol{k}_{B}=-2\boldsymbol{k}_{0}$, Bragg wavenumber $k_{B}=2k_{0}$ and
Bragg frequency $\omega_{B}=\sqrt{gk_{B}}$ play a prominent role in the
expression of the backscattered field. The complete expression of the ocean
Doppler cross-section appeared for the first time in [2]. It is given by a sum
of two terms, $\sigma(\omega)=\sigma_{1}(\omega)+\sigma_{2}(\omega)$, where
$\sigma_{1}$ contains the first-order Bragg peaks:

$\sigma_{1}(\omega)={\cal N\mit}\sum_{n_{1}=\pm
1}S_{d}(n_{1}\boldsymbol{k}_{B})\delta(\omega-\omega_{B}),$ (1)

and $\sigma_{2}$ is a second-order, continuous component:

$\begin{split}\sigma_{2}(\omega)&={\cal N\mit}\sum_{n_{1}=\pm
1}\sum_{n_{2}=\pm 1}\int
S_{d}(n_{1}\boldsymbol{k}_{1})S_{d}(n_{2}\boldsymbol{k}_{2})\\\
&\times\left|\Gamma\right|^{2}\delta(\omega-
n_{1}\omega_{1}-n_{2}\omega_{2})d\boldsymbol{k}_{1}.\end{split}$ (2)

where ${\cal N\mit}=2^{6}\pi K_{0}^{4}$ is a normalization constant and
$S_{d}(\boldsymbol{k})$ the directional ocean spectrum. The wave vectors
entering in this integral are related by the Bragg resonance condition and the
associated frequencies are given by the deep water gravity waves dispersion
relation:

$\boldsymbol{k}_{1}+\boldsymbol{k}_{2}=\boldsymbol{k}_{B},\
\omega_{1}=\sqrt{gk_{1}},\ \omega_{2}=\sqrt{gk_{2}},$ (3)

where $g=9.81\ ms^{-2}$ is the gravity constant. The kernel $\Gamma$ can be
decomposed into a hydrodynamic kernel $\Gamma_{H}$ and an EM kernel
$\Gamma_{EM}$:

$\Gamma(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\omega,\omega_{B})=\Gamma_{H}+\Gamma_{EM}$
(4)

with

$\Gamma_{H}=-\frac{\rm{i}}{2}\left(k_{1}+k_{2}-\frac{(k_{1}k_{2}-\boldsymbol{k_{1}}\cdot\boldsymbol{k_{2}})({\omega}^{2}+\omega_{B}^{2})}{n_{1}n_{2}\sqrt{k_{1}k_{2}}({\omega}^{2}-\omega_{B}^{2})}\right)$
(5)

and

$\Gamma_{EM}=\frac{1}{2}\frac{(\boldsymbol{k_{1}}\cdot\boldsymbol{k_{B}})(\boldsymbol{k_{2}}\cdot\boldsymbol{k_{B}})-2\boldsymbol{k_{1}}\cdot\boldsymbol{k_{2}}}{\sqrt{\boldsymbol{k_{1}}\cdot\boldsymbol{k_{2}}}-\Delta/2},$
(6)

where $\Delta=0.011-0.012i$ is an impedance term that takes into account the
finite conductivity of the sea surface in the High-Frequency domain. Each of
the four components in (2) corresponding to the different sign combinations
$(n_{1},n_{2})$ addresses a different frequency domain of the Doppler
spectrum. The combined geometrical constraints
$\omega=n_{1}\omega_{1}+n_{2}\omega_{2}$ and
$\boldsymbol{k}_{1}+\boldsymbol{k}_{2}=\boldsymbol{k}_{B}$ impose
$n_{1}=n_{2}=1$ if $\omega\geq\omega_{B}$, $n_{1}=n_{2}=-1$ if
$\omega\leq-\omega_{B}$ and $n_{1}=-n_{2}=\pm 1$ if
$\left|\omega\right|\leq\omega_{B}$. The different combinations of signs thus
cover the four distinct regions delimited by the first order Bragg lines.

The evaluation of the second-order integral (2) is difficult due to the delta
function constraint with respect to the $\omega_{1}$ and $\omega_{2}$
variables, which has been analyzed in detail in [10]. For a prescribed value
of $\omega$, the admissible values of $\boldsymbol{k}_{1}$ follow closed
contour lines in the $(k_{1x},k_{1y})$ plane. The contours have 2 disjoint
components for $\omega<\sqrt{2}\omega_{B}$ that merge into a single one for
$\omega>\sqrt{2}\omega_{B}$. The intensification of the contour lines at the
transition produces a singularity of the function $\sigma_{2}(\omega)$ at
$\omega=\sqrt{2}\omega_{B}$. A second singularity is caused by the EM kernel,
that has a sharp maximum for $\boldsymbol{k_{1}}\cdot\boldsymbol{k_{2}}=0$
(this would be a mathematical singularity without the finite impedance term).
The locus of this singularity is a circle which is found tangential to the
contour line $\omega=2^{3/4}\simeq 1.68\omega_{B}$; this results in a
secondary, sharp maximum at this frequency. Considering the finite integration
time $T_{i}$ and the resulting finite frequency resolution
$\Delta\omega=2\pi/T_{i}$, the theoretical Doppler spectrum is actually
smeared by some function $\Phi(\omega)$ of width $\Delta\omega$ [11]. This
amounts to replace the Dirac function $\delta(\omega-
n_{1}\omega_{1}-n_{2}\omega_{2})$ in the integrand of (2) by $\Phi(\omega-
n_{1}\omega_{1}-n_{2}\omega_{2})$ where the integration domain is now two-
dimensional with no further constraint on the integration variable
$\boldsymbol{k}_{1}$. This allows for a straightforward computation of
$\sigma_{2}$ but requires a tight 2D sampling of the k-space that is
numerically not efficient. An alternative computation can be obtained by
reducing the two-dimensional integral (2) to a single one by eliminating the
delta function through an appropriate change of variables. The standard way to
do this [10] is to use polar coordinates $(k_{1},\theta)$ for the wave vector
$\boldsymbol{k}_{1}$ and to operate a change a variable
$(k_{1},\theta)\to(\omega_{1}+\omega_{2},\theta)$. However, the change of
variable is only implicit and requires an extra numerical inversion for each
value of $\omega$ to calculate the new integrand. We here propose an
alternative change of variable that is more efficient and fully explicit. It
works only with deep water gravity waves, which follow the simple dispersion
relation $\omega^{2}=gk$. This is, however, the case of main interest for HF
radars. To simplify the calculation we use normalized frequencies and wave
numbers:

$\begin{split}\kappa_{1}=\boldsymbol{k}_{1}/k_{B},&\
\kappa_{2}=\boldsymbol{k}_{2}/k_{B},\
\boldsymbol{\kappa}_{B}=\boldsymbol{k}_{B}/k_{B}\\\
\nu_{1}=\omega_{1}/\omega_{B}&\ \nu_{2}=\omega_{2}/\omega_{B}\
\nu=\omega/\omega_{B}\end{split}$ (7)

This is relevant because the kernel $\Gamma$ depends only on the reduced
variables and the Bragg wavenumber.

$\Gamma(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\omega,\omega_{B})=k_{B}\Gamma(\boldsymbol{\kappa}_{1},\boldsymbol{\kappa}_{2},\nu,1)$
(8)

Assuming the reduced Bragg wave vector to be along the x-axis
($\boldsymbol{\kappa}_{B}=(1,0)$) we operate the change of variables:

$(\kappa_{1x},\kappa_{1y})\to(\nu_{1},\nu_{2})$ (9)

Note that we thus replace the two-dimensional wave vector $\kappa_{1}$ by a
set of 2 positive frequencies with the obvious advantage that these appear as
the argument of the delta function in the integrand.

Figure 1: Admissible values of the reduced frequency $\nu_{1}$ seen in the
$(\nu,\nu_{1})$ plane. For $\left|\nu\right|>1$, the integration domain of the
integral (12) corresponds to vertical cuts of the blue strips. For
$\left|\nu\right|<1$, the integration domains of the integrals (20) are given
by vertical cuts of the dark red ($I_{-1,1}$) and light red ($I_{1,-1}$)
regions.

On the upper half-space $\kappa_{1y}\geq 0$ we can derive a simple one-to-one
relationship between old and new integration variables:

$\begin{split}\kappa_{1x}(\nu_{1},\nu_{2})&=\frac{1}{2}\left(1+\nu_{1}^{4}-\nu_{2}^{4}\right)\\\
\kappa_{1y}(\nu_{1},\nu_{2})&=\left|\nu_{1}^{4}-\frac{1}{4}\left(1+\nu_{1}^{4}-\nu_{2}^{4}\right)^{2}\right|^{1/2}\end{split}$
(10)

and the same holds on the lower half-space by changing the sign of
$\kappa_{1y}$ in equation (10). The Jacobian of the transformation is found to
be:

$J(\nu_{1},\nu_{2})=\left|\frac{4\nu_{1}^{3}\nu_{2}^{3}}{\kappa_{1y}(\nu_{1},\nu_{2})}\right|$
(11)

We now evaluate the second-order Doppler spectrum in this new set of variables
by treating separately the frequency intervals between and outside the Bragg
lines.

### Case $\omega>\omega_{B}$ ($n_{1}=n_{2}=1$) or $\omega<-\omega_{B}$
($n_{1}=n_{2}=-1$)

By separating the integrals over the two half-spaces and performing the
integration over the $\nu_{2}$ variable to eliminate the delta function we
obtain:

$\sigma_{2}(\omega)={\cal N\mit}k_{B}^{4}\omega_{B}^{-1}\int_{I(\nu)}{\cal
S\mit}(\nu_{1})\gamma(\nu_{1}){\cal J\mit}(\nu_{1})\ d\nu_{1}$ (12)

where $\gamma$ is the reduced kernel expressed as a function of the sole
reduced frequencies $\nu=\omega/\omega_{B}$ and $\nu_{1}$:

$\gamma(\nu_{1})=\left|\Gamma(\boldsymbol{\kappa}_{1}^{+},\boldsymbol{\kappa}_{2}^{+},\nu,1)\right|^{2}=\left|\Gamma(\boldsymbol{\kappa}_{1}^{-},\boldsymbol{\kappa}_{2}^{-},\nu,1)\right|^{2},$
(13)

$\boldsymbol{\kappa}_{1}^{\pm},\boldsymbol{\kappa}_{2}^{\pm}$ are the reduced
wave vectors in the upper/lower half-space,

$\begin{split}\boldsymbol{\kappa}_{1}^{\pm}&=(\kappa_{1x}(\nu_{1},\left|\nu\right|-\nu_{1}),\pm{\kappa_{1y}(\nu_{1},\left|\nu\right|-\nu_{1})})\\\
\boldsymbol{\kappa}_{2}^{\pm}&=(1-\kappa_{1x}(\nu_{1},\left|\nu\right|-\nu_{1}),\mp{\kappa_{1y}(\nu_{1},\left|\nu\right|-\nu_{1})}),\end{split}$
(14)

${\cal S\mit}$ is the corresponding product of directional ocean spectra,

${\cal
S\mit}(\nu_{1})=S_{d}(n_{1}k_{B}\boldsymbol{\kappa}_{1}^{+})S_{d}(n_{1}k_{B}\boldsymbol{\kappa}_{2}^{+})+S_{d}(n_{1}k_{B}\boldsymbol{\kappa}_{1}^{-})S_{d}(n_{1}k_{B}\boldsymbol{\kappa}_{2}^{-})$
(15)

and ${\cal J\mit}$ is the Jacobian in the reduced frequency variable:

${\cal J\mit}=J(\nu_{1},\left|\nu\right|-\nu_{1})$ (16)

The integration domain $I(\nu)$ is a subset of the positive axis. It is
limited to the values of $\nu_{1}$ for which $\nu_{2}\geq 0$ and for which the
argument of the square root in (10) is positive, that is:

$\nu_{1}^{4}-\frac{1}{4}\left(1+\nu_{1}^{4}-(\left|\nu\right|-\nu_{1})^{4}\right)^{2}\geq
0$ (17)

The limits of the domain can easily be identified analytically with the roots
of the above equation. The integration domain is composed of one or two
disjoint intervals depending on the value of $\nu$:

$I(\nu)=\left[\frac{\nu^{2}-1}{2\left|\nu\right|},\frac{\left|\nu\right|-\sqrt{2-\nu^{2}}}{2}\right]\cup\left[\frac{\left|\nu\right|+\sqrt{2-\nu^{2}}}{2},\frac{1+\nu^{2}}{2\left|\nu\right|}\right]$
(18)

if $1\leq\left|\nu\right|\leq\sqrt{2}$ and

$I(\nu)=\left[\frac{\nu^{2}-1}{2\left|\nu\right|},\frac{\nu^{2}+1}{2\left|\nu\right|}\right]$
(19)

if $\left|\nu\right|>\sqrt{2}$. Note that the integration domain is splitted
in two disjoints intervals for $1<\left|\nu\right|<\sqrt{2}$. However, these
two intervals are interchanged by the transformation
$\nu_{1}\mapsto\nu_{1}-\nu$. Since this transformation amounts to swap
$\nu_{1}$ and $\nu_{2}$ in the integrand, this leaves $\gamma$, ${\cal J\mit}$
and ${\cal S\mit}$ unchanged. It follows that the 2 integrals on the lower and
upper interval are identical so that it suffices to calculate one of them and
multiply the result by two.

### Case $\left|\omega\right|<\omega_{B}$ ($n_{1}=-n_{2}=\pm 1$)

Similar calculations lead to the following expression for the second-order
Doppler spectrum inside the Bragg lines:

$\sigma_{2}(\omega)={\cal
N\mit}k_{B}^{4}\omega_{B}^{-1}\left\\{\int_{I_{1,-1}(\nu)}+\int_{I_{-1,1}(\nu)}\right\\}{\cal
S\mit}(\nu_{1})\gamma(\nu_{1}){\cal J\mit}(\nu_{1})\ d\nu_{1}$ (20)

with the following adaptation of the integrand. On the interval
$I_{1,-1}(\nu)$, the variable $\left|\nu\right|-\nu_{1}$ should be replaced by
$\nu_{1}-\nu$ in equations (14) and (16), and equation (15) should be replaced
by:

${\cal
S\mit}(\nu_{1})=S_{d}(k_{B}\boldsymbol{\kappa}_{1}^{+})S_{d}(-k_{B}\boldsymbol{\kappa}_{2}^{+})+S_{d}(k_{B}\boldsymbol{\kappa}_{1}^{-})S_{d}(-k_{B}\boldsymbol{\kappa}_{2}^{-})$
(21)

On the interval $I_{-1,1}(\nu)$, the variable $\left|\nu\right|-\nu_{1}$
should be replaced by $\nu_{1}+\nu$ in equations (14) and (16), and equation
(15) should be replaced by:

${\cal
S\mit}(\nu_{1})=S_{d}(-k_{B}\boldsymbol{\kappa}_{1}^{+})S_{d}(k_{B}\boldsymbol{\kappa}_{2}^{+})+S_{d}(-k_{B}\boldsymbol{\kappa}_{1}^{-})S_{d}(k_{B}\boldsymbol{\kappa}_{2}^{-})$
(22)

The integration domain $I_{1,-1}(\nu)$ is constrained by the conditions:

$\nu_{1}^{4}-\frac{1}{4}\left(1+\nu_{1}^{4}-(\nu-\nu_{1})^{4}\right)^{2}\geq
0,\ \nu_{1}\geq\nu,$ (23)

leading to:

$I_{1,-1}(\nu)=\left[\frac{\nu+\sqrt{2-\nu^{2}}}{2},\frac{1+\mbox{sign}(\nu)\nu^{2}}{2\left|\nu\right|}\right]$
(24)

Likewise the integration domain $I_{-1,1}(\nu)$ is constrained by the
conditions:

$\nu_{1}^{4}-\frac{1}{4}\left(1+\nu_{1}^{4}-(\nu+\nu_{1})^{4}\right)^{2}\geq
0,\ \nu_{1}\geq-\nu$ (25)

leading to:

$I_{-1,1}(\nu)=\left[\frac{-\nu+\sqrt{2-\nu^{2}}}{2},\frac{1-\mbox{sign}(\nu)\nu^{2}}{2\left|\nu\right|}\right]$
(26)

The complete integration domain defined by the union of $I(\nu)$ (for
$\left|\nu\right|\geq 1$), $I_{1,-1}(\nu)$ and $I_{-1,1}(\nu)$ (for
$\left|\nu\right|<1$) is shown in the $(\nu,\nu_{1})$ plane in Figure 1. Note
that the transformation $\nu_{1}\mapsto\nu_{1}-\nu$ maps the interval
$I_{1,-1}(\nu)$ to $I_{-1,1}(\nu)$. Again, since this transformation amounts
to swap $\nu_{1}$ and $\nu_{2}$ in the integrand, this leaves $\gamma$ and
${\cal J\mit}$ unchanged in the integral and transforms eq. (21) to eq. (22).
It follows that the 2 integrals in eq. (20) are identical so that it suffices
to calculate one of them:

$\sigma_{2}(\omega)=2{\cal
N\mit}k_{B}^{4}\omega_{B}^{-1}\int_{I_{-1,1}(\nu)}{\cal
S\mit}(\nu_{1})\gamma(\nu_{1}){\cal J\mit}(\nu_{1})\ d\nu_{1}$ (27)

When $\nu_{1}$ approaches the border of the domain, corresponding to one end
of the intervals $I(\nu),I_{1,-1}(\nu)$ or $I_{-1,1}(\nu)$, the Jacobian term
${\cal J\mit}[\nu_{1})$ has an integrable $1/\sqrt{\epsilon}$ singularity,
except at the particular values $\nu=\pm\sqrt{2}$, where its denominator has a
double root causing an non-integrable $1/\epsilon$ singularity. This recovers
the result established in [12] that the Doppler spectrum has a logarithmic
singularity at $\nu=\pm\sqrt{2}$. Apart from these 2 specific points, the
singularity of the Jacobian at the borders of the domain can be handled
properly by using the Gauss-Jabobi quadrature rule for numerical integration.
We recall that:

$\int_{-1}^{+1}(x-1)^{-1/2}(x+1)^{-1/2}f(x)dx\simeq\sum_{1}^{N}w_{j}f(x_{j}),$
(28)

where the nodes of the quadrature $x_{j}$ are the zeros of the Jacobi
polynomials and $w_{j}$ the associated weights, the approximation being exact
whenever $f$ is constant. Many open source numerical routines can be found on
the internet. This allows a fast and accurate calculation of the second-order
ocean Doppler spectrum (12) at the cost of an ordinary one-dimensional
integral.

The sea state dependence in the expression of the second-order Doppler
spectrum enters through the ocean wave spectrum. Following the standard
representation we write the directional wave spectrum in the form:

$S_{d}(\boldsymbol{k})=k^{-1}S_{o}(k)D(k,\theta)$ (29)

where $S_{o}(k)$ is the omnidirectional spectrum and $D(\theta)$ is an angular
spreading function with the normalization $\int_{0}^{2\pi}D(\theta)d\theta=1$
and the angle $\theta$ is measured with respect to the radar beam. It can also
be expressed in term of the classical frequency spectrum $S_{f}(\omega)$
(using $S_{o}(k)=S_{f}(\omega)d\omega/dk$ and $\omega=\sqrt{gk}$):

$S_{d}(\boldsymbol{k})=\frac{1}{2}g^{2}\omega^{-3}S_{f}(\omega)D(\omega,\theta),$
(30)

For the numerical simulations we have employed the classical Pierson-Moskowitz
(PM) frequency spectrum:

$S_{f}(\omega)=Ag^{2}\omega^{-5}\exp(-B(g/U_{10}\omega)^{4})$ (31)

corresponding to an omnidirectional wavenumber spectrum

$S_{o}(k)=\frac{A}{2}k^{-3}\exp(-Bg^{2}/(U_{10}^{4}k^{2})),$ (32)

where $A=0.0081$, $B=0.74$, $U_{10}$ is the wind speed at 10 meter above the
sea surface. For the angular spreading function we chose the classical
cardioïd shaped spreading function:

$D(\theta)=\alpha(\epsilon+(1-\epsilon)\cos^{4}((\theta-\theta_{w})/2)),$ (33)

where $\alpha$ is the normalization constant, $\theta_{w}$ is the direction of
wind with respect to the radar beam ($\theta_{w}=0$ upwind) and
$\epsilon=0.05$ is a small parameter that allows to have a fraction of waves
propagating against the wind, as it is often observed. An example of the
resulting computation for the upwind ocean Doppler spectrum at wind speed
$U_{10}=10$ m/s is shown in Figure 2 at the commonly used radar frequency
$f_{0}=16$ MHz. We superimposed the smoothed version of this spectrum, which
removes the delta-like singularity at $\omega=\pm\omega_{B}$ and the
logarithmic singularity at $\omega=\pm\sqrt{2}\omega_{B}$ and corresponds to a
practical situation. The computation was assessed by comparison (not shown
here) with the aforementioned direct method based on a two-dimensional
integration in the $\boldsymbol{k}_{1}$ space using a finite filter $\Phi$
instead of a Dirac function.

Figure 2: Second-order Doppler spectrum at radar frequency $f_{0}=16$ MHz
calculated with the improved formulation (12) as a function of the normalized
Doppler frequency $\nu=\omega/\omega_{B}$. The ocean wave spectrum is a
directional PM spectrum at 10 m/s wind speed with the radar looking upwind
(black solid lines). Also shown is a smoothed Doppler spectrum corresponding
to a finite frequency resolution of $0.01$ Hz (red solid lines). The black
dashed vertical lines mark the $\pm 2^{1/2}$ and $\pm 2^{3/4}$ singularities.

##  III The Zeta function

We now consider the following function, henceforth referred to as the “Zeta
function”, obtained by removing the coupling coefficient from the expression
of the second-order ocean Doppler spectrum and replacing it with $k_{B}^{2}$
to keep the same dimension:

$\begin{split}\zeta(\omega)&={\cal N\mit}k_{B}^{2}\sum_{n_{1}=\pm
1}\sum_{n_{2}=\pm 1}\int
S_{d}(n_{1}\boldsymbol{k}_{1})S_{d}(n_{2}\boldsymbol{k}_{2})\\\
&\times\delta(\omega-
n_{1}\omega_{1}-n_{2}\omega_{2})d\boldsymbol{k}_{1}\end{split}$ (34)

As first suggested in [9], the main parameters of the ocean wave spectrum can
be estimated with the first two moments of the Zeta function and we will
follow the same argumentation. We first integrate the Zeta function over the
frequency variable to get rid of the delta function in the integrand of
$\zeta$. The second step of the demonstration is obtained by noting that the
ocean wave spectrum has a dominant value near the peak wavenumber
$k_{p}<<k_{B}$ so that its leading contribution can be extracted from the
integral. This occurs when one of the variable $\boldsymbol{k}_{1}$ or
$\boldsymbol{k}_{2}$ in the integrand is close to zero, leading in the end to:

$\int_{-\infty}^{+\infty}\zeta(\omega)d\omega\simeq 4{\cal
N\mit}k_{B}^{2}\sum_{n_{1}=\pm 1}S_{d}(n_{1}\boldsymbol{k}_{B})\int
S_{d}(\boldsymbol{k})d\boldsymbol{k}$ (35)

The integral of the right-hand side is the variance of wave elevations:

$H_{0}^{2}=\int S_{d}(\boldsymbol{k})d\boldsymbol{k}$ (36)

while the discrete summation is the integral of the first-order Bragg peaks.
Hence the Significant Wave Height (SWH) $H_{s}=4H_{0}$ can be recovered from
the Zeta function:

$H_{s}^{2}=\frac{4\int_{-\infty}^{+\infty}\zeta(\omega)d\omega}{k_{B}^{2}\int_{-\infty}^{+\infty}\sigma_{1}(\omega)d\omega}$
(37)

To estimate the mean wave frequency we consider the positive frequency side of
the Zeta function spectrum (assuming it is stronger than the negative
frequency side) and take the following ratio:

${\Omega}=\frac{\int_{\omega_{B}}^{+\infty}(\omega-\omega_{B})\zeta(\omega)d\omega}{\int_{\omega_{B}}^{+\infty}\zeta(\omega)d\omega}$
(38)

After integration of the delta function this can be rewritten as:

${\Omega}=\frac{\int(\omega_{1}+\omega_{2}-\omega_{B})S_{d}(\boldsymbol{k}_{1})S_{d}(\boldsymbol{k}_{2})d\boldsymbol{k}_{1}}{\int
S_{d}(\boldsymbol{k}_{1})S_{d}(\boldsymbol{k}_{2})d\boldsymbol{k}_{1}}$ (39)

Again, the leading contribution of the integrand occurs when one of the
frequency $k_{1}$ or $k_{2}$ is close to zero and the other close to $k_{B}$,
so that one of the $\omega_{1},\omega_{2}$ variable cancel out with the
$-\omega_{B}$ term. This identifies $\Omega$ as the centroïd frequency:

${\Omega}\simeq\frac{\int\omega(\boldsymbol{k})S_{d}(\boldsymbol{k})d\boldsymbol{k}}{\int
S_{d}(\boldsymbol{k})d\boldsymbol{k}}=\frac{2\pi}{T}$ (40)

which defines the mean period $T$.

Figure 3: Estimation of the SWH (left panel) and the mean period (right panel)
from the Zeta function as a function of wind speed for a directional PM
spectrum with the radar aiming in the upwind (solid lines) or crosswind
(dashed lines) direction. The black squares show the theoretical value of the
SWH and mean period. The various radar frequencies correspond to different
color codes: 5 MHz (green); 10 MHz (magenta); 15 MHz (cyan); 20 MHz (red) ; 25
MHz (red).

We performed numerical simulations, shown in Figure 3, to estimate the SWH and
the mean wave period $T$ from the Zeta function for the directional PM
spectrum at various wind speed and direction at several radar frequencies
covering the typical range of HF radars (from 5 to 25 MHz). We found that the
SWH is increasingly well estimated by formula (36) as the radar frequency
varies from 10 to 25 MHz but is poorly estimated at the lowest frequency (5
MHz). The same conclusion applies to the mean period, which is overestimated
by formula (40-38) but becomes more accurate at higher radar frequencies and
larger wind speed. This is consistent with the expectation that the
approximation that was used to evaluate the integral of the Zeta function (eq.
(35)) necessitates a high ratio of the Bragg wavenumber $k_{B}$ to the peak
wavenumber $k_{p}$. The observed dependence of the estimated wave parameters
on the radar frequency in Figure (3) suggests that the estimated SWH
($\hat{H}_{s}$) should be corrected for by a multiplicative bias while the
estimated mean period ($\hat{T}$) should be corrected for by an additive bias:

$\hat{H}_{s}\to\alpha\hat{H}_{s},\ \hat{T}\to\hat{T}-T_{0}$ (41)

By adjusting the average of the upwind and cross-wind estimation to the
theoretical wave parameters for the ensemble of wind speeds we obtained
parameters $\alpha$ and $T_{0}$ that depend on the sole radar frequency. These
values are summarized in Table I. The corrective factor for the SWH ranges
from $0.93$ to $0.97$ when increasing the radar frequency to $10$ to $25$ MHz
(we discard the case $5$ MHz) while the bias to subtract from the mean period
ranges from $1.25$ to $0.4$ second.

Radar frequency (MHz) | 10 | 15 | 20 | 25
---|---|---|---|---
$\alpha$ | 0.93 | 0.95 | 0.96 | 0.97
$T_{0}$ (sec) | 1.25 | 0.76 | 0.53 | 0.4
TABLE I: Corrective factors to apply in the derivation of the sea state
parameters from the Zeta function

##  IV The weighting function

Figure 4: New weighting function (blue solid line) as compared to Barrick’s
weighting function (red solid line). The simple analytical fit (47) is also
shown (cyan line). The dashed vertical lines mark the 2 approximate frequency
intervals ($[0.35-0.8]$ and $[1.2-1.7]$ where the weighting function is
useful.

In order to make the expression of the ocean Doppler spectrum analytically
tractable for sea state inversion, Barrick [8] proposed to replace the
coupling coefficient $\left|\Gamma\right|^{2}$ by the so-called weighting
function, obtained by performing an angular average
$\overline{\left|\Gamma\right|^{2}}$ over the admissible wavenumbers:

$\frac{k_{0}^{2}}{8}{\cal W\mit_{B}}(\nu)=\overline{\left|\Gamma\right|^{2}}$
(42)

We found this weighting function difficult to implement from the original
publication and we used the numerical routine provided by the free RadarWIC
library [13], that uses a digitized version of the Figure 3 of [8]. Our
reformulation of the second-order integral (12) suggests the construction of a
different weighting function, using the same normalization factor as Barrick:

$\frac{k_{0}^{2}}{8}{\cal W\mit}(\nu)=k_{B}^{2}\langle\gamma\rangle_{I(\nu)}$
(43)

where $\langle{\gamma}\rangle_{I(\nu)}$ is an average value of $\gamma$ over
the reduced frequency domain $I_{\nu}$. For $\nu>1$, this domain is narrow and
it is justified to defined it with the mean value of the coupling coefficient
over this interval:

$\langle\gamma\rangle_{I(\nu)}=\frac{\int_{I(\nu)}\gamma(\nu_{1})d\nu_{1}}{\int_{I(\nu)}d\nu_{1}},\
\nu>1,$ (44)

For $0<\nu<1$, the values of both the kernel $\gamma$ and the product of
spectra in the integrands (20) are more contrasted, so that is not relevant to
retain the mean value. However, we observed numerically and checked
analytically that the integrand reach their dominant value at the edges of the
integration domain. To give more weight to the latter, we replaced the mean by
a simple 3 points average between the ends and the middle of the lower
interval $I_{-1,1}(\nu)$ (we recall that the integral on the upper interval
$I_{1,-1}(\nu)$ is identical):

$\langle\gamma\rangle_{I(\nu)}=\frac{1}{3}\left(\gamma(A_{\nu})+\gamma\left(\frac{A_{\nu}+B_{\nu}}{2}\right)+\gamma(B_{\nu})\right)$
(45)

with

$A(\nu)=\frac{\nu^{2}-1}{2\nu},\ B(\nu)=\frac{\nu-\sqrt{2-\nu^{2}}}{2}$ (46)

For $\nu<0$, the weighting function is extended by parity (i.e., ${\cal
W\mit}(-\nu)={\cal W\mit}(\nu)$). Barrick’s weighting function and ours are
shown in Figure 4. Both exhibit a singularity at $\nu=2^{3/4}$ and follow the
same trends but are quite different in magnitude. The difference in the
definitions of the two weighting functions lies in the domain and variables of
integration. Note that the new weighting function is no longer singular at
$\nu=2^{1/2}$ unlike Barrick’s weighting function. This is because we do not
include the Jacobian, which is at the origin of this singularity, in the
definition. The new weighting function can be easily implemented from (44) and
(LABEL:meangamma2) but it is more convenient to have an analytical expression
in the form of an empirical fit. We found the following elementary linear fit,
which is enough for our purposes:

${\cal W\mit}(\nu)=\left\\{\begin{split}\exp(13.87\nu^{2}-18.38\nu+7.72)&\
\mathrm{if}\ 0<\nu<0.63\\\ 4.64&\ \mathrm{if}\ 0.63<\nu<1\\\ -2.33\nu+5&\
\mathrm{if}\ 1<\nu<1.45\\\ 34.87\nu-48.93&\ \mathrm{if}\ \nu>1.45\\\
\end{split}\right.$ (47)

To simplify the fit we did not retain the peak of the $2^{3/4}$ singularity,
which has a negligible role in the integration of the Zeta function. Note also
that the new weighting function has a discontinuity at $\nu=1$, due to the
adaptation of the definition of the mean kernel for $0<\nu<1$. However, this
discontinuity has no consequence since the weighting function is only useful
in two small intervals of frequencies lying on each side of the Bragg
frequency, which are approximately $[0.35-0.8]$ and $[1.2-1.7]$. They
correspond to the intervals where the Zeta function has non-negligible values
(see Figure 5).

Figure 5: Comparison of the Zeta function (black lines) and its evaluation
with the ocean Doppler spectrum (eq. 48) using Barrick’s weighting function
(red line) and ours with the simple fit (47) (blue lines). The ocean Doppler
spectrum was simulated in the upwind (left column) and crosswind (right
column) direction using a directional PM spectrum by 7 m/s (top panel), 10 m/s
(middle panel) and 15 m/s (bottom panel) wind speed at a radar frequency $16$
MHz. The three functions have been smoothed with a window of width 0.01 Hz and
are shown in linear scale (i.e., not dB) to highlight the frequency intervals
where it the are non negligible.

Approximating the coupling coefficient with its average value in the
expression (12) of the second-order Doppler spectrum leads to the following
relation with the Zeta function:

$\zeta(\omega)\simeq 2^{5}\frac{\sigma_{2}(\omega)}{{\cal
W\mit}(\omega/\omega_{B})}$ (48)

The properties of the Zeta function therefore transfer to the ocean Doppler
spectrum provided the weighting function approximation is sufficiently
accurate. In particular this leads to the following formulas for the
estimation of the SWH and mean period:

$(k_{0}H_{0})^{2}=2\alpha^{2}\int_{-\infty}^{+\infty}R_{{\cal
W\mit}}(\omega)d\omega$ (49)

and

$T=\frac{2\pi\int_{\omega_{B}}^{+\infty}R_{{\cal
W\mit}}(\omega)d\omega}{\int_{\omega_{B}}^{+\infty}(\omega-\omega_{B})R_{{\cal
W\mit}}(\omega)d\omega}-T_{0}$ (50)

where we define the ratio of the weighted second-order to first order Doppler
spectrum following the notation of [14] as:

$R_{{\cal W\mit}}(\omega)=\frac{\sigma_{2}(\omega)/{\cal
W\mit}(\omega/\omega_{B})}{\int_{-\infty}^{+\infty}\sigma_{1}(\omega)d\omega}$
(51)

Apart from the corrective factors $\alpha$ and $T_{0}$ these expressions were
already established by Barrick (eq. (11) and (13) of [9]). Their accuracy is
related, on one hand on the ability of the weighting function to recover the
Zeta function (48), and on the other hand on the quality of the estimation
(36), (38) and (41) based on the Zeta function. To elucidate these points, we
performed numerical simulations at various wind speed and radar frequencies
using the directional PM spectrum. Figure 5 shows a comparison between the
exact expression of the Zeta function and its approximation (48) at the radar
frequency 16 MHz. We tested both Barrick’s weighting function and the fitted
version of the new weighting function (47). From a visual inspection (which is
confirmed numerically) it is clear than the latter is in general more
accurate, especially in the upwind direction. As a result, the estimation of
the sea states parameters from the second-order Doppler spectrum (eqs.
(49)-(50)) is improved when using this new weighting function, even in its
elementary fitted version (47). Figure 6 shows the percentage of error between
the estimated and actual values of the SWH and mean period as a function of
dimensionless parameter $k_{0}H_{s}$ for a directional PM wave spectrum at
various wind speeds and radar frequencies. For reference, the result of the
bias-corrected estimation with the Zeta function (eqs (36), (38) and (41)) is
also shown. An overall error reduction is obtained with the new weighting
function for $k_{0}H_{s}>0.5$. However, as already noted by other authors
(e.g. [15]), inconsistent estimations are obtained in the upwind and crosswind
direction; their discrepancy is reduced when using the new weighting function
but remains, at best, of the order of $9\%$ of the SWH for $k_{0}H_{s}>0.5$
and can reach $25\%$ for $k_{0}H_{s}<0.5$ while it is of the order of $10\%$
for the mean period. As seen, there is also a small discrepancy between the
upwind and crosswind estimation from the Zeta function, especially for lower
values of $k_{O}H_{s}$ and for the mean period.

Note that an empirical correction (49) for the SWH arising from Barrick’s
formula has already been employed by several authors to improve its
consistency with High-Frequency radar measurements [16, 17, 18]. At a radar
frequency of 25.4 MHZ, [18] found the best fit with $\alpha\simeq 0.58$ and
[16] with $\alpha\simeq 0.55$. However, our value $\alpha=0.97$ at 25 MHz has
been obtained from the Zeta function and not from its approximation (48) with
the weighting function, which can explain the difference.

Figure 6: Relative error (in percent) in the estimation of the SWH (top panel)
and the mean period (bottom panel) from the second-order Doppler spectrum
(eqs. (49)-(50)) using either Barrick’s weighting function (red symbols) or
fitted weighting function (47) (blue symbols). The errors are shown as a
function of the dimensionless roughness parameter $k_{0}H_{s}$, wherein the
results of the four frequency bands 10,15,20 and 25 MHz have been
concatenated. The estimation based on the Zeta function is shown for reference
in black symbols. The circles correspond to the upwind direction and the
crosses to the crosswind direction.

##  V A direct relation to the wave spectrum

In a series of papers [19, 20, 14, 21], attempts were made to find a simple
linear relationship between the ocean wave spectrum and the second-order
Doppler spectrum shift by the Bragg frequency, as originally proposed by [3].
An empirical relation of the type:

$k_{0}^{2}S_{f}(\omega)=2A(\omega)R_{{\cal W\mit}}(\omega\pm\omega_{B})$ (52)

which can be seen as the spectral form of Barrick’s formula (49), has been
sought with a frequency-dependent calibration coefficient $A(\omega)$ and
different weighting functions (taken to be one for [19, 20] and adapted to
wind-waves and swell for [14, 21]). Even though this method has shown good
results when employed with experimental data, it is not fully satisfactory as
it requires manual calibration for each radar site and bearing and is limited
to the range of frequencies where the second-order Doppler spectrum can be
linearized, i.e. not too far from the Bragg frequency. We propose here a
nonlinear approximation, inspired by the reformulation (12) of the second-
order Doppler spectrum as a frequency integral over the variable domain
$I(\nu)$. For $\left|\nu\right|>1$, the integration domain $I(\nu)$ is narrow
(see Figure 1) and we can therefore approximate the above integral by
extracting the mean value $\overline{{\cal S\mit}}$ of the function ${\cal
S\mit}$ over this interval, leading to:

$\sigma_{2}(\omega)\simeq{\cal N\mit}k_{B}^{4}\omega_{B}^{-1}\overline{{\cal
S\mit}}(\nu){\cal F\mit}(\nu),$ (53)

where:

${\cal F\mit}(\nu)=\int_{I(\nu)}\gamma(\nu_{1}){\cal J\mit}(\nu_{1})\
d\nu_{1}$ (54)

is a deterministic function that does not depend on sea state nor on the radar
frequency and can thus be calculated once for all. This function, which is
shown in Figure 7, carries the 2 singularities at $\nu=2^{1/2}$ and
$\nu=2^{3/4}$ but has otherwise a smooth polynomial behavior. For $\nu>1.7$ it
can be described by the following fit with an excellent accuracy:

${\cal F\mit}(\nu)=0.0592\nu^{3}-0.2935\nu^{2}+0.5038\nu-0.2958$ (55)

Figure 7: Function ${\cal F\mit}(\nu)$ (blue solid line) as defined in (54)
and expressed in decibels, together with its fit (55) for $\nu>1.7$ (dashed
red line).

To obtain a tractable expression we define $\overline{{\cal S\mit}}$ in the
following way. For $1\leq\nu\leq\sqrt{2}$, the integration domain $I(\nu)$
consists of 2 narrow intervals given by (18) and we may approximate ${\cal
S\mit}$ by its average between the endpoints of each interval; for
$\nu>\sqrt{2}$, the integration domain $I(\nu)$ consists of a single wider
interval on which we may approximate ${\cal S\mit}$ by its average between the
midpoint and the endpoints. Altogether, this leads to the following
definition:

$\overline{{\cal S\mit}}(\nu)=\left\\{\begin{array}[]{l}\frac{1}{4}({\cal
S\mit}_{1}^{+}+{\cal S\mit}_{1}^{-}+{\cal S\mit}_{2}^{+}+{\cal
S\mit}_{2}^{-}),\ \mathrm{if}\ 1\leq\nu\leq\sqrt{2}\\\ \\\ \frac{1}{4}(2{\cal
S\mit}_{0}+{\cal S\mit}_{2}^{+}+{\cal S\mit}_{2}^{-}),\ \mathrm{if}\
\nu>\sqrt{2}\end{array}\right.$ (56)

with

${\cal S\mit}_{0}={\cal S\mit}(\nu/2),\ {\cal S\mit}_{1}^{\pm}={\cal
S\mit}\left(\frac{\nu\pm\sqrt{2-\nu^{2}}}{2}\right),\ {\cal
S\mit}_{2}^{\pm}={\cal S\mit}\left(\frac{\nu^{2}\pm 1}{2\nu}\right)$ (57)

From the relation (30) we can convert this last expression in terms of
directional frequency spectrum:

$\begin{split}{\cal
S\mit}_{0}&=\frac{32g^{4}}{\omega^{6}}S^{2}_{f}(\omega/2)(D(\omega/2,\theta_{0})D(\omega/2,-\theta_{0})),\\\
{\cal
S\mit}_{1}^{\pm}&=\frac{32g^{4}}{(\omega_{1}^{+}\omega_{1}^{-})^{3}}S_{f}(\omega^{-}_{1})S_{f}(\omega^{+}_{1})D(\omega^{+}_{1},0)D(\omega^{-}_{1},0)\\\
{\cal
S\mit}_{2}^{\pm}&=\frac{32g^{4}}{(\omega_{2}^{+}\omega_{2}^{-})^{3}}S_{f}(\omega^{-}_{2})S_{f}(\omega^{+}_{2})D(\omega^{+}_{2},0)D(\omega^{-}_{2},0)\\\
\omega_{1}^{\pm}&=\frac{\omega\pm\sqrt{2\omega_{B}^{2}-\omega^{2}}}{2},\
\omega^{\pm}_{2}=\frac{\omega^{2}\pm\omega^{2}_{B}}{2\omega}\\\
\theta_{0}&=\mathrm{atan}\left(\frac{1}{2}\sqrt{\frac{\omega^{4}}{\omega_{B}^{4}}-4}\right)\end{split}$
(58)

Similar relations hold for negative frequencies ($\nu<-1$), with the argument
of the spreading function shifted by $\pi$ ($D(\omega/2,\pm\theta_{0})\to
D(\left|\omega\right|/2,\pm\theta_{0}+\pi)$, $D(\omega_{1,2}^{\pm},0)\to
D(\left|\omega_{1,2}^{\pm}\right|,\pi)$). Figure 8 shows a comparison of the
exact formula (12) and its approximated version using equations (53) for a PM
spectrum by 10 m/s wind speed for an upwind looking radar at 16 MHz. We found
the approximation good for $\nu>1.5$ (less than 2 dB error) and excellent for
$\nu>1.75$ (less than 0.5 dB error). We checked numerically other wind speeds
and radar frequencies with the same conclusion. This formula shows that
starting from $\omega=2^{1/2}\omega_{B}$, the part of the ocean wave frequency
spectrum that contributes to the Doppler spectrum at a given Doppler frequency
$\omega$ is a narrow interval around $\omega/2$. The directions of waves
involved in this process are those having an angle smaller than $\theta_{0}$
with respect to the radar beam. This last angle is null at
$\omega=2^{1/2}\omega_{B}$, reaches $45$ degree at $2^{3/4}\omega_{B}$, and
$60$ degree at $2\omega_{B}$ and converges to the cross-beam directions $90$
degree as the frequency keeps increasing. Note that in the vicinity of the
first singularity $2^{1/2}\omega_{B}$, the angle $\theta_{0}$ is close to zero
so that the radar dominantly “sees” the waves that are propagating along or
against the radar look direction.

##  VI Conclusion

We have proposed a formulation of the second-order ocean Doppler spectrum
based on an integration of the sole frequency variable. This new expression
allows for efficient and accurate numerical evaluation of the Doppler
spectrum. In addition, it leads naturally to an alternative weighting function
that can be used to improve the estimation of the main sea state parameters
after the classical method proposed by Barrick in 1977 ([9, 8]) and reduces
the discrepancy between the upwind and crosswind estimation. Nevertheless,
there are some limitations to this method, which are not related to the
quality of the weighting function but to the structure of the Zeta function. A
bias correction for both the SWH and mean period is necessary, with values
depending on the radar frequency, and there remains a small but non negligible
discrepancy between the upwind and crosswind estimation from a single radar
site. A second important consequence of our approach is the derivation of an
analytical nonlinear approximation for the ocean wave spectrum for Doppler
frequencies well beyond the Bragg frequency. This formula might be a useful
complement to classical techniques of ocean wave spectrum inversion which are
often limited to the frequencies close to the Bragg frequency (typically in
the range $0.4<\nu<1.6$) where the problem can be linearized (see e.g. the
discussions in [4] or [15]).

Figure 8: Exact numerical calculation (black lines) of the second-order
Doppler spectrum compared with its analytical approximation with equations
(53) (red lines) for a directional PM spectrum. The radar is looking upwind
(solid lines) or crosswind (dashed lines) at 16 MHz and the wind speed ranges
from 7 m/s to 15 m/s from left to right.

Acknowledgments: This research has been funded by the Agence Nationale de la
Recherche (ANR) under grant ANR-22-ASTR-0006-01 (ROSMED project: “Radar à
Ondes de Surface en MEDiterranée”)

## References

  * [1] Pablo Lorente et al.,  “Coastal high-frequency radars in the Mediterranean - Part 1: Status of operations and a framework for future development,”  Ocean Science, vol. 18, pp. 761–795, 2022.
  * [2] D E Barrick,  “Remote sensing of sea state by radar,”  in Remote sensing of the Troposphere. 1972, vol. 12, pp. 1–46, VE Derr, Editor, US Government,  Publisher: Washington, DC: US Government Printing Office.
  * [3] Klaus Hasselmann,  “Determination of ocean wave spectra from doppler radio return from the sea surface,”  Nature Physical Science, vol. 229, no. 1, pp. 16–17, 1971.
  * [4] Lucy R Wyatt,  “Limits to the inversion of HF radar backscatter for ocean wave measurement,”  Journal of Atmospheric and Oceanic Technology, vol. 17, no. 12, pp. 1651–1666, 2000.
  * [5] Belinda Lipa,  “Derivation of directional ocean-wave spectra by integral inversion of second-order radar echoes,”  Radio Science, vol. 12, no. 3, pp. 425–434, 1977.
  * [6] L Wyatt,  “A relaxation method for integral inversion applied to HF radar measurement of the ocean wave directional spectrum,”  Int J Remote Sens, vol. 11, no. 8, pp. 1481–1494, 1990.
  * [7] Yukiharu Hisaki,  “Nonlinear inversion of the integral equation to estimate ocean wave spectra from HF radar,”  Radio science, vol. 31, no. 1, pp. 25–39, 1996.
  * [8] Donald E Barrick,  “The ocean waveheight nondirectional spectrum from inversion of the HF sea-echo Doppler spectrum,”  Remote Sensing of Environment, vol. 6, no. 3, pp. 201–227, 1977,  Publisher: Elsevier.
  * [9] Donald E Barrick,  “Extraction of wave parameters from measured HF radar sea-echo Doppler spectra,”  Radio Science, vol. 12, no. 3, pp. 415–424, 1977,  Publisher: AGU.
  * [10] B J Lipa and D E Barrick,  “Extraction of sea state from HF radar sea echo: Mathematical theory and modeling,”  Radio Science, vol. 21, no. 1, pp. 81–100, 1986,  Publisher: Wiley Online Library.
  * [11] J Walsh and EW Gill,  “An analysis of the scattering of high-frequency electromagnetic radiation from rough surfaces with application to pulse radar operating in backscatter mode,”  Radio Science, vol. 35, no. 6, pp. 1337–1359, 2000.
  * [12] Dmitry V Ivonin, Victor I Shrira, and Pierre Broche,  “On the singular nature of the second-order peaks in HF radar sea echo,”  IEEE Journal of Oceanic Engineering, vol. 31, no. 4, pp. 751–767, 2006.
  * [13] Douglas Cahl, George Voulgaris, and Zaid Alattabi,  “Wave radar inversion code (waveric) v1.0.1 (version v1.0.1),” http://doi.org/10.5281/zenodo.2643696, April 2019.
  * [14] Zaid R Alattabi, Douglas Cahl, and George Voulgaris,  “Swell and wind wave inversion using a single very high frequency (VHF) radar,”  Journal of Atmospheric and Oceanic Technology, vol. 36, no. 6, pp. 987–1013, 2019.
  * [15] L. Wyatt,  Ocean Remote Sensing Technologies—High-Frequency, Marine and GNSS-Based Radar, edited by Huang, Weimin and Gill, Eric W, chapter Ocean wave measurement, pp. 145–178,  SciTech Publishing Raleigh, NC, USA, 2021.
  * [16] SF Heron and ML Heron,  “A comparison of algorithms for extracting significant wave height from HF radar ocean backscatter spectra,”  Journal of Atmospheric and oceanic technology, vol. 15, no. 5, pp. 1157–1163, 1998.
  * [17] Joseph W Maresca Jr and TM Georges,  “Measuring rms wave height and the scalar ocean wave spectrum with HF skywave radar,”  Journal of Geophysical Research: Oceans, vol. 85, no. C5, pp. 2759–2771, 1980.
  * [18] Rafael J Ramos, Hans C Graber, and Brian K Haus,  “Observation of wave energy evolution in coastal areas using HF radar,”  Journal of atmospheric and oceanic technology, vol. 26, no. 9, pp. 1891–1909, 2009.
  * [19] Klaus-Werner Gurgel, Heinz-Hermann Essen, and Thomas Schlick,  “An empirical method to derive ocean waves from second-order Bragg scattering: Prospects and limitations,”  IEEE Journal of Oceanic Engineering, vol. 31, no. 4, pp. 804–811, 2006,  Publisher: IEEE.
  * [20] Guiomar Lopez, Daniel C Conley, and Deborah Greaves,  “Calibration, validation, and analysis of an empirical algorithm for the retrieval of wave spectra from HF radar sea echo,”  Journal of Atmospheric and Oceanic Technology, vol. 33, no. 2, pp. 245–261, 2016.
  * [21] Zaid R Al-Attabi, George Voulgaris, and Daniel C Conley,  “Evaluation and validation of HF radar swell and wind-wave inversion method,”  Journal of Atmospheric and Oceanic Technology, vol. 38, no. 10, pp. 1747–1775, 2021.