## Abstract

Observation of polarization modulation at the output of a submarine link, extracted from a standard coherent telecom receiver, can be used to monitor geophysical events such as sea waves and earthquakes occurring along the cable. We analyze the effect of birefringence perturbations on the polarization at the output of a long-haul submarine transmission system, and provide analytical expressions instrumental to understanding the dependence of the observed polarization modulation on the amplitude and spatial extension of the observed events. By symmetry considerations, we show that in standard single mode fibers with random polarization coupling, if polarization fluctuations are caused by strain or pressure, the relative birefringence fluctuations are equal to the relative fluctuations of the polarization averaged phase. We finally show that pressure induced strain is a plausible explanation of the origin of polarization modulations observed in a long submarine link. The presented analysis paves the way for the transformation of transoceanic fiber optic links during operation into powerful sensing tools for otherwise inaccessible geophysical events occurring in the deep ocean.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

The depths of the oceans are the places on Earth most immune to anthropic perturbations. They are ideal places where coherent Earth motions can be studied and monitored. However, ocean depths are hardly accessible. Buoys and cabled observatories provide monitoring and sensing of limited regions, but give only point-wise information on the ocean and seafloor conditions. The only man-made artifacts that could be encountered in the deep ocean are fiber optics cables belonging to long-haul optical communication systems. Equipped with sensors, they would be of great help in understanding the Earth motions because they lie in the most quiet environment available on Earth [1]. Unfortunately, equipping undersea cable with sensors and collecting data from them entail not trivial and costly modifications of the cable design, which would unavoidably interfere with the cable mission of transmitting information.

Coherent optical communication systems require for their operation the identification in real time of the matrix that characterizes the input–output relation of field polarization, the so-called Jones matrix [2,3]. The Jones matrix of a long fiber link depends on the cumulative effect of the birefringence of the fiber that connects the transmitter to the receiver. While in terrestrial systems the temporal modulation of the Jones matrix is significant, and mainly affected by anthropic activities and electromagnetic perturbations such as lightning, in submarine systems, the modulation of the Jones matrix is minimal, and mostly affected by the surrounding ocean environment. These data can be extracted with virtually no extra hardware, and would be of utmost value in the study of seafloor dynamics if the effect of environmental perturbations on the Jones matrix is properly characterized.

The knowledge of the Jones matrix enables the operation of virtual interferometric measurements reproducing the injection of a continuous wave with fixed polarization and observation of polarization modulation at the output of the line. This procedure effectively turns a long-haul transmission system into a polarization interferometer [4]. The use of undersea cables to monitor geophysical events, based on the detection of perturbations induced by geophysical events in the deep sea on the absolute phase of an ultra-stable laser of sub-hertz linewidth, was already demonstrated by Marra *et al.* in [5]. Polarization has the advantage over phase of being intrinsically stable, and decoupled from the noise of the absolute phase of the laser. The intrinsic stability of polarization enables high-precision sensing by use of the transmit and local oscillator lasers of the coherent transceiver, which are telecom-grade lasers of tens of kilohertz linewidth. Distributed acoustic sensing (DAS), a technique based on the detection of the phase of Raleigh backscattering originated from distributed fiber locations, is also an extremely valuable technique because it provides high spatial resolution [6]. However, its spatial range is limited and its use mainly restricted to dark fiber plans.

The use of polarization, however, poses significant complications to the interpretation of the results. First of all, while the dynamic range of phase measurements is in principle unlimited, to preserve linearity, the dynamic range of polarization measurements is limited to observation of differential phase shifts between the two polarizations much smaller than $2\pi$. In terms of optical path length, this means that the measurement is limited to events that produce a difference in the optical path lengths on the two polarizations smaller than the laser wavelength, 155 µm. Being the length of the submarine systems of the order of thousands of kilometers, this implies that to be detected with no ambiguity, the events need to produce changes in the two polarizations’ optical path length difference much smaller than one part over ${10^{12}}$. This means that geophysical sensing is possible only with a transmission system where the Jones matrix of the system or, equivalently, the output polarization for a fixed input polarization, is extremely stable. Furthermore, an added complication is that even in the absence of perturbation, the two independent polarizations of the transmitted light randomly couple along the fiber because of static random birefringence [7]. Finally, the relation between polarization modulation at the fiber output and the environmental perturbations that affect birefringence is not an obvious one. The development of a detailed model of the dependence of polarization fluctuations on birefringence perturbations in a long fiber link is therefore a crucial step in turning observation of polarization modulations into a powerful monitoring and sensing tool, because it will enable both the interpretation of output polarization and understanding of the conditions under which sensing is feasible.

In this paper, we present a model for the state of polarization (SOP) modulation at the output of a fiber link. We show that the mean square deviation of the Stokes vector is proportional to the square of the local polarization mode dispersion (PMD) coefficient of the fiber. We also show that the random orientation of the fiber birefringence causes an incoherent addition of polarization modulation along the fiber path, and this in turn produces a linear growth of the mean square deviation of the Stokes vector with the length of the portion of the cable affected by the perturbation. This linear dependence increases the dynamic range of the measurement, and explains the recent observations [4] that even earthquakes of magnitude larger than ${M_w} = 7$ did not produce on a 10,000 km long submarine link polarization modulations that covered the entire Poincaré sphere. This also suggests that the dynamic range of these measurements is large enough to be used in practical sensing of geophysical events of a wide energy range. As a byproduct of this investigation, we show by symmetry considerations that the way birefringence in standard single mode fibers (SSMFs) with randomly coupled birefringence is affected by strain and pressure is not arbitrary, but is unequivocally determined by the way the polarization averaged phase changes. This property enables in principle the prediction of the observed SOP modulation from DAS measurements. We finally show that pressure induced strain is at the origin of SOP modulations caused in a long submarine cable by sea waves, and it is also a plausible explanation for those observed during an earthquake of magnitude ${M_w} = 5.3$.

## 2. ANALYSIS

The equation describing the evolution of the PMD vector $\vec \tau$ with propagation distance $z$ is [7]

There is a strong analogy between ${\vec \beta _\omega}$ and $\Delta \vec \beta (z,t)$. The first parameter quantifies the sensitivity of the birefringence to a change of the light center frequency $\omega$, the second to a change of environmental conditions. Being $\beta = \omega \delta n/c$, if we neglect the frequency dependence of the difference $\delta n$ of the effective refractive index between the two polarizations, the vector ${\vec \beta _\omega}$ is to a good approximation parallel to $\vec \beta$, and hence we can set $\vec \beta = \omega {\vec \beta _\omega}$ [8]. When the birefringence perturbation is caused by hydrostatic pressure or strain, $\Delta \vec \beta (z,t)$ is also parallel to $\vec \beta$ for the axial symmetry of pressure and strain.

If we use a frame rotating with the birefringence at $\omega$ as in [7], each vector is replaced by $\vec v^\prime (z) = {\textbf{R}^{- 1}}(z,0)\vec v(z)$, where $\textbf{R}(z,0)$ is a rotation operator defined by ${\rm d}\textbf{R}(z,0)/{\rm d}z = \vec \beta (z) \times \textbf{R}(z,0)$. Doing so, we remove the static, $z$ dependent rotations caused by $\vec \beta (z)$, and obtain at frequency $\omega$

and where we removed the prime in the symbol of the rotated vectors, meaning from now on that $\vec \beta$, ${\vec \beta _\omega}$, $\Delta \vec \beta (z,t)$, $\vec \tau$, and $\vec s$ are defined in the new frame rotating with $\vec \beta = \omega {\vec \beta _\omega}$, the static birefringence at frequency $\omega$. In the following, we will use arguments based on symmetry to identify the direction of various vectors. These symmetries involve relative orientations between vectors, and are valid in the rotated frame as well as in the original frame, because the relative angle between vectors is preserved by rotations. Finally, the action of the rotation operator $\textbf{R}(z,0)$, which is the concatenation of many independent rotations, makes the rotated vectors $\vec \beta$, ${\vec \beta _\omega}$, and $\Delta \vec \beta (z,t)$ isotropically distributed, although they originally described linear birefringence vectors, hence belonging to the plane ${s_3} = 0$.The use of a reference frame rotating with the static birefringence is equivalent to using a perfect compensator for the static birefringence at center frequency $\omega$, and implies that in the absence of perturbations, the output polarization at this frequency is equal to the input polarization, namely, that $\vec s(z) = \vec s(0) = {\vec s_0}$. If we define $\vec s = \Delta \vec s + {\vec s_0}$, we obtain to first order for $\Delta s \ll 1$ (from now on, we will use for any vector $\vec v$ the standard notation $v = |\vec v|$)

which shows that only the component of $\Delta \vec \beta (z,t)$ orthogonal to ${\vec s_0}$ is effective, and that $\Delta \vec s$ belongs to the plane orthogonal to ${\vec s_0}$. Solving this equation, we obtain where $\Delta {\vec \beta _ \bot}(z,t) = \Delta \vec \beta (z^\prime ,t) \times {\vec s_0}$ is the component of $\Delta \vec \beta (z,t)$ perpendicular to ${\vec s_0}$, rotated clockwise by 90° around ${\vec s_0}$. The time dependence of the vector $\Delta \vec s$ characterizes SOP modulation. Equation (6), valid for small perturbations, shows that the deviation $\Delta \vec s$ from the unperturbed position of $\vec s$ is equal to the component of the (isotropically distributed) birefringence perturbations orthogonal to the input SOP ${\vec s_0}$, integrated from zero to $z$. Consequently, the spectrum of $\Delta \vec s$ is equal to the spectrum of the integrated birefringence perturbations, and filtering the first is equivalent to filtering the latter. This reality has the important consequence that it enables spectral analysis of the monitored processes (earthquakes or sea waves) by a spectral analysis of SOP modulations if the coupling between the cable and the seafloor is linear. Of course, the analysis of the relation between polarization deviation and birefringence perturbations is complicated by the random nature of the magnitude and orientation of $\Delta \vec \beta (z^\prime ,t)$, which we are now going to characterize.The correlation function of the frequency derivative of the birefringence is stationary in space:

Equations (7) and (8) may also describe cases in which, besides $\langle \Delta {\beta ^2}\rangle$, also the PMD coefficient $\langle {\beta ^2}\rangle /{\omega ^2}$ is $z$ dependent over a much longer length scale than the length scale of $g(z - z^\prime)$. From now on, we will implicitly assume that both $\langle {\beta ^2}\rangle /{\omega ^2}$ and $\langle \Delta {\beta ^2}\rangle$ are weakly dependent on $z$.

Equation (6) yields

Using the relation between the average PMD square and the average square PMD $\langle {\tau ^2}\rangle = (3\pi /8){\langle \tau \rangle ^2}$, and defining ${\kappa ^2} = {\langle \tau \rangle ^2}/z$, we may write ${\langle {\tau ^2}\rangle _z} = (3\pi /8){\kappa ^2}$. Using also $\omega = 2\pi \lambda /c$, where $c$ is the speed of light in vacuum and $\lambda$ the wavelength, Eq. (14) can be rewritten as

Let us now show, in the next section, a property of randomly coupled SSMFs that is the consequence of its average circular symmetry, which will be extremely useful in the identification of possible sources of SOP modulation at the output of a long fiber link.

## 3. PHASE AND SOP SENSITIVITY TO STRAIN OF SSMF

One of the most successful fiber-optic based technique for geophysical sensing is DAS [6]. This technique is based on the time-resolved measurement of the phase of the Rayleigh backscattered signal, which is modulated by the time-dependent interaction of the fiber with the environment. Recently, a number of measurement campaigns have been reported that used DAS on standard (armored or lightwave) submarine cables [10–12]. To compare with DAS measurements, we need to establish a relation between the average phase accumulated during propagation and the SOP modulations at the output. This is the goal of the following analysis.

The propagation of the electric field $\vec E$ over a length $\ell \ll {L_F}$ of a birefringent fiber is described by ${\rm d}\vec E/{\rm d}z = i \textbf{H} \vec E$, where, in the basis of the local linear eigenpolarizations, $\textbf{H}$ is a $2 \times 2$ diagonal matrix with entries the wave vectors of the two eigenpolarizations ${\beta _{1,2}} = 2\pi {n_{1,2}}/\lambda$, with ${n_{1,2}}$ the corresponding effective refractive indices [3]. If we use the expansion $\textbf{H} = {\beta _0}\textbf{I} + (\beta /2) {\vec \sigma _1}$, where $\textbf{I}$ is the identity matrix and ${\sigma _1}$ the diagonal Pauli matrix with $1$ and ${-}1$ on the main diagonal, then ${\beta _0} = ({\beta _1} + {\beta _2})/2$ is the polarization averaged wave vector, and $\beta = {\beta _1} - {\beta _2}$ is the modulus of the birefringence vector $\vec \beta = (\beta ,0,0)$, parallel to the first axis of the Stokes space because we have used the basis of fiber eigenpolarizations [3]. We have ${\beta _0} = 2\pi n/\lambda$ and $\beta = 2\pi \delta n/\lambda$, where $n = ({n_1} + {n_2})/2$ is the polarization averaged effective refractive, and $\delta n = {n_1} - {n_2}$ is the difference of the effective refractive indices between the two eigenpolarizations.

Assume now that axial strain is applied to a birefringent fiber. The key observation is that axial symmetric strain cannot rotate the birefringence axes. Strain can only change the eigenvectors ${\beta _1}$ and ${\beta _2}$ of the quantities $\Delta {\beta _1} = 2\pi \Delta {n_1}/\lambda$ and $\Delta {\beta _2} = 2\pi \Delta {n_1}/\lambda$, where $\Delta {n_{1,2}}$ are caused by the strain optic effect and by a change in the mode profile. Strain also affects the propagating distance. Both effects contribute to a change of the polarization averaged phase and of the differential phase. Let us analyze the effect on the average phase first.

The phase accumulated by light propagating over a length $\ell$ is ${\varphi _0} = {\beta _0}\ell$. Strain affects ${\varphi _0}$ through a change of $\ell$ and of ${\beta _0}$ and, hence it produces the phase shift $\Delta {\varphi _0} = \Delta {\beta _0}\ell + {\beta _0}\Delta \ell$, where $\Delta {\beta _0} = (\Delta {\beta _1} + \Delta {\beta _2})/2$. This phase shift can be written as $\Delta {\varphi _0} = {\epsilon _0}{\varphi _0}$, where ${\epsilon _0} = \Delta {\beta _0}/{\beta _0} + \Delta \ell /\ell = \Delta ({\beta _0}\ell)/({\beta _0}\ell)$. If $n\ell$ is the polarization averaged optical path length, than ${\epsilon _0} = \Delta (n\ell)/(n\ell)$ so that ${\epsilon _0}$ is also the relative variation of the polarization-averaged optical path length. The quantity ${\epsilon _0}$ can therefore be considered as an *optical* strain, as opposed to the material strain ${\epsilon _0} = \Delta \ell /\ell$, which is the relative variation of the fiber’s geometrical length. The optical and material strains are proportional, and their relation is characterized in the DAS literature [12] by ${\epsilon _0} = \xi {\epsilon _0}$, where $\xi$ is the photoelastic scaling factor.

The three-dimensional Stokes vector $\Delta \vec \beta = (\Delta \beta ,0,0)$, where $\Delta \beta = \Delta {\beta _1} - \Delta {\beta _2} = 2\pi \Delta (\delta n)/\lambda$, quantifies the effect of strain on the birefringence vector $\vec \beta$. The differential phase between the two eigenpolarizations is $\varphi = \beta \ell$, and its change is $\Delta \varphi = \Delta \beta \ell + \beta \Delta \ell = \epsilon \varphi$, where $\epsilon = \Delta \beta /\beta + \Delta \ell /\ell = \Delta (\beta \ell)/(\beta \ell)$. Also in this case, we have that $\epsilon = \Delta (\delta n \ell)/(\delta n \ell)$ is the relative variation of the difference between the optical path lengths of the two eigenpolarizations.

A remarkable property is that in SSMFs with randomly coupled birefringence, $\epsilon = {\epsilon _0}$. This property is equivalent to $\Delta \beta /\beta = \Delta {\beta _0}/{\beta _0}$, which is in turn equivalent to $\Delta {\beta _1}/{\beta _1} = \Delta {\beta _2}/{\beta _2}$, as it may be verified by direct substitution. The last property is valid if strain affects equally the two eigenpolarizations, which is a well-founded assumption for the average circular symmetry of SSMFs with randomly coupled birefringence (it is not valid for high birefringence fibers, where circular symmetry is intentionally broken).

In the following analysis, it is convenient to keep the propagation length constant, adding to $\Delta \beta$ a contribution that accounts for the change of $\ell$, defining $\Delta \beta ^\prime = \Delta \beta + (\Delta \ell /\ell)\beta$. $\Delta \vec \beta ^\prime $ being parallel to $\vec \beta$, we also have $\Delta \vec \beta ^\prime = \Delta \vec \beta + (\Delta \ell /\ell)\vec \beta$. Using this definition, we have $\Delta \vec \beta ^\prime = \epsilon \vec \beta$. In the following, we will always deal with $\Delta \vec \beta ^\prime $ and not use $\Delta \vec \beta$ any longer, so that we will remove the prime in $\Delta \vec \beta ^\prime $ for simplicity of notation.

The equality ${\epsilon _0} = \epsilon$ is an important byproduct of this analysis. It allows to deduce, in a SSMF with random polarization coupling, the sensitivity of polarization to strain and pressure from the phase response to strain and pressure of the same fiber.

## 4. INTERPRETATION OF THE SOP MODULATION

Using $\Delta \vec \beta = \epsilon \vec \beta$ in Eq. (6), and assuming that the time-dependent strain $\epsilon$ varies slowly and deterministically with $z$, we obtain

where ${\vec \beta _ \bot}(z) = \vec \beta (z) \times {\vec s_0}$ is the component of $\vec \beta (z)$ perpendicular to ${\vec s_0}$, rotated clockwise by 90° around ${\vec s_0}$. Using $\langle \Delta {\beta ^2}\rangle = {\epsilon ^2}\langle {\beta ^2}\rangle$ in Eq. (15), we also obtainLet us now use Eq. (16) to investigate the temporal correlations of $\Delta \vec s$. In that equation, ${\vec \beta _ \bot}(z)$ is $z$ dependent, static, and random, while $\epsilon (t)$ is deterministic, time dependent, and weakly dependent on $z$. The birefringence deviations $\Delta \vec \beta (z,t) = \epsilon (t)\vec \beta (z)$ at a given position $z$ and times $t$ and $t^\prime $ differ only for their lengths and for being parallel or antiparallel to one another. Therefore, the left-hand side of Eq. (8) can be readily generalized to yield $\langle \Delta \vec \beta (z,t) \,\cdot$$\Delta \vec \beta (z^\prime ,t^\prime)\rangle = g(z - z^\prime) \epsilon (t) \epsilon (t^\prime)\langle {\beta ^2}\rangle$, where $g(z - z^\prime)$ is the spatial correlation function defined by Eq. (7). Using this expression, and again the delta function approximation for $g(z - z^\prime)$, we may derive the following equation for the (non-stationary) temporal correlation function of SOP deviations, which generalizes Eq. (17):

A possibly significant source of birefringence fluctuations not related to strain is fiber twist (curvature can also add birefringence, but the curvatures that can be realistically generated by external perturbations are too small to have an effect in submarine cables, which are almost straight). In this case, to first order, $\Delta \vec \beta = 2\vec \alpha \times \vec \beta$, where $\alpha$ is the angle of rotation in real space (the angle of rotation in Stokes space is twice the angle of rotation in real space [3]). Neglecting the small stress-optic rotation coefficient, in the original frame, $\vec \alpha$ transforms linear birefringence into a rotated linear birefringence; hence it is orthogonal to $\vec \beta$, so that $\Delta \beta = 2\alpha \beta$. In the rotating frame, angles are preserved so that, assuming $\vec \alpha$ deterministic and slowly varying with $z$, we have $\langle \Delta {\beta ^2}\rangle = 4{\alpha ^2}\langle {\beta ^2}\rangle$, and hence Eq. (15) becomes

## 5. DISCUSSION

Equation (15) (and the ones derived from them) has important implications. First of all, it shows that the modulation of output polarization, for a given birefringence modulation, is proportional to the local PMD of the fiber, which is the only significant parameter. This implies that PMD quantifies the ability of a system to be used for sensing. If it is too large, earthquakes of larger magnitudes are likely to saturate the polarization modulation for energetic events. The fact that PMD is a parameter usually well characterized in installed systems facilitates the interpretation of the measurements.

Equation (15) gives also useful information on the expected dependence of polarization modulations on the spatial extension of the perturbed region of the fiber. When a fiber of length $z$ and constant birefringence axis (i.e., constant direction of $\Delta \vec \beta$) is perturbed, to first order, the output polarization changes by $\Delta \vec s = \int_0^z \Delta \vec \beta {\rm d}z^\prime \times {\vec s_0}$, so that $\Delta {s^2} = {({\int_0^z \Delta {\beta _ \bot}{\rm d}z^\prime})^2}$, where $\Delta {\vec \beta _ \bot}$ is the component of the birefringence perturbation $\Delta \vec \beta$ orthogonal to the input polarization, grows quadratically with distance. Conversely, Eq. (12) shows that the random orientation of the fiber birefringence axes produces a linear growth with distance of the mean square polarization modulation. This reality on one hand implies a sub-linear increase in the amplitude of SOP modulation with the extension of the region affected by perturbation, and on the other hand, being a function less peaked than its square, tends to enhance in the integral over $z$ the contribution of regions with more intense birefringence perturbation, making the sensing more point-wise. The smoother increase with the length of the section of the fiber affected by the perturbations enhances the stability of the link against birefringence perturbations and increases the dynamic range of the measurement when this deviation is used for sensing. This observation may also explain why, in the measurement campaign reported in [4], even earthquakes of the greatest magnitudes (${M_w} = 7$ and more) never showed such a strong polarization modulation as to cover the whole Poincaré sphere.

The incoherent nature of the interaction has also another important implication. Let us assume for the sake of illustration that the perturbation is caused by strain, and Eq. (17) applies. Equation (17) shows that the signal power is proportional to the integral of the strain variation square, not to the integral square of the strain variation, as it would be in the case of coherent interaction. The dependence of polarization fluctuations on the square of the local strain implies that positive and negative strains excited at different locations do not average out. This property makes sensing based on polarization somehow complementary to that in [5] based on phase. With that technique, the measurement is sensitive to the integral of the strain over the fiber length $\Delta \varphi = \int_0^z \Delta {\beta _0} {\rm d}z^\prime = (2\pi /\lambda)\int_0^z n \epsilon {\rm d}z^\prime.$ When there is a strong spatial inhomogeneity of $\epsilon$ with positive and negative regions, then $\Delta \varphi$ gets averaged. This may happen, for instance, when detecting earthquakes with epicenters at close distance to the cable. In this case, in the expression of the power of the signal $\Delta {\varphi ^2} = (2\pi /\lambda {)^2}\int_0^z \int_0^z n(z^\prime) \epsilon (z^\prime)n(z^{\prime \prime}) \epsilon (z^{\prime \prime}){\rm d}z^\prime {\rm d}z^{\prime \prime} $, the values outside the main diagonal average, out and the only surviving terms will be those with $z^\prime \simeq z^{\prime \prime} $. Under these circumstances, the detection of phase and polarization produces similar responses. Because of this effect, the detection with phase of earthquakes with epicenters at close distances will have an attenuated response, whereas those at larger distances, exciting a uniform strain on the fiber, will have an enhanced response because of the coherence of the phase accumulation.

## 6. COMPARISON WITH OBSERVATIONS

Polarization fluctuations generated by ocean waves and earthquakes have been recently detected in an $L = 10{,}500\; {\rm km}$ long submarine cable from Los Angeles to Valparaiso, the Curie cable [4]. The PMD coefficient of the cable ranges from $0.01\;{\rm ps}/\sqrt {{\rm km}}$ and $0.07\;{\rm ps}/\sqrt {{\rm km}}$, with an average of $0.03\;{\rm ps}/\sqrt {{\rm km}}$. As mentioned in the Introduction, the transponder of a coherent transmission system tracks the output polarization at tens of gigahertz rate. We extracted from the transponder, downsampled at 56 ms sampling period, the reconstructed Jones matrices, and we derived from them the output polarizations that correspond to a fixed input polarization. We then processed the obtained sequence of virtual output Stokes vectors with the following procedure. First, the Stokes parameters were averaged over a moving time window of 200 s, and then renormalized to unit length. Second, a slowly varying time-dependent rotation was applied to the reference frame in Stokes space such that the moving average was always centered on the North Pole of the Poincaré sphere, of coordinates (0,0,1). In this slowly rotating reference frame, the measured Stokes parameters represent the polarization deviations caused by perturbations faster than 200 s. This procedure was key to filter out the slow, long-term drift of the polarization caused by disturbances, e.g., slow thermal drifts, in the terminal stations, although added also an effective high pass filter that prevented us from being sensitive to processes whose dynamics are below about 5 mHz. The rotation that moves the point representing the average polarization to the North Pole is not unique, and we chose the one moving the point on a big circle of the Poincaré sphere. The class of rotations that transform one point of the sphere into another can be decomposed as the concatenation of a rotation on a big circle and a rotation around the final point. Therefore, the SOP points obtained with different choices all share the same absolute values of SOP deviations from the North Pole; they differ only for the orientation around the average value.

Figure 1 shows the spectrogram of $\Delta {s^2}$ obtained from 1 Jun 2020 to 12 Jul 2020, displaying dispersive wave packets around 0.06 Hz, each lasting for a few days. The timings of the wave packets coincide well with the primary and secondary microseism pairs observed at coastal seismic stations located in the North America region of the cable, and are attributed to excitations in the primary microseism band caused from ocean swells in the Los Angeles region [4]. Microseism signals at coastal sites are related to ocean swells produced by distant storms. The fact that the double-frequency secondary microseism, which is the seismic waves produced by wave–wave interactions, is not observed on the spectrum of $\Delta \vec s$, suggests that the dispersive wave packets are caused by seafloor pressure perturbations from ocean swells in shallow water, in the vicinity of the Los Angeles terminal of the cable. The absence of low-frequency noise down to 0.02 Hz should be noted. Again, this makes the use of polarization for sensing somewhat complementary to interferometric methods based on the absolute phase [5,13], which are instead noisier in the low-frequency region, because of being sensitive to the $1/{f^2}$ noise caused by the random walk of the laser phase induced by amplified spontaneous emission noise.

Figure 1 (and also Fig. S6 in [4]) shows that the power spectral density of $\Delta \vec s$ caused by the ocean swells in the Curie cable is of the order of ${10^{- 3}}\;{\rm Hz}^{- 1}$. This value, once integrated over a bandwidth of about 0.01 Hz, gives $\Delta {s^2} \simeq {10^{- 5}}$. Modulation of the SOP is likely, in this case, to be directly caused by pressure. Linear gravity wave theory dictates that the pressure decreases with depth $h$ as $\Delta P = \Delta {P_0}/\cosh ({k_w}h)$, where ${k_w} = 2\pi /{\lambda _w}$ is the wave vector of the ocean wave. The depth of the first 12 km of the Curie cable is about $h = 100\;{\rm m}$. The frequency of the wave is $\omega = \sqrt {g{k_w}\tanh ({k_w}h)}$ so that, being $\omega = 2\pi 0.06\; {\rm rad/s}$ (see Fig. 1), we have ${k_w} = 0.0158\;{\rm m}^{- 1}$, corresponding to ${\lambda _w} = 2\pi /{k_w} = 405\;{\rm m}$. Assuming $h = 100\;{\rm m}$ and $\Delta {P_0} = 20\;{\rm kPa}$ (corresponding to 2 m high sea waves), we obtain $\Delta P = \Delta {P_0}/\cosh ({k_w}h) = 8\;{\rm kPa}$. Using in Eq. (19) the observed value for the average of polarization modulation $\langle \Delta {s^2}\rangle {= 10^{- 5}}$, for the PMD coefficient the value $\kappa = 0.03\;{\rm ps}/\sqrt {{\rm km}}$, and assuming that the pressure modulation is applied over the 12 km length where the depth is less than 100 m, gives for the coefficient ${\epsilon _P}$ the estimate ${\epsilon _P} = 3.5 \cdot {10^{- 9}}\;{\rm Pa}^{- 1}$.

The strain induced by pressure on an installed armored cable has been characterized with DAS in the measurement campaign reported in [12]. In that paper, the curves reported in Fig. 4(d) show that the pressure spectral profile (gray curve) and the material strain spectral profile measured by DAS (red curve) are nearly proportional from ${10^{- 3}}\;{\rm Hz}$ to 10 Hz, with a coefficient of proportionality (calculated from the peak values $|{\varepsilon _0}({f_{{\rm max}}}{)|^2}{= 10^{- 2.2}}\unicode{x00B5} {\epsilon ^2}/{\rm Hz}$ and $|\Delta P({f_{{\rm max}}}{)|^2}{= 10^{14.7}}\unicode{x00B5}{\rm Pa}^2/{\rm Hz}$) ${\varepsilon _{0,P}} = {\varepsilon _0}/\Delta P \simeq 3.6 \cdot {10^{- 9}}\;{\rm Pa}^{- 1}$. The optical strain is obtained through multiplication by the photoelastic scaling factor $\xi = 0.78$ [12], giving ${\epsilon _{0,P}} \simeq 2.8 \cdot {10^{- 9}}\;{\rm Pa}^{- 1}$. This value is equal within the experimental uncertainties to that characterizing the sensitivity of polarization to pressure of the Curie cable, ${\epsilon _P} = 3.5 \cdot {10^{- 9}}\;{\rm Pa}^{- 1}$. These values are about one order of magnitude higher than those typical of fibers not embedded in a cable [14,15]. The mechanical characteristics of the petroleum jelly, and especially its low bulk modulus [14], may play a role in the enhancement of the pressure sensitivity of a cabled fiber. The increase with pressure of the local deviations from circular symmetry due to random inhomogeneities of the width of the polymeric coating may enhance the sensitivity of the fiber birefringence to pressure modulations. Overall, the fact that the estimated value of ${\epsilon _P}$ is similar to the measured one of ${\epsilon _{0,P}}$, although in a different cable, validates the equality, theoretically predicted in Section 3, ${\epsilon _0} = \epsilon$. It also shows that if the polarization modulation is caused by strain (that is, excluding twist), the result of measurements with polarization can be predicted by using in the integral of Eq. (17) $\epsilon = {\epsilon _0}$, where ${\epsilon _0}$ is obtained from DAS measurements.

We note that the measurements of Fig. 4(d) in [12] were obtained from a single armored light cable laid on the seafloor, whereas the first 12 km of the Curie cable are buried. However, it was found in [10] that the armoring of the cable or its burial have little effect (less than 15%) on the cable response (see Fig. S3 in that paper; the first 2.1 km of the cable are buried and the rest have different types of armoring). These findings, and other evidences reported in [11], also show that submarine cables are indeed permeable to hydrostatic pressure, and that it is likely that petroleum jelly plays an important role in transmitting the outside pressure modulation to the optical fibers.

Figure 2 shows the two components of $\Delta \vec s$, i.e., $\Delta {s_1}$ and $\Delta {s_2}$, after realignment of the Stokes vector parallel to the third axis of the Stokes space, for the ${M_w} = 5.3$ submarine earthquake that occurred on 4 August 2020 at 09:30:50 offshore Mexico with epicenter at 35 km to the closest point of the Curie cable. The polarization modulation is $\Delta {s^2} \simeq 0.01$. If we assume that the source of the birefringence perturbation is strain, and that strain is applied over an extension of 100 km and again use $\kappa = 0.03\;{\rm ps}/\sqrt {{\rm km}}$, the observed value of $\langle \Delta {s^2}\rangle = 0.01$ used in Eq. (17) requires $\epsilon = 3.1 \cdot {10^{- 4}}$. If the source of perturbation is the direct action of pressure on the fiber, setting $\epsilon = {\epsilon _P}\Delta P = 3.1 \cdot {10^{- 4}}$ and using the value ${\epsilon _P} = 3.5 \cdot {10^{- 9}}\;{\rm Pa}^{- 1}$ returns for the amplitude of the pressure modulation $\Delta P \simeq 88\;{\rm kPa}$, which is a reasonable value for pressure waves excited by earthquakes of similar magnitudes. For earthquake excitations, however, other mechanisms for SOP modulation cannot be excluded.

## 7. IMPLICATIONS OF SOP STOCHASTICITY

Let us now discuss the implications of the nondeterministic response of the SOP technique. We have already mentioned in Section 2 that $\Delta {s_1}$ and $\Delta {s_2}$, the components of $\Delta \vec s$, are zero-mean independent identically distributed Gaussian variables in the ensemble of all static fiber realizations. This ensemble of fibers is not accessible, however, as we can follow the time evolution of the only fiber realization available. A similar problem arises with the PMD vector, and it is solved by repeating the measurement of the PMD vector over a wide bandwidth and using the ergodicity in frequency of the PMD vector to gather sufficient statistics [16,17]. In principle, we can follow the same approach here by the simultaneous detection of the same event over multiple channels. Leaving the analysis of this option to future studies, let us concentrate on the time traces of the SOP modulation available on a single channel. In this case, the temporal correlation function given by Eq. (20) is non-stationary, and hence the SOP is not time ergodic, so that the we cannot replace ensemble averages with time averages. Some information on the expected statistics of the signal can, however, be inferred from the previous analysis.

Let us consider again the small modulation case, when Eq. (16) is valid. This equation shows that $\Delta \vec s$ has the same bandwidth as $\epsilon (t)$ integrated over the whole link and weighted by random fiber birefringence. If the spatial dependence of $\epsilon (t)$ were a rectangular function of $z$, i.e., constant over the section interested by the perturbation and zero elsewhere, then $\Delta \vec s$ would be a randomly oriented two-dimensional vector proportional to $\epsilon (t)$ with a Rayleigh distributed amplitude. In this case, the temporal trace of polarization on the Poincaré sphere would be a straight line, and the two traces of $\Delta {s_1}$ and $\Delta {s_2}$ one the replica of the other with in general unequal amplitudes related to the orientation of the straight line. In general, however, we observe traces of $\Delta \vec s$ that are different and cover an entire, slightly asymmetrical, region of the Poincaré sphere around the unperturbed Stokes vector. This is because $\epsilon (t)$ is never an ideal square function of $z$. The traces reported in Fig. 2, for instance, show that in the initial stage, the average amplitude of $\Delta {s_2}$ is slightly smaller than that of $\Delta {s_1}$ and then tends to equalize on average. Most earthquake detections observed in the campaign in [4] show traces of $\Delta {s_1}$ and $\Delta {s_2}$ that have similar but nonidentical amplitudes. This is because $\epsilon (t)$, although possibly dominated by the region closer to the earthquake epicenter and with the highest coupling, is spatially dependent in the region interested by the perturbation.

This behavior is quantitatively described by the temporal correlation function of Eq. (20), which shows that the faster $\epsilon (t)$ changes with $z$ in the region interested by the perturbation, the shorter the correlation time, which is the width of the temporal correlation function. Samples of $\Delta \vec s$ spaced by more than the correlation time are independent realizations of the Rayleigh variable. Consequently, when the system effectively samples a statistically significant number of independent realizations of the Rayleigh distribution of $\Delta \vec s$ over the time of the earthquake, the random nature of the SOP over the time scale of the processes of interest (whose spectrum goes from hundredths of hertz to a few hertz and lasts for hundreds of seconds) gets partially alleviated by self-averaging, especially for longer times. The self-averaging does not affect the spectrum of $\Delta \vec s$, which therefore reproduces quite faithfully the spectrum of $\epsilon (t)$. These considerations explain why the spectra of microseism and of earthquakes extracted from SOP modulation (see [4]) compare well with those obtained by techniques with a deterministic response such as DAS [6,12].

## 8. SENSITIVITY

In Fig. 1, we resolve deviations $\Delta \vec s$ from the unperturbed output polarization $\Delta {s^2} \simeq {10^{- 5}}$, which correspond to an angular deviation of $\vec s$ over the Poincaré sphere by $\theta \simeq 3 \cdot {10^{- 3}}\;{\rm rad}$. A full rotation of the Poincaré sphere is equivalent to a phase-shift between two eigenpolarizations of $2\pi\;{\rm rad}$, which is in turn generated by an optical path length difference of one wavelength, $\lambda = 1.55 \cdot {10^{- 6}}\;{\rm m}$. Our virtual interferometric measurement is therefore sensitive to an optical path length difference of $\Delta (\delta n \ell) = \lambda \theta /(2\pi) = 8 \cdot {10^{- 10}}\;{\rm m}$. If we define the sensitivity as the ratio of the minimum detectable variation of the optical path length difference to the total optical path length, $n\ell = 1.5 \cdot {10^7}\;{\rm m}$ for the Curie cable, the sensitivity is about $\Delta (\delta n \ell)/(n\ell) = 5 \cdot {10^{- 17}}$. Notice that this sensitivity is achieved with transmit and local oscillator lasers of hundreds of kilohertz linewidth. This is possible because polarization is decoupled from the laser phase noise, which is common to the two polarizations, and consequently, the assessment of the phase difference between the two polarizations is not limited by the linewidth of either the transmit or local oscillator laser.

## 9. CONCLUSION

In this paper, we analyzed the effect of birefringence perturbations on the output polarization of long-haul submarine transmission systems. The observation of SOP modulation at the output of a system, extracted from standard coherent telecom receivers, can be used to monitor geophysical events occurring in the ocean. We showed that the strain caused by modulation of the hydrostatic pressure caused by sea waves is at the origin of the polarization modulation observed in a long fiber link. We analyzed the polarization modulation induced by a magnitude 5.3 earthquake, and showed that the pressure excited by the earthquake is also a possible explanation of the observed SOP modulation. The analytical expressions obtained shed light on the role of the local PMD of the fiber in determining the sensitivity to birefringence modulations, showing that portions of the link with higher PMD give a higher contribution to the observed modulation of the SOP. We gave the scaling of the observed modulation with the amplitude and spatial extension of the detected events. We obtained an expression of the temporal correlation function of SOP modulation, which was used to discuss the implications of the random nature of SOP modulation induced by random static fiber birefringence. Finally, by symmetry considerations, we showed that in SSMFs with random polarization coupling, if polarization fluctuations are caused by strain or pressure, the relative birefringence fluctuations are equal to the relative fluctuations of the polarization averaged phase. This result has been used to give an order of magnitude estimate of the expected SOP fluctuations from standard submarine cables whose phase responses have been characterized using DAS.

## Funding

Ministero dell’Istruzione, dell’Università e della Ricerca (PRIN 2017 - FIRST); CIPE (INCIPICT); Gordon and Betty Moore Foundation.

## Acknowledgment

A. Mecozzi acknowledges financial support from the Italian Government under Cipe resolution n. 135 (21 Dec. 2012), project INnovating City Planning through Information and Communication Technologies (INCIPICT) and the Ministero dell’Istruzione, dell’Università e della Ricerca. Z. Zhan and J. Castillo Castellanos are funded by the Gordon and Betty Moore Foundation.

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

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

## REFERENCES

**1. **B. M. Howe, B. K. Arbic, J. Aucan, C. R. Barnes, N. Bayliff, N. Becker, R. Butler, L. Doyle, S. Elipot, G. C. Johnson, F. Landerer, S. Lentz, D. S. Luther, M. Müller, J. Mariano, K. Panayotou, C. Rowe, H. Ota, Y. T. Song, M. Thomas, P. N. Thomas, P. Thompson, F. Tilmann, T. Weber, and S. Weinstein, “Smart cables for observing the global ocean: science and implementation,” Front. Mar. Sci. **6**, 424 (2019). [CrossRef]

**2. **S. J. Savory, “Digital filters for coherent optical receivers,” Opt. Express **16**, 804–817 (2008). [CrossRef]

**3. **J. P. Gordon and H. Kogelnik, “PMD fundamentals: polarization mode dispersion in optical fibers,” Proc. Natl. Acad. Sci. USA **97**, 4541–4550 (2000). [CrossRef]

**4. **Z. Zhan, M. Cantono, V. Kamalov, A. Mecozzi, R. Müller, S. Yin, and J. C. Castellanos, “Optical polarization-based seismic and water wave sensing on transoceanic cables,” Science **371**, 931–936 (2021). [CrossRef]

**5. **G. Marra, C. Clivati, R. Luckett, A. Tampellini, J. Kronjäger, L. Wright, A. Mura, F. Levi, S. Robinson, A. Xuereb, B. Baptie, and D. Calonico, “Ultrastable laser interferometry for earthquake detection with terrestrial and submarine cables,” Science **361**, 486–490 (2018). [CrossRef]

**6. **M. R. Fernández-Ruiz, M. A. Soto, E. F. Williams, S. Martin-Lopez, Z. Zhan, M. Gonzalez-Herraez, and H. F. Martins, “Distributed acoustic sensing for seismic activity monitoring,” APL Photon. **5**, 030901 (2020). [CrossRef]

**7. **J. P. Gordon, “Statistical properties of polarization mode dispersion,” in *Polarization Mode Dispersion* (Springer, 2005), pp. 52–59.

**8. **A. Galtarossa, L. Palmieri, M. Schiano, and T. Tambosso, “Measurement of birefringence correlation length in long, single-mode fibers,” Opt. Lett. **26**, 962–964 (2001). [CrossRef]

**9. **M. Cantono, V. Curri, A. Mecozzi, and R. Gaudino, “Polarization-related statistics of Raman crosstalk in single-mode optical fibers,” J. Lightwave Technol. **34**, 1191–1205 (2016). [CrossRef]

**10. **A. Sladen, D. Rivet, J. P. Ampuero, L. De Barros, Y. Hello, G. Calbris, and P. Lamare, “Distributed sensing of earthquakes and ocean-solid earth interactions on seafloor telecom cables,” Nat. Commun. **10**, 5777 (2019). [CrossRef]

**11. **N. J. Lindsey, T. C. Dawe, and J. B. Ajo-Franklin, “Illuminating seafloor faults and ocean dynamics with dark fiber distributed acoustic sensing,” Science **366**, 1103–1107 (2019). [CrossRef]

**12. **H. Matsumoto, E. Araki, T. Kimura, G. Fujie, K. Shiraishi, T. Tonegawa, K. Obana, R. Arai, Y. Kaiho, Y. Nakamura, T. Yokobiki, S. Kodaira, N. Takahashi, R. Ellwood, V. Yartsev, and M. Karrenbach, “Detection of hydroacoustic signals on a fiber-optic submarine cable,” Sci. Rep. **11**, 2797 (2021). [CrossRef]

**13. **C. Clivati, A. Tampellini, A. Mura, F. Levi, G. Marra, P. Galea, A. Xuereb, and D. Calonico, “Optical frequency transfer over submarine fiber links,” Optica **5**, 893–901 (2018). [CrossRef]

**14. **B. Budiansky, D. C. Drucker, G. S. Kino, and J. R. Rice, “Pressure sensitivity of a clad optical fiber,” Appl. Opt. **18**, 4085–4088 (1979). [CrossRef]

**15. **L. Schenato, A. Galtarossa, A. Pasuto, and L. Palmieri, “Distributed optical fiber pressure sensors,” Opt. Fiber Technol. **58**, 102239 (2020). [CrossRef]

**16. **M. Karlsson and J. Brentel, “Autocorrelation function of the polarization-mode dispersion vector,” Opt. Lett. **24**, 939–941 (1999). [CrossRef]

**17. **M. Shtaif and A. Mecozzi, “Study of the frequency autocorrelation of the differential group delay in fibers with polarization mode dispersion,” Opt. Lett. **25**, 707–709 (2000). [CrossRef]