In ocean acoustics, shallow water propagation is conveniently described using normal mode propagation. This article proposes a framework to describe the polarization of normal modes, as measured using a particle velocity sensor in the water column. To do so, the article introduces the Stokes parameters, a set of four real-valued quantities widely used to describe polarization properties in wave physics, notably for light. Stokes parameters of acoustic normal modes are theoretically derived, and a signal processing framework to estimate them is introduced. The concept of the polarization spectrogram, which enables the visualization of the Stokes parameters using data from a single vector sensor, is also introduced. The whole framework is illustrated on simulated data as well as on experimental data collected during the 2017 Seabed Characterization Experiment. By introducing the Stokes framework used in many other fields, the article opens the door to a large set of methods developed and used in other contexts but largely ignored in ocean acoustics.

## I. INTRODUCTION

Sound can be described with thermodynamic variables, the most common being pressure *p*, or kinematic flow variables such as particle acceleration **a** and velocity **v**. All of these continuous field variables represent, in the linear acoustic regime, small perturbations from their background values. For example, in a static fluid (a common approximation for sea water), *p* and **v** are related through the linearized Euler equation,

with *ρ* the ocean density and *t* the time.

Although the ocean acoustics literature tends to focus on the pressure *p*, which is clearly easier to measure and encumbered with fewer systematic errors, experimental techniques for inclusion of vector acoustic fields exist (e.g., Dahl and Dall'Osto, 2020a; D'Spain *et al.*, 1991; Gray *et al.*, 2016; Shchurov *et al.*, 2011). Along with defining sound energy flow for which *p* and **v** are in phase (active intensity), particle velocity notably carries information about the source and the environment and thus can be used as input for source localization (Hawkes and Nehorai, 2003; Thode *et al.*, 2010) and/or environmental estimation (Bonnel *et al.*, 2019a; Dahl and Dall'Osto, 2020b; Dall'Osto *et al.*, 2012; Ren and Hermand, 2012; Shi *et al.*, 2019). Particle velocity also emerges as an important field for fishes and crustaceans, given their use of **v** either exclusively or possibly in addition to *p* in auditory and environmental sensing processes (Popper and Hawkins, 2018).

The aim of the current article is to provide a physical and signal processing framework adapted to the study of particle velocity, with a specific focus on low-frequency (*f* < 500 Hz) propagation in shallow water (depth *D* < 200 m). In this context, the propagation is described by normal mode theory: the environment acts as a dispersive waveguide. The sound field is conveniently described as a sum of (dispersive) modes, each mode propagating with its own frequency-dependent group velocity [Jensen *et al.* (2011), Chap. 5]. A common assumption associated with normal mode propagation is that the waveguide depends exclusively on range *r* and depth *z* (it is rotationally invariant around the source). As a result, the particle velocity can be seen as a bi-dimensional vector and can be expressed using its horizontal and vertical components: $v=[vr,vz]T$. Note that the vector field is usually measured by a three-dimensional (3D) sensor with a vertical (*v _{z}*) and two horizontal components (usually denoted

*v*and

_{x}*v*). The latter (

_{y}*v*,

_{x}*v*) can be projected onto the source/receiver axis to create

_{y}*v*.

_{r}In this context, the particle velocity defines a *bivariate signal* (Flamant, 2018), where the interrelations between horizontal and vertical components encode relevant information about the oceanic medium. Unfortunately, conventional approaches to process bivariate signals such as rotary components (Gonella, 1972; Mooers, 1973) or classic bivariate analysis (Priestley, 1981) do not allow for straightforward interpretation of their geometric and physical properties. To circumvent this issue, a general framework for the analysis and processing of bivariate signals has been recently proposed (Flamant *et al.*, 2017, 2019). It unravels the key role played by the physical notion of polarization, usually encountered in optics (Born and Wolf, 1980) and radar (Lee and Pottier, 2017), to describe and understand the inner properties of bivariate signals. This framework also establishes the general relevance of Stokes parameters (Stokes, 1851)—a set of four real-valued observables widely used to describe polarization properties of light (Rubin *et al.*, 2019; Schaefer *et al.*, 2007; Stenflo, 2013)—for bivariate signals. Their extension to the case of elastic waves (denoted *elastodynamic Stokes parameters*) was proposed (Turner and Weaver, 1994) for the study of multiple scattering of ultrasounds in elastic media, but they were never used, to the best of our knowledge, for the description of waterborne sound particle velocity. Importantly, Stokes parameters are energetic quantities, meaning that they can be easily estimated from experimental bivariate signal measurements.

In ocean acoustics, sound polarization has mostly been described through complex intensity $I=pv*$ (Dahl and Dall'Osto, 2020a; Dall'Osto *et al.*, 2012; D'Spain *et al.*, 1991). Here, we propose to study the polarization of normal modes using the Stokes framework. Stokes parameters are directly derived from **v**, such that experimental measures do not require a coherent and collocated measure of *p* and **v**. Insightful polarization properties of normal modes are introduced. It is notably demonstrated that the polarization of individual modes does not depend on source position. It is also shown that the angle of the main axis of the polarization ellipse is proportional to modal attenuation; in a lossless waveguide, the polarization ellipse is horizontal. Basic tools to assess modal polarization, both in the frequency and in the time-frequency (TF) domains, are introduced: Stokes parameters can be easily derived from standard signal processing quantities [Fourier transforms (FTs) and spectrograms]. The proposed framework is illustrated on simulated as well as experimental particle velocity data, collected at sea during the 2017 Seabed Characterization Experiment (SBCEX17) (Wilson *et al.*, 2020).

The remainder of the article is organized as follows. Section II presents the physical background required for the article. It overviews the basics of normal mode propagation, with the main focus on particle velocity that follows the development in Dahl and Dall'Osto (2020a) prior to introducing the concept of the bivariate signal. Section III covers the signal processing background required for this article. It presents the bivariate signal context, overviews polarization concepts including the polarization ellipse, and introduces the Stokes parameters. Section IV constitutes the heart of the article. It makes the link between Secs. II and III. Modal Stokes parameters are derived and illustrated on a simulated scenario that mimics SBCEX17. Section V introduces the concept of polarization spectrograms, which extend the Stokes parameter to the TF domain and enable the visualization of modal Stokes parameters using a single vector sensor and an impulsive source. Polarization spectrograms for modal propagation are illustrated on simulated data as well as on SBCEX17 experimental data in Sec. VI. The article is concluded in Sec. VII, which provides a summary and discusses future opportunities. Two appendixes supplement the article. Appendix A presents a comparison between the proposed Stokes framework and vector metrics recently defined in Dahl and Dall'Osto (2020a). Appendix B provides an experimental observation of the degree of polarization (i.e., polarization variability).

## II. PARTICLE VELOCITY AND NORMAL MODES

Low-frequency acoustic propagation in shallow water is conveniently described using normal mode theory. Given a broadband source signal $\Omega (f)$ emitted at depth *z _{s}* in a range-independent waveguide, the pressure field

*p*received at depth

*z*after propagation of a range

*r*is given by [Jensen

*et al.*(2011), Chap. 5]

where *M* is the number of propagating modes, $\Psi n$ and *k _{n}* are, respectively, the modal depth function and the horizontal wavenumber of mode

*n*, $Q=ej\pi /4/8\pi \rho (zs)$ is a factor depending only on the density $\rho (zs)$ at depth

*z*, and $j2=\u22121$. In general, the wavenumber

_{s}*k*is a complex number and can be written as $kn=kn(r)\u2212j\beta n$, with $kn(r)$ and

_{n}*β*the real and imaginary parts of

_{n}*k*. In a lossless waveguide (i.e., with no attenuation), the modal wavenumber is real: $kn=kn(r)$ and $\beta n=0$.

_{n}To simplify notations, we define

so that the pressure associated with mode *n* is

The modal particle velocity is obtained by combining Eqs. (1) and (4). The horizontal component of particle velocity for mode *n* is

while its vertical component is

In the frequency domain, the particle velocity vector associated with mode *n* is thus $Vn(f)=[Vrn(f),Vzn(f)]T$. One obtains the modal particle velocity in the time domain, $vn(t)$, using an inverse FT. Formally,

where $F$ is the FT operator.

As stated previously, **v**(*t*)$=\Sigma vn(t)$ is said to be a bivariate signal, with components made up as a linear combination of modes.

Figure 1 illustrates the interpretation of the bivariate signal as a time-evolving vector. Figure 1(a) depicts a bivariate monochromatic signal, the elementary brick underpinning our analysis. Such a signal is said to be elliptically polarized, as its time-evolving trace in the two-dimensional (2D) (*v _{r}*,

*v*) plane is an ellipse, whose shape is governed by the interrelations between the amplitudes and phases of radial and vertical harmonic components. Figure 1(b) gives an example of a more sophisticated narrowband bivariate signal. Although the signal is still narrowband, its polarization is no longer constant: the signal's trace in the 2D (

_{z}*v*,

_{r}*v*) plane is now an ellipse with changing shape, size, and orientation. This is due to complex relationships between the amplitudes and phases of the signal's components. The main goal of bivariate signal analysis is to extract polarization information from such signals. The Stokes parameter framework presented in this paper is a convenient and physics-based solution to this problem.

_{z}## III. BIVARIATE SIGNAL DESCRIPTION

This section provides a short introduction to the framework for bivariate signal processing introduced and developed in Flamant (2018). Unlike existing approaches, this new framework enables straightforward interpretations of well-established signal processing tools in terms of the physical concept of (wave) polarization. Notably, Stokes parameters, a set of four real-valued energetic parameters widely used in polarization optics (Born and Wolf, 1980), lie at the core of the framework. The numerous physical interpretations enabled by the framework greatly ease the design, analysis, and processing of bivariate signals.

Here, for brevity, we only review the necessary ingredients and results of the framework, and we refer to the original papers for further details (Flamant *et al.*, 2017, 2018, 2019).

### A. The geometry of monochromatic bivariate signals

Let us consider first the most elementary bivariate signal, i.e., a monochromatic bivariate signal $q(t)$ of frequency *f*. This vector signal writes the following way:

with $ax,ay\u22650$ and $\phi x,\phi y\u2208[0,2\pi )$ the amplitude and phase of *x* and *y*, respectively. The two components *x* and *y* of **q** are univariate monochromatic signals with the same frequency *f*, but they may have different amplitude and phase. As a result, the full description of $q(t)$ in (8) requires four parameters: the amplitudes *a _{x}*,

*a*, and phases $\phi x$, $\phi y$ of the two components. The interrelations between these quantities precisely govern the geometric properties of the signal $q(t)$.

_{y}Figure 2(a) depicts the ellipse drawn over time by the bivariate signal **q**(*t*) in the *x*–*y* 2D plane. The ellipse trajectory is described by four real-valued parameters:

the

*amplitude*or*size*$\kappa \u22650$ of the ellipse;the

*phase*at origin $\phi \u2208[0,2\pi )$;the ellipse

*orientation*$\theta \u2208[\u2212\pi /2,\pi /2]$, representing the angle between the main axis of the ellipse and the horizontal axis;the

*ellipticity*angle $\chi \u2208[\u2212\pi /4,\pi /4]$, which characterizes the shape of the ellipse as well as the rotation direction in the ellipse: counterclockwise if $\chi \u22650$ and clockwise when $\chi <0$.

The first two parameters are classical as they correspond to the standard notion of amplitude and phase. The two remaining parameters are geometric and encode the *polarization* of the bivariate monochromatic signal. In particular, when *χ* = 0, the ellipse becomes a line segment: we say that $q(t)$ is linearly polarized. For $\chi =\xb1\pi /4$, the ellipse degenerates into a circle, so that $q(t)$ is said to be circularly polarized. Moreover, the signal is said to be horizontally polarized when *θ* = 0, whereas if $\theta =\xb1\pi /2$ it is said to be vertically polarized.

The *canonical* parameters $[\kappa ,\phi ,\theta ,\chi ]$ neatly encode the trajectory of the bivariate monochromatic bivariate signal $q(t)$. In other terms, they allow for a joint description of the properties of the univariate components *x* and *y*. In particular, we can rewrite $q(t)$ in Eq. (8) directly in terms of $\kappa ,\phi ,\theta ,\chi $ as (Flamant *et al.*, 2019),

We can also express the canonical parameters in terms of amplitudes *a _{x}*,

*a*and phases $\phi x,\phi y$ as

_{y}when $ax\u2260ay$. The case *a _{x}* =

*a*corresponds to circular polarization: here $\chi =\xb1\pi /4$ and

_{y}*θ*is undefined.

### B. Stokes parameters and the Poincaré sphere

To describe polarization, a popular alternative to the canonical geometric parameters introduced above consists of the Stokes parameters denoted by $S0,S1,S2,S3$. These four real-valued parameters benefit from being easily measured experimentally from intensity measurement, which made them particularly convenient for applications in optics (Rubin *et al.*, 2019; Schaefer *et al.*, 2007; Stenflo, 2013).

Considering again the bivariate monochromatic signal $q(t)$ given in Eq. (8), Stokes parameters can be expressed in terms of the ellipse parameters as

where $\Phi \u2208[0,1]$ is the *degree of polarization* of the signal, further discussed below. First, note that expressions (14)–(17) do not include the overall phase $\phi $, since Stokes parameters are energetic quantities, which, by definition, are insensitive to a global phase shift.

For the sake of generality, Eqs. (15)–(17) include the degree of polarization $\Phi $, a statistical measure of the dispersion of the polarization ellipse across multiple realizations (Flamant *et al.*, 2017) (see Fig. 3). If one assumes that the bivariate signal is generated by a stochastic (i.e., random) process, then the degree of polarization is used to quantify the stability of the polarization ellipse from one realization to another. When Φ = 1, the polarization ellipse is always strictly the same, and the signal is said to be fully polarized. When Φ = 0, the ellipse is fully random, and the signal is unpolarized. In intermediate cases ($0<\Phi <1$), the signal is said to be partially polarized.

Note that deterministic signals, such as $q(t)$ given by (8), are fully polarized signals (Φ = 1). Indeed, all the realizations of a deterministic signal are strictly identical; therefore, the polarization ellipse is the same for all the realizations. In this paper, we assume that bivariate signals are fully polarized unless stated otherwise. Section VII discusses briefly the case of partial polarization, and an example from field observations is given in Appendix B.

Stokes parameters are homogeneous to intensities, such that *S*_{0} describes pure energetic information while the three remaining parameters *S*_{1}, *S*_{2}, and *S*_{3} give geometric polarization properties in intensity units. It is convenient to remove the intensity dependence in $S1,S2,S3$ by defining normalized Stokes parameters as

Thanks to the Poincaré sphere of polarization states shown in Fig. 2(b), normalized Stokes parameters have a natural interpretation as the Cartesian coordinates of the vector described by spherical coordinates $(\Phi ,2\theta ,2\chi )$. Each point on the sphere is associated with a single polarization state. Thus, one can easily switch from one representation to another: for instance, the amplitude *a*, orientation *θ*, ellipticity *χ*, and degree of polarization $\Phi $ are obtained from the Stokes parameters using

### C. Spectral Stokes parameters for general bivariate signals

Until now, we have only considered the definition of Stokes parameters for a single bivariate monochromatic signal of frequency *f*. To generalize the Stokes formalism to generic, broadband bivariate signals, one needs to define frequency-dependent Stokes parameters, as explained below.

Spectral Stokes parameters simultaneously describe the energetic and polarization properties of a bivariate signal with respect to frequency. They are defined easily in terms of power spectral densities (PSDs). Let $Pxx(f)$ be the PSD of *x*(*t*), $Pyy(f)$ the PSD of *y*(*t*), and $Pxy(f)$ the cross-spectral density (CSD) between *x*(*t*) and *y*(*t*). Note that in practice, for a deterministic signal with a finite support (length *T*), the PSDs are trivially obtained using the FT. Let $X(f)=Fx(t)$ and $Y(f)=Fy(t)$. Then $Pxx(f)=X(f)X*(f)/T,\u2009Pyy(f)=Y(f)Y*(f)/T$ and $Pxy(f)=X(f)Y*(f)/T$. More generally, the PSD and CSD can always be defined as the FT of the autocorrelations and cross-correlations of the univariate signals *x*(*t*) and *y*(*t*).^{1} Spectral Stokes parameters are then defined as (Flamant *et al.*, 2017; Schreier and Scharf, 2010)

where $\u211c[.]$ and $\u2111[.]$ stand for real and imaginary parts. Just as before, for frequencies such that $S0(f)\u22600$, one defines the normalized Stokes parameters, with values between –1 and 1, as $s1(f)=S1(f)/S0(f),\u2009\u2009s2(f)=S2(f)/S0(f),\u2009\u2009s3(f)=S3(f)/S0(f)$.

## IV. MODAL POLARIZATION

### A. Theory

This section makes a link between the polarization parameters (Sec. III) and the modal propagation (Sec. II). The propagation environment is assumed to be known and non-fluctuating, and the modes are noiseless. As a result, considered signals are deterministic: the modes are fully polarized (Φ = 1). The context of partially polarized modes, which arises when the received signal is noisy and/or when the propagation environment is fluctuating, will be discussed in Sec. VII and in Appendix B.

For a given mode *n*, we make a link between the particle velocity [Eq. (7)] and the vector model for bivariate signals [Eq (8)]. This link is done through $x(t)=vrn(t)$ and $y(t)=vzn(t)$. In the frequency domain, assuming a particle velocity signal of length *T*,

One can easily show that

The normalized Stokes parameters are thus

Equations (34)–(36) are important results that characterize the modal polarization. First of all, it is reassuring to see that the (normalized) Stokes parameters are real numbers and that they do not depend on the signal length *T*. Also, they depend on the environment (through *k _{n}* and $\Psi n$) and on the receiver depth

*z*. However, they do not depend on the range

*r*or on the source depth

*z*. In other words, for a given receiver position, they are fully independent from the source position. This important behavior is obtained thanks to the normalization process (division by

_{s}*S*

_{0}). Remember that the results presented here have been derived in a range-independent waveguide. The derivation of the Stokes parameters for range-dependent waveguides, particularly when mode coupling occurs, is an exciting research question, but it is beyond the scope of this paper.

Finally, remember that in the specific case of a lossless waveguide, $kn=kn(r)$, so that $Pxy(f)$ is purely imaginary. This indicates that the spectral components of *x* and *y* are in phase quadrature (90° phase shift). Further, in the same scenario, $\beta m=0$ leads to $s2=0$. Using (20), this implies that *θ* = 0 and the major axis of the polarization ellipse is horizontal. Equivalently, this also means that its polarization state is located on the prime meridian of the Poincaré sphere (Fig. 2). In a lossless waveguide, the particle velocity of individual normal modes is polarized horizontally.^{2} In a general waveguide with attenuation, *s*_{2} is directly proportional to *β _{m}*.

### B. Example

This section illustrates the Stokes parameters on a simulated scenario that mimics SBCEX17 (Wilson *et al.*, 2020). The experiment took place on the “New England Mud Patch,” about 100 km south of Cape Cod, Massachusetts. A specificity of the environment is that the seabed features a first, thick layer of mud over more consolidated sediments. In this article, a notional geoacoustic model of the environment is considered:

water column: depth

*D*= 74.4 m, isovelocity sound speed profile $cw=1468.3$ m/s;layered seabed with a mud layer (thickness

*h*= 11 m, sound speed gradient from $cmudTOP=1460$ m/s to $cmudBOT=1550$ m/s, constant density $\rho mud=1.6$, and constant attenuation $\alpha mud=0.05$ dB/_{mud}*λ*);semi-infinite basement (sound speed

*c*= 1800 m/s, density_{b}*ρ*= 2, and attenuation $\alpha b=0.1$ dB/_{b}*λ*);source range

*r*= 16 km, source depth $zs=18.3$ m, receiver depth $zr=72.2$ m.

This specific simulated scenario has been chosen to reproduce the conditions of a particle velocity study by Dahl and Dall'Osto (2020a). In their article, Dahl and Dall'Osto derived four quantities from coherent combinations between pressure and particle velocity. Although they are based on the complex intensity $I=pv*$ rather than on just **v**, those quantities are closely related to the modal polarization. A thorough comparison between these quantities and the Stokes parameters is presented in Appendix A.

The simulated Stokes parameters for the first five modes are shown in Fig. 4 (continuous lines). The figure shows *S*_{0} and the three normalized Stokes parameters (*s*_{1}, *s*_{2}, *s*_{3}) for frequencies between 0 and 200 Hz. One clearly sees the effect of modal dispersion on the Stokes parameters: they depend both on mode number and frequency. Remember that $s2(f)\u221d\beta m$ was an important theoretical result from Sec. IV A. Since *β _{m}* is typically a very small number, so is

*s*

_{2}. Note here that the vertical scale for

*s*

_{2}is different from the one used for

*s*

_{1}and

*s*

_{3}.

To illustrate the polarization sensitivity to geoacoustic parameters, the Stokes parameters are also computed in an environment where the mud layer is replaced by sand, which follows an example given in Dahl and Dall'Osto (2020a). The only difference from the previous simulation is that $cmudTOP=cmudBOT=1600$ m/s. The resulting data are also shown in Fig. 4 (dotted lines), and these differ markedly from the previous dataset. This suggests that the Stokes parameters have a high sensitivity to the environment and may be good inputs for geoacoustic inversion. This is confirmed by looking at the group speed associated with the two different environments, also plotted in Fig. 4. Group speeds are usually derived from pressure data and are known to be good input data for environmental inversion (e.g., Bonnel *et al.*, 2021; Potty *et al.*, 2000). However, Stokes parameters are clearly more sensitive to environmental changes. Provided they can be correctly estimated from particle velocity data, the Stokes parameters should become excellent input for geoacoustic inversion.

## V. POLARIZATION SPECTROGRAMS

Figure 4 shows that the modal polarization is frequency-dependent. This, obviously, is intrinsically related to modal dispersion. In the following, we use TF analysis to study modal dispersion. Indeed, spectrograms are commonly used to visualize and estimate the dispersion of individual modes with application for geoacoustic inversions (Ballard *et al.*, 2014; Bonnel *et al.*, 2019a; Potty and Miller, 2020) or marine mammal localization (Bonnel *et al.*, 2014; Thode *et al.*, 2017).

One may thus wonder if the modal separation provided by TF analysis could be helpful to assess the modal polarization. The recent paper by Dahl and Dall'Osto (2020a) (mentioned in Sec. IV B and detailed in Appendix A) suggests that TF analysis is indeed a good tool to study polarization-related properties of modes. In this section, the concept of the polarization spectrogram, based on the Stokes framework presented above, is introduced and illustrated on simulated data. An experimental example will be provided in Sec. VI. The current section shows that the polarization spectrograms of particle velocity data, as measured by a single vector sensor, allow the visualization of modal polarization properties of individual modes.

### A. Theory

The Stokes framework, presented above, has been derived for monochromatic signals. However, it can easily be extended to study non-stationary deterministic signals in the TF domain. Indeed, one can build *polarization spectrograms* showing the TF distribution of the (normalized) Stokes parameters.

The polarization spectrogram theory is fully described in Flamant *et al.* (2019). In practice, polarization spectrograms can be computed very simply. The procedure is similar to what is presented in Sec. III C, except that the FT operator $F$ needs to be replaced by a short-time FT operator. The TF Stokes parameter $S0(t,f)$ is thus obtained by summing the spectrograms of the vertical and horizontal particle velocity, while $S1(t,f)$ is obtained by computing their difference. The parameters $S2(t,f)$ and $S3(t,f)$ are obtained by evaluating the real and imaginary part of the cross-spectrograms. Normalized TF Stokes parameters $s1(t,f),s2(t,f)$ and $s3(t,f)$] are obtained through normalization by $S0(t,f)$. The four TF Stokes parameters $S0(t,f)$, $s1(t,f)$, $s2(t,f)$, and $s3(t,f)$ will now be called “polarization spectrograms.” Polarization spectrograms can be trivially computed in most programming languages using any off-the-shelf TF toolbox. Note that the normalization by $S0(t,f)$ may be problematic if $S0(t,f)\u22430$. This issue is easily solved by setting a lower threshold *ϵ* on $S0(t,f)$. To do so, one finds all the TF points such that $S0(t,f)<\u03f5$ and then sets $S0(t,f)=\u03f5$ for all those points. The parameter *ϵ* is typically fixed as a small percentage (typically 0.1% or less) of the total energy of the signal and/or the maximum value attained by *S*_{0} in the TF plane. Note that a python toolbox dedicated to the TF analysis of bivariate signals is available with Flamant *et al.* (2019).

An important property of polarization spectrograms is that, for a single-component noise-free signal and no interference, the polarization spectrogram values at the signal's TF location (i.e., the ridge) exactly give the values of the underlying Stokes parameter (Flamant *et al.*, 2019). In more complex scenarios, if the polarization spectrograms are not overly contaminated by interference and/or noise, this property is still approximately true. Note that this behavior is fully similar to that of traditional spectrograms, which give the signal's energy along the TF ridge (up to a constant multiplicative factor that depends on spectrogram parameters).

As a reminder, in our modal propagation context, the theoretical TF location of mode *m* with frequency *f* is given by (Bonnel *et al.*, 2020)

with *r* the source/receiver range, $vm(f)$ the group speed of mode *m*, and $ts(f)$ the emission time of frequency *f* by the source. If the source is impulsive, all the frequencies are emitted at the same time, and $ts(f)=t0$ is constant. Note that modal separation in the TF domain increases when range *r* increases.

### B. Example

The simulated scenario from Sec. IV B is now used to illustrate polarization spectrograms. Simulated time series are computed using the first five propagating modes, assuming a perfectly impulsive source with frequencies between 0 and 200 Hz. Polarization spectrograms are computed using the procedure explained in Sec. V A. The threshold *ϵ* is chosen as 0.1% of the maximum of $S0(t,f)$. The obtained polarization spectrograms are presented in Fig. 5. To facilitate reading, the modal TF locations are shown as black curves on all the spectrograms.

The polarization spectrogram $S0(t,f)$ shows the TF distribution of the bivariate signal's energy. It can be interpreted as a traditional spectrogram for a univariate signal. Here, the simulated range is large enough for $S0(t,f)$ to show cleanly separated modes. The other polarization spectrograms $s1(t,f)$, $s2(t,f)$, and $s3(t,f)$ show the TF distribution of the normalized Stokes parameters. They enable the evaluation of individual mode polarization properties. From a signal processing point of view, the polarization properties fully encode the interdependencies between amplitudes and phases of the vector field, and the polarization spectrograms enable doing so for individual modes.

The polarization spectrogram $s1(t,f)$ shows that *s*_{1} is larger than 0.5 for all modes at all frequencies, with the notable exception of mode 5, for which $s1(t,f)$ drastically drops for frequencies below 150 Hz. This behavior is fully consistent with normalized Stokes parameters computed for individual modes, as shown in the top-right panel of Fig. 4. From Fig. 4, note that negative values of $s1(f)$ are predicted for mode 5 and at frequencies between $\u224375$ and 100 Hz; this is not visible on $s1(t,f)$ (Fig. 5) because this frequency band is highly attenuated.

The polarization spectrogram $s2(t,f)$ is quite puzzling at first look. Although we know $|s2|\u226a1$ (see Fig. 4), $s2(t,f)$ clearly shows high values between modes. This is due to interference resulting from the TF uncertainty, which spreads the modes outside of their theoretical locations (i.e., this interference is an artifact from the TF processing). However, for frequencies higher than $\u224375$ Hz, one sees that $s2(t,f)$ is indeed very small along the modal TF position. On the other hand, for frequencies less than $\u224375$ Hz, $s2(t,f)$ is relatively large even at the modal TF position. This is because, at those frequencies, actual interference between modes exists. This is shown in Fig. 5 by the black curves that cross each other. Physically, the Airy phase of mode *m*_{0} interferes with modes *m*, with $m<m0$. Note that the acoustic energy associated with the modal Airy phase is highly attenuated so that this behavior is unlikely to be resolved on most noisy and/or experimental data.

Finally, the polarization spectrogram $s3(t,f)$ shows clear polarization difference between modes. As an example, $s3(t,f)<0$ for mode 1 at all frequencies, while $s3(t,f)>0$ for modes 3–5. More interestingly, $s3(t,f)$ changes sign for mode 2 around 50 Hz. This detailed behavior is fully consistent with $s3(f)$, as computed for individual modes (see the bottom-right panel of Fig. 4).

## VI. EXPERIMENTAL APPLICATION

This section presents an experimental application of the Stokes parameter framework for data collected during SBCEX17.

### A. SBCEX17

SBCEX17 was a multi-institutional, multi-ship, multi-disciplinary effort that took place on the New England Mud Patch, about 110 km south of Cape Cod, in March/April 2017. Its main objective was to advance understanding of the acoustic properties of fine-grained sediments with clay, i.e., mud. An overview of SBCEX17 is provided in Wilson *et al.* (2020).

During SBCEX17, multiple acoustic receivers and sources were deployed, covering frequencies from about 10 Hz to 10 kHz. Of particular interest here is the Intensity Vector Autonomous Recorder (IVAR) deployed by the Applied Physics Laboratory, University of Washington (Dahl and Dall'Osto, 2020a) at location (40.48655° N; 70.63831° W). IVAR is a bottom moored vector sensor: it records the 3D particle velocity vector field at a given location in the water column, $\u223c1$ m above the seafloor.

In this article, we consider an excerpt of IVAR data that contains the signal from a distant combustive sound source (CSS) deployed by the Applied Research Laboratory, University of Texas. The specific CSS transmission occurred on March 18 around 16:19 UTC at a way-point called “station 35,” located at (40.49881° N; –70.45240° W), $\u223c16$ km away from IVAR. Three shots were successively emitted at depth $zs\u224320$ m; this paper considers only the first. The CSS signal is a low-frequency high-energy impulse, followed by several weaker replicas called “bubble pulses” (McNeese *et al.*, 2014). The source signal was measured during the experiment using a hydrophone hard-mounted on the CSS deployment frame, $\u223c1$ m away from the center of the source chamber. This signal will be used to perform source deconvolution.

### B. Data

As stated in Sec. VI A, IVAR measures the 3D particle velocity vector field. More specifically, IVAR measures two horizontal velocity components and a vertical one. To build the bivariate signal $v(t)=[vr(t),vz(t)]T$, it is required to project the two IVAR horizontal components into a single $vr(t)$, with the *r* axis pointing toward the source. This is done using simple geometrical rules, as explained in Eq. (12) of Dahl and Dall'Osto (2020a).

The received signal is further processed using a band stop filter between 86 and 88 Hz to remove contamination from internal interference (note that a high-pass filter starting at 25 Hz is also embedded in the recording system to prevent signal saturation). Last, source deconvolution is performed using the simple “water level” method (Clayton and Wiggins, 1976). The method description and its application to SBCEX17 CSS data are presented in depth in Bonnel *et al.* (2020), along with data examples and companion matlab code.

The preprocessing chain (horizontal projection, filtering, and source deconvolution) leads to an experimental bivariate signal $v(t)=[vr(t),vz(t)]T$ that represents the waveguide impulse response [i.e., as if $\Omega (f)=1$]. This signal is presented in Fig. 6. Since the signal contains several modes, its time-domain representation is not legible. Still, one notes that the peak value of $vz(t)$ is roughly an order of magnitude smaller than the peak value of $vr(t)$. This leads to a 2D trace in the *r*–*z* domain that looks like a horizontal elongated ellipse. This will be further quantified, on an individual mode basis, using polarization spectrograms.

### C. Polarization spectrograms

Experimental polarization spectrograms are shown in Fig. 7. The theoretical TF positions of the modes, computed using the environmental parameters from Sec. IV B, are also plotted as black curves. The polarization spectrograms were computed using an *ϵ* threshold equal to 1% of the maximum of $S0(t,f)$. This value was empirically chosen to provide visually good results. It is notably higher than the one used in Sec. V B to cope with experimental noise.

The experimental polarization spectrograms enable the visualization of the polarization properties of individual modes. These spectrograms are qualitatively similar to the simulated ones. Important features include high values of $s1(t,f)$ for modes 1–4 at all frequencies. Although $s3(t,f)$ is clearly contaminated by noise and interference, it shows a pattern similar to the one predicted in Fig. 5. Indeed, $s3(t,f)>0.5$ at most frequencies for modes 4 and 5, while it has smaller but positive values for mode 3. Modes 1 and 2 are not as legible, but mode 2 shows $s3(t,f)<0$ for frequencies between 50 and 100 Hz, as predicted by the simulation.

Note that a clear mismatch exists between experiment and simulation for fine details of the polarization spectrograms. As an example, the experimental data do not show $s1(t,f)<0$ for mode 3 and *f* < 100 Hz. This is to be expected since the simulation has been performed using a notional environmental model. Obtaining a true match between simulation and experiment would require a dedicated environmental inversion, which is beyond the scope of this paper. Nonetheless, the qualitative agreement between simulation and experiment clearly demonstrates the relevance of the proposed method.

## VII. DISCUSSION

The Stokes parameter framework presented in this article enables a full description of the polarization of underwater sound. It differs from existing underwater vector acoustics works mentioned previously in that it focuses exclusively on the particle velocity **v** rather than on the complex intensity $I=pv*$ that requires a concurrent and collocated measure of both pressure *p* and particle velocity **v**. Since *p* and **v** have to be measured with two different sensors (usually a hydrophone for *p* and an accelerometer for **v**), the proposed framework simplifies the experimental task associated with measuring underwater sound polarization.

Restricting the scope of the study to modal propagation, Dahl and Dall'Osto recently defined a set of “modal vector metrics” that describes modal polarization (Dahl and Dall'Osto, 2020a). As stated above, the Dahl and Dall'Osto metrics are different from the Stokes parameters because they are derived from $I=pv*$. Another important difference is that the Dahl and Dall'Osto metric definition requires some (limited) knowledge about the environment, notably the water sound speed and density at the receiver location. Last, these metrics are fully defined using physics-based arguments. While they are very informative about the modal polarization, there is no path to determine if those metrics form a complete description of the polarization. On the other hand, since the Stokes framework is based on signal processing arguments, it fully describes the polarized signal. We note that indeed the Dahl and Dall'Osto metrics are very similar to the Stokes parameters, and the formal link between the two is presented in Appendix A.

Whether considering Dahl and Dall'Osto metrics or the Stokes parameters, the reader may wonder if modal polarization has practical applications. In Sec. IV B, we suggested that the Stokes parameters are highly sensitive to the environment (see Fig. 4) and thus may be used as input for geoacoustic inversion. This idea was explored for the Dahl and Dall'Osto metrics in Bonnel *et al.* (2019a). A similar study (not presented here for the sake of concision) was run for the Stokes parameters, and similar results were obtained. The modal polarization parameters (Stokes or Dahl and Dall'Osto metrics) are more sensitive to the seafloor properties than the modal group speeds, a classical input for geoacoustic inversion (Ballard *et al.*, 2014; Bonnel *et al.*, 2019b; Potty *et al.*, 2000). As a result, Stokes parameters appear to be promising input data for upcoming inversion studies. An example is the recent study by Dahl and Dall'Osto (2021) involving what they refer to as circularity, or *s*_{3} in this study, for geoacoustic inversion of underwater ship noise from SBCEX17.

Also, it is important to come back to the notion of degree of polarization $\Phi $, which was introduced in Sec. III B. The degree of polarization can be defined as a statistical measure of dispersion of the polarization ellipse across multiple realizations (Flamant *et al.*, 2017; Schreier and Scharf, 2010). In the context considered in this article, the signal is noise-free and the environment is not fluctuating. As a result, the signal is fully deterministic and modal polarization does not change. Therefore, the signal is fully polarized and Φ = 1. This, obviously, will never be true in an experimental context. Let us consider a classical tomographic setup with a fixed receiver and a fixed source. If the source emits recurrent signals, each received signal can be considered as a realization of the underlying stochastic oceanic process (Colosi, 2016). One can then estimate the Stokes parameters and derive the degree of polarization $\Phi $ using Eq. (22). It is expected that $\Phi $ would be representative of environmental fluctuations. The exact link between $\Phi $ and environmental fluctuation regimes (unsaturated, partially saturated, and fully saturated) needs to be determined. However, interesting perspectives arise in using the Stokes parameters to quantify environmental fluctuations. An important result from Dahl and Dall'Osto (2020a) is revisited in Appendix B to experimentally illustrate the concept of degree of polarization of individual modes. More broadly, the degree of polarization can be used to track fluctuations of either the whole signal or individual arrivals (modes or rays). This opens new research avenues, notably for long range propagation in deep water. Indeed, in this context, water column fluctuations highly impact the signal, and new methods are required to better link the physical oceanography and the acoustics.

Finally, using the Stokes parameters for practical marine applications calls into question our ability to properly estimate them from a particle velocity signal. For generic bivariate signals, Eqs. (23)–(26) show that spectral Stokes parameters can be efficiently obtained by combination of conventional non-parametric spectral density estimators (e.g., periodogram, multitaper) [see Flamant *et al.* (2017) for details]. In the context of particle velocity signals, estimation methods for individual modes should further take into consideration the near-horizontal polarization, i.e., $|s2(f)|\u226a1$. Other physical constraints may be also included, such as parametric modeling of the spectral dispersion of Stokes parameters. Also, polarization properties of modal interference can be examined [see Dahl and Dall'Osto (2020b, 2021) for examples with ship noise], which opens the door to interesting questions about the link between polarization and the waveguide invariant [Jensen *et al.* (2011), Chap. 5]. As a result, the practical estimation of Stokes parameters defines key challenges, in terms of both physical interpretability and robustness to noise. The development of dedicated estimation procedures is required to ensure the full exploitation of the information gathered in Stokes parameters. The proposed framework opens the door to novel physics-based bivariate signal processing methods for ocean acoustics.

## ACKNOWLEDGMENTS

This work was supported by the Direction Générale de l'Armement (France) and the Office of Naval Research. The authors warmly thank the anonymous “Reviewer 2” for excellent comments and suggestions. The reviewer contribution notably includes accounting for waveguide attenuation with notable impact on Eqs. (32) and (35) as well as suggestions that led to the derivation of Eqs. (A5), (A6), (A8), and (A9).

### APPENDIX A: OTHER MODAL VECTOR METRICS

This appendix makes a comparison between the Stokes parameters and the modal vector metrics defined by Dahl and Dall'Osto (2020a).

##### 1. Definitions

Using coherent combination between pressure and velocity channels, Dahl and Dall'Osto (2020a) define four quantities derived from the vector acoustic measurements: modal phase speed, circularity, normalized vertical reactive intensity, and the depth-dependent mode speed of energy. All these quantities depend on the mode number *n* and on the frequency. The first quantity, the phase speed, is not related to polarization and thus not interesting here. The other three quantities (hereinafter referred to as Dahl and Dall'Osto metrics) are defined as follows:

- circularity,(A1)$\Theta n(f)=2\u2111{Vrn(f)Vzn*(f)}|Vrn(f)|2+|Vzn(f)|2;$
- normalized vertical reactive intensity,(A2)$Qzn(f)=\rho cw\u2111{Pn(f)Vzn*(f)}|Pn(f)|2;$
- depth-dependent mode speed of energy,(A3)$uen(f)=2\u211c{Pn(f)Vrn*(f)}|Pn(f)|2/(\rho cw2)+\rho (|Vrn(f)|2+|Vzn(f)|2).$

In the above formula, *ρ* and *c _{w}* are, respectively, the density and sound speed in the water column. Both are assumed to be constant. Note that $uen(f)$ has previously been defined as “transport velocity” in other papers (e.g., D'Spain

*et al.*, 1991; Schiffrer and Stanzial, 1994).

##### 2. Comparison with Stokes parameters

An interesting question arises in comparing the Dahl and Dall'Osto metrics with the Stokes parameters. First of all, it is obvious that

The link between the other two Dahl and Dall'Osto metrics and the Stokes parameters requires further derivation. By combining Eqs. (4) and (5), one can show that $pn(f)=\rho cn(f)Vrn(f)$, with $cn(f)=2\pi f/kn$. In a lossless waveguide, $kn=kn(r)$ is a real number, and *c _{n}* is the modal phase speed. In a general lossy waveguide,

*c*is merely a quantity convenient for derivation, $cn=(2\pi f/|kn|2)(kn(r)+j\beta n)=2\pi f(kn*/|kn|2)$, with $*$ the complex conjugation operator. One can then show that

_{n}As stated above, in a lossless waveguide, $cn(f)$ is the modal phase speed. It can be rewritten as $cn(f)=cw/\u2009cos[\u03d1n(f)]$, with $\u03d1n(f)=arccos[kn(f)cw/(2\pi f)]$ the mode angle. In such a case, one can show that

which relates the depth-dependent mode speed of energy and the Stokes parameters through the mode angle. If one further makes a low-grazing angle approximation, i.e., $\u03d1n(f)\u22430$ and $cn(f)=cw$, the relationship becomes even simpler,

Similarly, one can derive that, in a general lossy waveguide, the normalized vertical reactive intensity

Following the same approximation as above, one can show that, in a lossless waveguide,

A notable difference between Eqs. (A8) and (A9) is that the latter does not contain $s2(f)$. This is because, in a lossless waveguide, $s2(f)=0$ for all modes at every frequency. Further, under a low-grazing angle approximation, one obtains the simple relationship

Interestingly, Eqs. (A7) and (A10) are very good approximations even for actual propagating modes ($\u03d1n(f)>0$) in a lossy waveguide. The three Dahl and Dall'Osto metrics (A1)–(A3) are simulated in the environment described in Sec. IV B. The lossless low-grazing angle Stokes approximations [Eqs. (A7) and (A10)] as well as *s*_{3} [Eq. (A4)] are also simulated in the same environment. A comparison between the two sets of metrics is presented in Fig. 8: the match is very good for all modes at all frequencies, except for the depth-dependent speed of energy of modes 1 and 2 at low frequencies. This is probably due to the low-grazing angle approximation, which does not hold for these modes at those frequencies.

### APPENDIX B: EXPERIMENTAL OBSERVATION OF THE DEGREE OF POLARIZATION

This appendix revisits an important result from Dahl and Dall'Osto (2020a) and uses it to experimentally illustrate the concept of degree of polarization $\Phi $ (see Sec. III B and Fig. 3). In their paper, Dahl and Dall'Osto consider experimental data collected during SBCEX17. The dataset consists of 22 different source signals recorded on a vector sensor (IVAR). The 22 acoustic sources are at different locations, all with locations at least 5 km away from the IVAR to yield sufficient TF dispersion to identify separate mode arrivals. The received signals are used to estimate the polarization metrics $\Theta n(f)$, $Qzn(f)$, and $uen(f)$ (see Appendix A).

The whole process leads to 22 different experimental estimates of the polarization metrics. However, since the metrics are independent from the source location (exactly like the normalized Stokes parameters), the 22 estimates are roughly similar. This is illustrated for the circularity $\Theta n(f)$ in Fig. 9. This specific metric is chosen here because $\Theta n(f)=s3(f)$. Similar results are presented in Dahl and Dall'Osto (2020a) for the other metrics $Qzn(f)$ and $uen(f)$.

Figure 9 shows the estimation statistics in terms of median and 25th and 75th percentiles. While some dispersion exists, the estimated circularity is very consistent for all the source signals. The variability may be due to two factors: the ambient noise and environmental fluctuations. These two factors contribute in making estimated polarization properties different when using different source signals. If one assumes ergodicity (source signals recorded at different times have the same underlying statistics) and spatial homogeneity (the environment is perfectly range- and azimuth-independent), this constitutes an experimental observation of the degree of polarization $\Phi $.

Formally, note that estimating $\Phi $ from the dataset used in Dahl and Dall'Osto (2020a) is not straightforward and is beyond the scope of this paper. In short, it would require assumptions about the variability of $\Phi $. If $\Phi $ is assumed to be spatially homogeneous, then it could be estimated by averaging normalized Stokes parameters (further assuming ergodicity). On the other hand, if the environment is spatially variable, then it is likely that $\Phi $ depends on source position. In this case, it must be estimated using Stokes parameters rather than the normalized ones. In such a scenario, normalized Stokes parameters (as well as Dall and Dall'Osto metrics) would depend on source position.

^{1}

This works both for deterministic and stochastic signals.

^{2}

It is possible for *s*_{2} = 0 along with $\theta \u22600$ while $\chi =\xb1\pi /4$, but in this case, the polarization ellipse is a circle, and orientation is meaningless.