Polarization sensing using submarine optical cables

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 ﬁbers with random polarization coupling, if polarization ﬂuctuations are caused by strain or pressure, the relative birefringence ﬂuctuations are equal to the relative ﬂuctuations of the polarization averaged phase. We ﬁnally 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 ﬁber optic links during operation into powerful sensing tools for otherwise inaccessible geophysical events occurring in the deep ocean. ©2021OpticalSocietyofAmericaunderthetermsoftheOSAOpenAccessPublishingAgreement


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π . 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.

ANALYSIS
The equation describing the evolution of the PMD vector τ with propagation distance z is [7] where β(z) is the static fiber birefringence. The equation that describes the evolution with z of the unit Stokes vector s representing the polarization state of the field is where β(z, t ) is a time dependent perturbation of the birefringence, and t = t − nz/c is the retarded time, with c the speed of light in vacuum and n the refractive index of glass. Since the birefringence dynamics we are interested in are seconds or fractions of a second, much slower than the transit time of links up to 10,000 km long, which is about 50 ms, in the following, we will approximate everywhere t with t.
There is a strong analogy between β ω and β(z, t). The first parameter quantifies the sensitivity of the birefringence to a change of the light center frequency ω, the second to a change of environmental conditions. Being β = ωδn/c , if we neglect the frequency dependence of the difference δn of the effective refractive index between the two polarizations, the vector β ω is to a good approximation parallel to β, and hence we can set β = ω β ω [8]. When the birefringence perturbation is caused by hydrostatic pressure or strain, β(z, t) is also parallel to β for the axial symmetry of pressure and strain.
If we use a frame rotating with the birefringence at ω as in [7], each vector is replaced by v (z) = R −1 (z, 0) v(z), where R(z, 0) is a rotation operator defined by dR(z, 0)/dz = β(z) × R(z, 0). Doing so, we remove the static, z dependent rotations caused by β(z), and obtain at frequency ω and where we removed the prime in the symbol of the rotated vectors, meaning from now on that β, β ω , β(z, t), τ , and s are defined in the new frame rotating with β = ω β ω , the static birefringence at frequency ω. 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 R(z, 0), which is the concatenation of many independent rotations, makes the rotated vectors β, β ω , and β(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 ω, and implies that in the absence of perturbations, the output polarization at this frequency is equal to the input polarization, namely, that s (z) = s (0) = s 0 . If we define s = s + s 0 , we obtain to first order for s 1 (from now on, we will use for any vector v the standard notation v = | v|) which shows that only the component of β(z, t) orthogonal to s 0 is effective, and that s belongs to the plane orthogonal to s 0 . Solving this equation, we obtain where β ⊥ (z, t) = β(z , t) × s 0 is the component of β(z, t) perpendicular to s 0 , rotated clockwise by 90 • around s 0 . The time dependence of the vector s characterizes SOP modulation. Equation (6), valid for small perturbations, shows that the deviation s from the unperturbed position of s is equal to the component of the (isotropically distributed) birefringence perturbations orthogonal to the input SOP s 0 , integrated from zero to z. Consequently, the spectrum of 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 β(z , t), which we are now going to characterize. The correlation function of the frequency derivative of the birefringence is stationary in space: where we used the normalization g (0) = 1, and in the last equality, we used the relation where L F is the birefringence coherence length, with L F of the order of meters [8]. Although the birefringence perturbation is, rigorously, a non-stationary process because the strength of the perturbation is generally nonuniform along the link, the length scale of its variations is of the order of hundreds of meters, much longer than L F , so that the correlation function of the birefringence perturbation β(z, t) can be written as where β 2 is assumed to be weakly dependent on z. We assume that the normalized correlation function g (z − z ) of β and β are the same, because the rotation around the same static fiber birefringence randomizes the directions of both β and β. Equations (7) and (8) may also describe cases in which, besides β 2 , also the PMD coefficient β 2 /ω 2 is z dependent over a much longer length scale than the length scale of g (z − z ). From now on, we will implicitly assume that both β 2 /ω 2 and β 2 are weakly dependent on z.
Equation (6) yields If we now use the isotropy of β(z, t) and Eq. (8), we obtain and from Eqs. (3) and (7), Replacing g (z − z ) with its delta function approximation g (z − z ) = 2L F δ(z − z ), Eqs. (10) and (11) become Taking z derivative of Eq. (13), τ 2 z = (2L F /ω 2 ) β 2 , solving for L F and inserting into Eq. (12), we obtain where τ 2 z is the local average square PMD coefficient in ps 2 /km. For small deviations, where the first order theory is valid, s is a Rayleigh distributed random variable because it is the length of a two-dimensional vector that is the sum of many uncorrelated isotropically distributed vectors. The components s 1 and s 2 of s are, for the central limit theorem, zero-mean independent and identically distributed Gaussian variables. We will return to the implications of the statistics of s 1 and s 2 in Section 7.
Using the relation between the average PMD square and the average square PMD τ 2 = (3π/8) τ 2 , and defining κ 2 = τ 2 /z, we may write τ 2 z = (3π/8)κ 2 . Using also ω = 2π λ/c , where c is the speed of light in vacuum and λ the wavelength, Eq. (14) can be rewritten as Equation (15) is the main result of this paper. Its foundations lie on the well-established PMD theory, which has been confirmed by decades of experience. Once inverted, it allows for a given measured s 2 the estimation of the integral at the right-hand side, which depends on the strength of the perturbation. This equation shows that the points along the link where the PMD coefficient is larger have a higher sensitivity to birefringence perturbations. Equation (15) becomes less accurate for s 2 approaching 2/3 (which corresponds to an isotropic distribution over the Poincaré sphere). When s 2 is of the order of 2/3, although the probability density function of Eq. (4) can be analytically derived within the delta function approximation of g (z − z ) [9], the absence of a linear relation with the perturbations of the fiber birefringence makes the use of SOP modulations as a sensing tool questionable.
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.

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][11][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 E over a length L F of a birefringent fiber is described by d E /dz = iH E , where, in the basis of the local linear eigenpolarizations, H is a 2 × 2 diagonal matrix with entries the wave vectors of the two eigenpolarizations β 1,2 = 2π n 1,2 /λ, with n 1,2 the corresponding effective refractive indices [3]. If we use the expansion H = β 0 I + (β/2) σ 1 , where I is the identity matrix and σ 1 the diagonal Pauli matrix with 1 and −1 on the main diagonal, then β 0 = (β 1 + β 2 )/2 is the polarization averaged wave vector, and β = β 1 − β 2 is the modulus of the birefringence vector β = (β, 0, 0), parallel to the first axis of the Stokes space because we have used the basis of fiber eigenpolarizations [3]. We have β 0 = 2π n/λ and β = 2π δn/λ, where n = (n 1 + n 2 )/2 is the polarization averaged effective refractive, and δ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 β 1 and β 2 of the quantities β 1 = 2π n 1 /λ and β 2 = 2π n 1 /λ, where 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 is ϕ 0 = β 0 . Strain affects ϕ 0 through a change of and of β 0 and, hence it produces the phase shift ϕ 0 = β 0 + β 0 , where β 0 = ( β 1 + β 2 )/2. This phase shift can be written as . If n is the polarization averaged optical path length, than 0 = (n )/(n ) so that 0 is also the relative variation of the polarization-averaged optical path length. The quantity 0 can therefore be considered as an optical strain, as opposed to the material strain 0 = / , 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 0 = ξ 0 , where ξ is the photoelastic scaling factor.
The three-dimensional Stokes vector β = ( β, 0, 0), where β = β 1 − β 2 = 2π (δn)/λ, quantifies the effect of strain on the birefringence vector β. The differential phase between the two eigenpolarizations is ϕ = β , and its change is . Also in this case, we have that = (δn )/(δn ) 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, = 0 . This property is equivalent to β/β = β 0 /β 0 , which is in turn equivalent to β 1 /β 1 = β 2 /β 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 β a contribution that accounts for the change of , defining β = β + ( / )β. β being parallel to β, we also have β = β + ( / ) β. Using this definition, we have β = β. In the following, we will always deal with β and not use β any longer, so that we will remove the prime in β for simplicity of notation.
The equality 0 = 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.

INTERPRETATION OF THE SOP MODULATION
Using β = β in Eq. (6), and assuming that the time-dependent strain varies slowly and deterministically with z, we obtain where β ⊥ (z) = β(z) × s 0 is the component of β(z) perpendicular to s 0 , rotated clockwise by 90 • around s 0 . Using β 2 = 2 β 2 in Eq. (15), we also obtain It is worth to specialize the above equations to the case in which strain is caused by hydrostatic pressure. Let P be the variation of hydrostatic pressure P of the medium in which the cable is immersed, which is transmitted to the fiber by the cable structure and directly by the medium surrounding the fiber, usually petroleum jelly. In this case, = P P , where P = ∂ /∂ P . Equation (16) becomes Equation (18) shows that the spectrum of s reproduces the spectrum of the modulation of the hydrostatic pressure acting upon the fiber. Equivalently, the strain sensed by the polarization averaged wave vector is 0 = 0,P P , with 0,P = ∂ 0 /∂ P . The equality 0 = valid in SSMFs with randomly coupled birefringence of course implies that 0,P = P . Equation (17) also becomes Let us now use Eq. (16) to investigate the temporal correlations of s . In that equation, β ⊥ (z) is z dependent, static, and random, while (t) is deterministic, time dependent, and weakly dependent on z. The birefringence deviations β(z, t) = (t) β(z) at a given position z and times t and t differ only for their lengths and for being parallel or antiparallel to one another. Therefore, the lefthand side of Eq. (8) can be readily generalized to yield β(z, t) · β(z , t ) = g (z − z ) (t) (t ) β 2 , where g (z − z ) is the spatial correlation function defined by Eq. (7). Using this expression, and again the delta function approximation for g (z − z ), we may derive the following equation for the (non-stationary) temporal correlation function of SOP deviations, which generalizes Eq. (17): If strain is caused by pressure, the replacement (t) = P P (t) in Eq. (20) leads to the generalization of Eq. (19). This equation shows that the dynamics of the SOP, and hence its spectrum, does not contain information only on the square modulus of but also on its sign. 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, β = 2 α × β, where α 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, α transforms linear birefringence into a rotated linear birefringence; hence it is orthogonal to β, so that β = 2αβ. In the rotating frame, angles are preserved so that, assuming α deterministic and slowly varying with z, we have β 2 = 4α 2 β 2 , and hence Eq. (15) becomes and Eq. (20) where α(t) is positive if the angle between α and β is π/2 and negative if it is −π/2. A small amount of twist is sufficient to give significant polarization fluctuations, and hence twist may be the dominant mechanism of SOP modulation in terrestrial systems where strong perturbations of anthropic origin act upon the fiber. Different from strain and pressure, this mechanism does not affect the polarization averaged phase because SOP modulation is caused by a rotation of the local birefringent vector. We will show in the following that pressure induced strain is the most likely mechanism for the SOP modulation generated by sea waves in the cable in [4]. In addition, the mechanisms that may apply twist to the cable are several, and may be related to cable configurations on the seafloor that are difficult to be experimentally distinguished. For instance, significant twist may be generated where the cable is unintentionally suspended over the seafloor because of rough bathymetry, but in this case, it is likely that the response to external perturbations is so strong that the SOP fluctuations would cover the whole Poincaré sphere, and this circumstance never occurred in the events recorded during the measurement campaign reported in [4]. For all these reasons, we will not consider twist further in this paper, and concentrate instead on strain and pressure induced SOP modulations.

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 β) is perturbed, to first order, the output polarization changes by where β ⊥ is the component of the birefringence perturbation β 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 ϕ = z 0 β 0 dz = (2π/λ) z 0 n dz . When there is a strong spatial inhomogeneity of with positive and negative regions, then ϕ 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 ϕ 2 = (2π/λ) 2 z 0 z 0 n(z ) (z )n(z ) (z )dz dz , the values outside the main diagonal average, out and the only surviving terms will be those with z z . 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.

COMPARISON WITH OBSERVATIONS
Polarization fluctuations generated by ocean waves and earthquakes have been recently detected in an L = 10,500 km long submarine cable from Los Angeles to Valparaiso, the Curie cable [4]. The PMD coefficient of the cable ranges from 0.01 ps/ √ km and 0.07 ps/ √ km, with an average of 0.03 ps/ √ 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 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 doublefrequency secondary microseism, which is the seismic waves produced by wave-wave interactions, is not observed on the spectrum of 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 s caused by the ocean swells in the Curie cable is of the order of 10 −3 Hz −1 . This value, once integrated over a bandwidth of about 0.01 Hz, gives s 2 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 P = P 0 / cosh(k w h), where k w = 2π/λ w is the wave vector of the ocean wave. The depth of the first 12 km of the Curie cable is about h = 100 m. The frequency of the wave is ω = g k w tanh(k w h) so that, being ω = 2π0.06 rad/s (see Fig. 1), we have k w = 0.0158 m −1 , corresponding to λ w = 2π/k w = 405 m. Assuming h = 100 m and P 0 = 20 kPa (corresponding to 2 m high sea waves), we obtain P = P 0 / cosh(k w h) = 8 kPa. Using in Eq. (19) the observed value for the average of polarization modulation s 2 = 10 −5 , for the PMD coefficient the value κ = 0.03 ps/ √ 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 P the estimate P = 3.5 · 10 −9 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 Hz to 10 Hz, with a coefficient of proportionality (calculated from the peak values |ε 0 ( f max )| 2 = 10 −2.2 µ 2 /Hz and | P ( f max )| 2 = 10 14.7 µPa 2 /Hz) ε 0,P = ε 0 / P 3.6 · 10 −9 Pa −1 . The optical strain is obtained through multiplication by the photoelastic scaling factor ξ = 0.78 [12], giving 0,P 2.8 · 10 −9 Pa −1 . This value is equal within the experimental uncertainties to that characterizing the sensitivity of polarization to pressure of the Curie cable, P = 3.5 · 10 −9 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 P is similar to the measured one of 0,P , although in a different cable, validates the equality, theoretically predicted in Section 3, 0 = . 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) = 0 , where 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 s , i.e., s 1 and 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 s 2 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 κ = 0.03 ps/ √ km, the observed value of s 2 = 0.01 used in Eq. (17) requires = 3.1 · 10 −4 . If the source of perturbation is the direct action of pressure on the fiber, setting = P P = 3.1 · 10 −4 and using the value Fig. 1. Spectrogram of s (sum of the power spectral density of its two components) obtained from 1 Jun 2020 to 12 Jul 2020 in the Curie cable [4]. The figure shows dispersive wavepackets around 0.06 Hz, in the primary microseism band, from ocean swells, each lasting for a few days.  [4]. The red, blue, and green dashed lines are earthquake origin time and predicted P and S arrival times, respectively, from the earthquake to the closest point along the Curie cable. The traces are filtered from 0.5 to 2.5 Hz. P = 3.5 · 10 −9 Pa −1 returns for the amplitude of the pressure modulation P 88 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.

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 s 1 and s 2 , the components of 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 s has the same bandwidth as (t) integrated over the whole link and weighted by random fiber birefringence. If the spatial dependence of (t) were a rectangular function of z, i.e., constant over the section interested by the perturbation and zero elsewhere, then s would be a randomly oriented two-dimensional vector proportional to (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 s 1 and 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 s that are different and cover an entire, slightly asymmetrical, region of the Poincaré sphere around the unperturbed Stokes vector. This is because (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 s 2 is slightly smaller than that of s 1 and then tends to equalize on average. Most earthquake detections observed in the campaign in [4] show traces of s 1 and s 2 that have similar but nonidentical amplitudes. This is because (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 (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 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 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 s , which therefore reproduces quite faithfully the spectrum of (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].

SENSITIVITY
In Fig. 1, we resolve deviations s from the unperturbed output polarization s 2 10 −5 , which correspond to an angular deviation of s over the Poincaré sphere by θ 3 · 10 −3 rad. A full rotation of the Poincaré sphere is equivalent to a phase-shift between two eigenpolarizations of 2π rad, which is in turn generated by an optical path length difference of one wavelength, λ = 1.55 · 10 −6 m. Our virtual interferometric measurement is therefore sensitive to an optical path length difference of (δn ) = λθ/(2π ) = 8 · 10 −10 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 = 1.5 · 10 7 m for the Curie cable, the sensitivity is about (δn )/(n ) = 5 · 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.

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.