P. L. Makowski and A. W. Domanski
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
1825
Approach for fast numerical propagation of uniformly polarized random electromagnetic fields in dispersive linearly birefringent systems Piotr L. Makowski* and Andrzej W. Domanski Faculty of Physics, Warsaw University of Technology, Koszykowa 75, 00-662 Warsaw, Poland *Corresponding author: [email protected] Received May 8, 2013; revised August 1, 2013; accepted August 1, 2013; posted August 2, 2013 (Doc. ID 190016); published August 22, 2013 An efficient simulation technique is proposed for computing propagation of uniformly polarized statistically stationary fields in linear nonimage-forming systems that includes dispersion of linear birefringence to all orders. The method is based on the discrete-time Fourier transformation of modified frequency profiles of the spectral Stokes parameters. It works under the condition that all (linearly) birefringent sections present in the system are described by the same phase birefringence dispersion curve, being a monotonic function of the optical frequency within the bandwidth of the light. We demonstrate the technique as a supplement for the Mueller–Stokes matrix formalism extended to any uniformly polarized polychromatic illumination. Accuracy of its numerical implementation has been verified by using parameters of a Lyot depolarizer made of a highly birefringent and dispersive monomode photonic crystal fiber. © 2013 Optical Society of America OCIS codes: (030.0030) Coherence and statistical optics; (260.1440) Birefringence; (260.2030) Dispersion; (260.5430) Polarization; (060.2420) Fibers, polarization-maintaining. http://dx.doi.org/10.1364/JOSAA.30.001825
1. INTRODUCTION In order to calculate the polarization state of a nonmonochromatic random electromagnetic beam passing through any system, one has to perform integration over spectral components of the field. In the case of one-dimensional propagation in systems with linear birefringence it can be done in several different ways. In the literature concerning description of the Lyot depolarizer we find the Mueller–Stokes matrix equations integrated over the power spectrum density (PSD) of the passing light [1,2] and depolarization formulas based on the complex degree of coherence being the Fourier transform of the PSD [3–5]. In [6] the problem of the optimal placement of a Lyot depolarizer within the fiber-optic gyroscope has been discussed with the use of the spectral coherency matrix formalism [7]. In our previous paper [8] we show how to approximate integrals of this kind by a summation [discretetime Fourier transform (DTFT)] in a way that accounts for an arbitrary dispersion dependence of linear birefringence in the system. That formulation, however, is restricted to fields originating from a single optical disturbance (as are those discussed in [1–5]) and thus always fully polarized in the spectral domain [9]. In the present paper this approach is extended for collimated beams or modal fields in single-mode optical fibers with any frequency distribution of spectral degree of polarization [9] as long as they are uniformly polarized. We also discuss practical aspects of numerical implementation of the concept. The presented technique, based on a frequencyintegrated variant of the Mueller–Stokes matrix formalism, can be useful in numerical modeling of a class of bulk or fiber-optic devices, such as polarimetric sensors [10,11]. It does not cover the general case of a partially coherent and 1084-7529/13/091825-07$15.00/0
partially polarized stationary beam in a uniaxial medium [12], for which the numerical integration of the (two-point) crossspectral density (CSD) matrix over frequency has to be performed in a different way. It is also not suitable for modeling light in inhomogeneous birefringent media other than optical systems formed by a series of homogeneous sections. In case of a well collimated field propagating in an anisotropic system with continuous or stochastic spatial variation of the material parameters the differential version of the Jones calculus or the Mueller–Stokes calculus should be used [13–15].
2. PHYSICAL PRINCIPLE In this section we derive the polarization (coherence) matrix of any uniformly polarized field transmitted through a single piece of uniaxial medium in a way to account for its phase birefringence dispersion curve. This will provide physical grounds for a numerical approach that can be extended for multisection systems with linear birefringence. The reasoning conducted here applies both to properly cut anisotropic crystals (for collimated beams in the paraxial region) and to weakly guiding nontwisted single-mode polarization-maintaining fibers. The second-order correlations of the electric vector of a statistically stationary (in the wide sense) beam-like field propagating in the paraxial region in the positive z direction is completely described by the CSD matrix Wρ1 ; ρ2 ; z; ω [16], where ρ x; y indicates a point in a given z constant plane. If the field is uniformly polarized (in the spacefrequency domain) throughout the region A in the z 0 plane and its transverse intensity profile is essentially kept unchanged along certain propagation distance z L, the CSD matrix, evaluated in coinciding points ρ1 ρ2 ρ, © 2013 Optical Society of America
1826
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
Wρ; ρ; ω hE i ρ; ωEj ρ; ωi;
i; j x; y;
P. L. Makowski and A. W. Domanski
(1)
is sufficient to describe propagation along any line parallel to the z axis bounded by the region A and the propagation distance z L. The angle brackets in Eq. (1) denote the ensemble average over all possible realizations of the complex amplitudes Ex ρ; ω, Ey ρ; ω of the random but strictly monochromatic waves that constitute the polychromatic beam, according to the coherent-mode (scalar-mode) expansion of the CSD matrix [17–19]. The matrix Wρ; ρ; ω Wω may be also called the spectral polarization (coherence) matrix of the field because of its strong resemblance to the spectral coherency matrix defined by Barakat [7]. In order to combine all the frequencies present in the beam at z 0 we switch to the complex mutual coherence (covariance) matrix using the generalized Wiener–Khintchine theorem [20]: Z Γτ
∞
−∞
Wω exp−iωτdω;
(2)
which is also defined by Γτ hE i tE j t τi;
i; j x; y:
(3)
Here, the correlated random functions Ex ρ; t, Ey ρ; t are the complex analytic signals associated with the Cartesian components of the fluctuating electric vector. For the time delay τ 0 the covariance matrix is reduced to the Wolf’s polarization (coherence) matrix: Γ0 J. We notice that transformation (2) can describe onedimensional propagation of an uniformly polarized field along a single section of a homogeneous lossless uniaxial medium with its optical axis perpendicular to the z axis. It follows directly from the physical model and from definition (3) that the field transmitted through such medium of length z is described by the polarization matrix of the form [9] Jz
Γxx 0 Γxy τz ; Γxy τz Γyy 0
(4)
with τz z∕cB being the physical time delay introduced between the birefringence axes of the medium; c is the speed of light in vacuum, and B is the phase birefringence. It follows from the conventions used in Eqs. (3) and (4) that one gets τz > 0 when y is the fast axis. This implies that the phase birefringence in τz z∕cB has to be defined eff as B neff x − ny . The result, Eq. (4), is valid only under the assumption that the delay τz is a constant for all spectral components of the field. In order to account for the birefringence dispersion we can simply acknowledge the frequency dependence of τ τz in Eq. (2), τz ω z∕cBω, which makes z the constant parameter of the integral: Z Γij z
∞
−∞
W ij ω exp−iωτz ωdω:
(5)
The resulting expression (5) is meant as a replacement for Eq. (2) in Eq. (4) and can raise doubts, as it is clearly no longer a Fourier integral and thus it is not supported by the Wiener– Khintchine theorem. Definition (3) also loses its meaning for
τ τω. There is a physical justification for Eq. (5), however, which we will provide in Section 4. For now it is of interest to notice that Eq. (5) can be reduced back to the Fourier integral form by applying a variable substitution: ω ~ Fω
Bω ω Bω0
(6)
with an arbitrarily chosen effective frequency ω0 . By applying Eq. (6) to Eq. (5) we define a new quantity, which can be called the effective complex mutual coherence function: Γij τz 0
Z
∞
−∞
~ ij ω ~ exp−iωτ ~ z ~ W 0 dω;
−1 dF ω ~ ~ ij ω ~ ~ W F −1 ω; W ~ ij dω
(7a)
(7b)
where τz 0 z∕cBω0 is the effective time delay between the polarization components in the medium at distance z. ~ indicates the inverse of the transformation function deF −1 ω ~ ~ fined by Eq. (6). Both τz 0 and the transformed profiles W ij ω depend on the particular choice of ω0 in a way that makes the overall result independent of ω0 . Transformation formula (7b) is independent of the propagation distance. Representations (7) make sense only if Fω is invertible in the range of optical frequencies where jW ij ωj > 0, i; j x; y. If this condition is met, relation (7a), being a Fourier integral again, carries the same information as Eq. (5). The physical meaning of this correspondence can be expressed as follows: In the case of one-dimensional propagation of a well collimated uniformly polarized field, the influence of the birefringence dispersion on the polarization state of light passing through a linearly birefringent medium can be reproduced by means of a frequency-domain transformation of the spectral polarization (coherence) matrix of the illuminating field. In other words, in place of a medium of length z and birefringence Bω and illuminated by a field described by the matrix Wω, one can take a medium of the same length z with constant birefringence Bω0 and illuminated by light ~ ω, ~ and the state of light at described by a different matrix W the output would be the same. This interpretation, however, has limited use. It follows from Eq. (6) that, when Bω switches its sign within the passband of the field, the domain ~ ij ω ~ will naturally include negof the transformed profiles W ~ ij ω ative frequencies. Therefore, in this case the profiles W ~ can no longer be regarded as standard CSD distributions, which implies that the discussed correspondence is no longer physically realizable. Consequently, the correlation function at the left-hand site of Eq. (7a) is not an analytic signal in general; thus it loses one of the basic properties of the ordinary complex mutual coherence function ([20], part 3.1.3). Nevertheless, it is still a valid expression of correlations in the vectorial field propagating through a medium of the considered type. This is ensured by the equivalence of the integrals of Eqs. (7a) and (5). In Section 5 we illustrate usefulness of the form (7a) in numerical simulations in cases where the negative frequencies occur. The effective complex mutual coherence functions (7a) substituted into Eq. (4) describe propagation in a single birefringent and dispersive element. It follows from the formulas
P. L. Makowski and A. W. Domanski
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
derived in [3–5] that in case of a two-section system, such as the Lyot depolarizer, we need mutual coherence functions of the form Γxy τz2 τz1 . Taking into account the birefringence dispersion, we can express them as Γxy τz2 ω τz1 ω Z ∞ z z W ij ω exp −iω 2 B2 ω 1 B1 ω dω: c c −∞
(8)
It is apparent from Eq. (8) that in order to enforce the linear frequency dependence of the term in square brackets by means of a variable substitution, as it was done for Eq. (5), it is necessary to have B1 ω B2 ω Bω for all ω. In conclusion, in general case of multisection systems the direct usage of the representation in Eqs. (6) and (7) of the effective mutual coherence functions is limited by the two conditions: (a) All the (linearly) birefringent sections have to be described by the same phase birefringence curve Bω and (b) The difference of effective propagation constants Δβω ω∕cBω being proportional to Fω defined by Eq. (6) must not possess local extrema within the bandwidth of the passing light, because Eq. (7b) requires the existence of the inverse of Fω. In cases where condition (b) is not met, an indirect use of representations (7) is still possible. It is done by splitting the initial integral of Eq. (5) into a sum of integrals with neighboring subdomains in which Fω is invertible separately. This sum can be transformed into a sum of integrals of the form ~ intervals (overlapping in general). of Eq. (7a) with finite ω Because the local extremal points of Fω will correspond to singularities in the transformed cross-spectral densities, ~ intervals), one can preEq. (7b) (at the related edges of the ω dict that every local extremum in Fω will cause an oscillatory component in the effective mutual coherence function Γij τz 0 ([20], part 2.4.2). It can be easily verified numerically by calculating the integral of Eq. (5) that a local extremum in Bω results in a correlation function whose absolute value does not approach zero for increasing medium length z.
3. NUMERICAL IMPLEMENTATION As the form of Eq. (7a) represents a Fourier integral again, we can easily approximate it for fast numerical computation. Since the real and imaginary parts of the elements Γij τ essential for determining the state of polarization of the beam often exhibit rapid oscillations within the range of the coherence time of the light wave, the DTFT is a convenient solution: ~ Γij τz 0 ≈ jΔωj
N X
~ ij ω ~ n exp−iω ~ n τz W 0 ;
(9)
n1
where N is the total number of equidistant sampling points ~ n g with the sampling step Δω ~ ω ~N −ω ~ 1 ∕N − 1. The fω ~ n g should be chosen in such a way to make range of fω ~ 1 − Δω∕2; ~ ~ FDω ω ω~ N Δω∕2, where Fω is defined by Eq. (6) and Dω fω∶jW ij ωj > 0g. This ensures full accuracy of expression (9) (exact equality) if the following conditions are strictly fulfilled: (a) we consider τz 0 ∈ −T 0 ∕2; T 0 ∕2, ~ is the period of the transform, and where T 0 2π∕jΔωj (b) the value T 0 is large enough (Δω ~ is set small enough)
1827
~ ij ω ~ n g to eliminate overlapping of the for given elements fW replicas of Γij τz . Because of the continuity in the transform 0 domain there is no periodicity in the signal domain, which implies that the sum of Eq. (9) is actually performed over the ~ but only infinite set of samples covering the entire domain ω, ~ n g contribute to the sum. This fact elimthose within the set fω inates spectral leakage errors. It is formally required that one ~ 0, which restricts the of the samples occupy the point ω choice of ω1 , ωN , and N, but violation of this rule does not affect the shape of the transform within the range of interest: −T 0 ∕2; T 0 ∕2. In order to utilize Eq. (9) in calculation of Jz one has to ~ ij ω ~ n first. obtain the transformed complex distributions W Using the DTFT approximation directly on Eq. (5) in the presence of the birefringence dispersion leads to errors discussed in our previous work [8] for the special case of light linearly polarized at input. These errors originate from the variable distance between the subsequent phases inside the term, exp−iωn τz ωn in Eq. (5), which is inconsistent with the definition of the DTFT. This is why transformation [Eq. (7b)] of the CSD profiles is needed before discretization on a regular grid. ~ ij ω ~ n through direct implementation of Obtaining W Eqs. (6) and (7b) is highly inefficient. However, both evalu~ −1 ω ~ and numerical inversion ation of the derivative d∕dωF of Fω itself can be avoided with little loss of precision if the following transformation procedure is applied. Suppose that Ωω denotes the real or the imaginary part of the original CSD distribution W ij ω. We define a set of samples ωm , where m 1; 2; …; M. If the points fωm g are equidistant, the values fΩωm g can be treated as “frequencies” of a histogram approximating the continuous distribution Ωω, with boundaries of the bins positioned at ωm − Δω∕2 and ωm Δω∕2 for each sample m 1; 2; …; M. Next, we ~ domain and associate center map these boundaries into the ω ~ 0m g, with new bins of a points of the obtained intervals, say fω transformed histogram whose frequencies are of the form ~ 0 ω ~ 0m κm Ωωm g. The coefficients κm Δω∕jFωm fΩ Δω∕2 − Fωm − Δω∕2j ensure an unchanged total area ~ 0 ω ~ 0m g are therefore covered by the bins. The new values fΩ (irregularly distributed) samples of the transformed density ~ ω. ~ They should be resampled on a regular distribution Ω grid with the desired number of points N. The resulting ~ ω ~ n g, n 1; 2; …; N, obtained for the real two data sets fΩ and the imaginary part of W ij ω determine the discrete com~ ij ω ~ n for formula (9). plex distribution W Although the described procedure does not seem to require Fω to be invertible, this restriction still holds. A local ~ ω extremum of Fω implies a singularity in Ω ~ (visible as κm → ∞), which is impossible to represent accurately by a finite number of samples (in the sense of the Nyquist–Shannon sampling theorem). The algorithm would also generate over~ 0 ω ~ 0m g. lapping bins in the new histogram fΩ The procedure is not exact due to the assumption that each sampling point ω ~ 0m corresponds to the central point between the boundaries ωm − Δω∕2 and ωm Δω∕2 mapped to the ~ domain. This is strictly correct only if Fω can be approxiω mated by a linear function within every region covered by every single bin of the histogram fΩωm g. Equation (9) supplemented by the above transformation procedure constitutes the basis of the proposed computational
1828
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
P. L. Makowski and A. W. Domanski
technique. It includes the linear birefringence dispersion dependence into the (effective) mutual coherence functions of the uniformly polarized field passing through a system with one or many linearly birefringent elements meeting the conditions listed in the previous section. These effective mutual coherence functions can be employed directly as replacements for the corresponding terms in some known descriptions based on the complex degrees of coherence [3–5,21]. Those solutions, however, are not in principle dedicated to numerical modeling. It is then of interest to ask how the presented technique relates to the Mueller–Stokes matrix formalism in polychromatic illumination [1,2,5,8], which is equivalent to the spectral coherency matrix formalism [6,7].
4. APPLICATION IN MUELLER–STOKES CALCULUS Every polychromatic field characterized by the frequency profile of the spectral polarization matrix Wω can be equivalently described by the spectral Stokes vector sω [22]. The ordinary (temporal) Stokes vector S of the field is obtained through the integral Z S
∞ −∞
sωdω;
3 3 2 J xx J yy Z ∞ W xx ω W yy ω 6 W xx ω − W yy ω 7 6 J xx − J yy 7 7 7 6 6 4 2ReJ xy 5 −∞ 4 2ReW xy ω 5dω; 2ImJ xy 2ImW xy ω
(10a)
2
(10b)
which reflects the generalized Wiener–Khintchine theorem (5) R evaluated for τ 0: J Wωdω. On the other hand, it follows from representation (1) of the matrix Wω that each spectral Stokes vector sω can be regarded as a second-order statistic of an ensemble of strictly monochromatic waves: fEρ; ω E x ρ; ω; Ey ρ; ωg. It has been confirmed experimentally [23] that the complex amplitudes of these waves are measurable and that they propagate in accordance with Maxwell’s equations. Therefore, it is justified to model propagation of each realization, 2 3 E x ωEx ω E y ωEy ω 6 E x ωEx ω − E y ωEy ω 7 7; sω 6 (11) 4 5 2ReE x ωE y ω 2ImEx ωEy ω such that hsωi sω by use of the Mueller–Stokes matrix formalism: s
out
ω Mωs ω: in
(12)
If we now take the ensemble average of both sides in Eq. (12), we get sout ω
hMωsin ωi
Mωsin ω;
(13)
where Mω is a Mueller matrix of a static medium or an averaged Mueller matrix of a stochastic medium that fluctuates independently from the field. Relation (13) combined with (10a) defines the Mueller–Stokes matrix formalism for any polychromatic illumination uniformly polarized in the space-frequency domain:
Z Sout
∞
−∞
Meff ωsin ωdω;
(14)
where Meff is the averaged effective Mueller matrix of a linear system. Because the resulting vector Sout is meant to describe the whole cross section of the field, we require that it be uniformly polarized in the space–time domain as well. The spectral Stokes vector sin ω represents monochromatic components whose spectral degrees of polarization Pω [9] can take any value between zero and unity. It can be shown that the Mueller–Stokes equations integrated over the optical frequency that we find in works by Billings [1] and Loeber [2] in the context of the Lyot depolarizer are special cases of Eq. (14). This correspondence takes place under the assumption that sin represents the Stokes vector common for all spectral components, which are also spectrally fully polarized and thus expressible by deterministic plane waves. If Meff is additionally derivable from a Jones matrix (a Mueller–Jones matrix [24]), Eq. (14) becomes equivalent to the coherency matrix formalism defined by Barakat [7]. In this technique, the polarization state of the field is represented directly by the coherency matrix Φω, which is essentially the same as Wω defined by the inverse form of Eq. (2), apart from the symmetry of the frequency profiles of the elements; the optical system is represented by the effective Jones matrix, which can express a single train of optical elements or a sum of two matrices corresponding to the two round-trip paths in an interferometer [6]. In order to switch from this technique to the formalism defined by Eq. (14) one has to replace the matrix Win ω by the equivalent spectral Stokes vector sin ω and transform the effective Jones matrix of the system into the corresponding Mueller– Jones matrix [24]. Because of the use of real-valued quantities only and because of the simple algebra of the polarization state transformation, the method of Eq. (14) is well suited for numerical modeling. As for the theoretical analysis, the formalism by Barakat [7] and the Jones formalism are more adequate [25]. Let us now consider propagation through a linearly birefringent element, as discussed in Section 1. We will prove that the output Stokes vector obtained by Eq. (14) is also expressed by elements of polarization matrix (4) for τz τz ω. This will provide physical justification for Eqs. (5)–(7) and allow us to approximate Eq. (14) on the same terms that apply to the DTFT form [Eq. (9)] of the coherence functions Γij τ0z . The Mueller matrix of an uniaxial medium with its optical axis perpendicular to the propagation z axis has the form 2 6 Mδ 6 4
1
3 1
7 7; cosδ sinδ 5 − sinδ cosδ
(15)
where δω Δβωz ωτz ω is the retardation for frequency ω. The placement of the minus sign in Eq. (15) is determined by conventions used in definitions of τz , Γτ, and sω and by the proportionality E ω t ∝ exp−iωt. The magnitude of δω is unlimited, as the spectral degree of polarization [16] of the monochromatic wave sω is unaffected by time shifts introduced between its Cartesian components [9]. We also notice that matrix (15) includes birefringence dispersion to all orders in a natural way, but
P. L. Makowski and A. W. Domanski
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
for the moment let us assume that the birefringence in the medium is nondispersive, meaning τz constant for all ω. If we substitute Meff ω Mδω into Eq. (14), then according to definitions (10) the two nontrivial output Stokes parameters take the forms Z S out 2 2
∞
−∞ ∞
Z
−2
−∞
Z S out 3 2
∞
ImW xy ω sin−ωτz dω;
(16a)
ReW xy ω sin−ωτz dω
−∞ ∞
Z
2
ReW xy ω cos−ωτz dω
−∞
ImW xy ω cos−ωτz dω:
(16b)
On the other hand, it follows from Eq. (2) that Γxy τz
Z
∞
−∞
ReW xy ω iImW xy ω
× cos−ωτz i sin−ωτz dω:
(17)
Comparing Eqs. (16) with Eq. (17) leads to the immediate result 3 2 Γxx 0 Γyy 0 6 Γxx 0 − Γyy 0 7 7 (18) Sout z 6 4 2ReΓxy τz 5; z 2ImΓxy τ which is equivalent to Eq. (4), as expected. In the case of a out multisection system we would obtain a vector P S composed of linear combinations of functions Γxy k ak τk with ak 0; 1 and τk τzk , where k numbers the birefringent sections. It is important to note that the established correspondence is satisfied both for the nondispersive τz constant and for the dispersive case, τz τz ω, where the latter case entails replacement of Γxy τz by Γxy z defined through Eq. (5). This is because for τz ω z∕cBω the integrals in Eqs. (16) can be collected to the form of the real and imaginary parts of the right-hand site of Eq. (5). Two essential conclusions follow from this fact. First, definition (5) of the modified coherence function Γxy z is given clear physical sense, as it leads to the same results as the physically well justified Mueller–Stokes equation (14). Second, because the spectral vector sω is composed of the functions W ij ω, the numerical approach for propagation of polychromatic light in an uniaxial medium expressed by Eq. (9) translates into ~ Sout ≈ jΔωj
N X
~ n τz ~ n ; Mω sin ω 0 ~
(19)
n1
which, unlike Eq. (9), can be easily extended for the case of multisection systems composed of linearly birefringent parts and possibly other nonbirefringent elements: ~ Sout ≈ jΔωj
N X n1
z2 zK in ~ n ; τz1 Meff ω ~s ω~ n ; 0 ; τ0 ; …; τ 0
(20)
1829
where k 1; …; K numbers the linearly birefringent elements in the optical system. Despite the fact that Eq. (20) does not resemble any kind of discrete Fourier transformation, it is subject to the same conditions as Eqs. (9) and (7). We keep ~ defined through Eqs. (7), will not always in mind that s~ in ω, represent the standard spectral Stokes vector, as formula (7b) can introduce negative frequencies to the transformed profiles. In the following section we examine the accuracy of expression (20) supplemented by the transformation algorithm described in Section 3. The simulation errors will be evaluated by comparing the output (temporal) degree of polarization predicted by Eq. (20) with the one obtained through direct numerical integration using Eq. (14).
5. NUMERICAL EXAMPLE For the sake of clarity we choose a simple two-section system (the fiber-optic Lyot depolarizer [26]) and uniformly coherent light [9] with unit spectral degree of polarization Pω for all ω, in particular a linearly polarized light wave described completely by its power spectral density distribution Ψω W xx ω W yy ω and the initial azimuth angle α. In this case Eq. (20) takes the form Sout τ
0
~ ≈ jΔωj
X N
~ ~ n ; τ0 ; 2τ0 Ψω ~ n sˆ in ; Meff ω
(21)
n1
~ ω where Ψ ~ is the original PSD Ψω after the domain transformation defined by Eqs. (6) and (7), sˆ in sin α∕sin 0 α is the normalized spectral Stokes vector of the input field linearly ~ n ; τ0 ; 2τ0 polarized with the azimuth angle α, and Meff ω ~ n τ0 is the effective Mueller matrix of M2ω~ n τ0 R45°Mω the depolarizer. In order to obtain the result in the coordinate system associated with the first section, the second matrix has to be surrounded by two rotation matrices: ~ n τ0 R45°. The set of sampling points fω~ n g R−45°M2ω meets analogical conditions as discussed for Eq. (9). The simulations were performed using the PSD of the superluminescent Broadlighter diode Q1350-HP from Superlum, Ltd. with spectral width of about 200 nm. The two sections of the depolarizer are made of a highly birefringent holey photonic crystal fiber [27], single-mode within the spectrum width of the Q1350-HP diode. Its phase birefringence, depicted in Fig. 1, switches sign within the spectrum range of the diode. Figure 2 depicts the result of the transformation defined by Eq. (7b) and the procedure described in Section 3. The original power spectrum samples fΨωm g, m 1; …; M (black points), are replaced by samples of the transformed profile ~ ω ~ n g, n 1; …; N (gray points), which encodes both the fΨ diode spectrum and the phase birefringence curve of the fiber. At the same time, the original curve Bω (black line) is replaced by a constant effective birefringence Bω0 (gray line). The effective frequency has been chosen arbitrarily at the right border of the Q1350-HP diode spectrum, ω0 ωM ; hence Bω0 takes the highest value of Bω. The number of the transformed samples N has been chosen in such a way to make Γ11 τz 0 fade almost completely for τ 0 < T 0 ∕2, where T 0 is the period of DTFT transforms (10). Because the phase birefringence Bω switches sign within the passband of Ψω, ~ ω the transformed samples fΨ ~ n g occupy both positive and
1830
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
Fig. 1. Input data: phase and group birefringence of the highly birefringent holey photonic crystal fiber plotted in the region covered by the Q1350-HP diode spectrum. The phase birefringence Bλ is represented by a fifth-degree polynomial fitting to the samples of the original phase birefringence curve [27].
Fig. 2. Input data: PSD samples of the superluminescent diode source along with associated phase birefringence curve of the highly birefringent photonic crystal fiber (black) and the corresponding functions after domain transformation (7b) (gray).
negative frequencies. This is an expected property of Eqs. (6) ~ simply switches the phase and (7b); the negative frequency ω sign inside the exponential term in (7a), just like the negative sign of the phase birefringence Bω switches the phase sign inside integral (5). The two integrals have to be equivalent, hence the necessity for the negative frequencies in the domain ~ ω ~ in this example. of Ψ Figure 3(a) presents the residual degree of polarization P 1∕S 0 S 21 S 22 S 23 1∕2 [obtained from Eq. (21)] of light emerging from the fiber-optic depolarizer of total length L L1 L2 , where L2 2L1 . The horizontal axis α represents the azimuth angle of the input linear polarization. The white arrows indicate the total length L L1 L2 such that τL2 T 0 ∕2; therefore they determines the end of the 0
P. L. Makowski and A. W. Domanski
expected range of accurate results (the analogue of the Nyquist frequency). In the upper region of the diagram in Fig. 3(a) we observe a nonphysical compensation of depolarization, which corresponds to τL2 T 0 ; it is then an expected 0 manifestation of the periodicity of DTFT (9), which is implicitly present in expression (21), according to the correspondence established in Section 4. Figures 3(b) and 3(c) express the inaccuracy of formula (21) for M 30 and M 50, while the number of transformed samples N 30 remains unchanged (along with the value of the Nyquist delay). The reference degree of polarization values were obtained by direct numerical evaluation of integral (14) corresponding to the sum of Eq. (21) by use of the quad_qags procedure in Maxima CAS (v5.24.0, GCL 2.6.7). ~ ω ~ n g for The discrete spectra fΨωm g transformed to fΨ proposed solution (21) and the corresponding curves Ψω for reference method (14) were prepared in such way to allow evaluation of errors originating from the approximate character of the domain transformation procedure described in Section 3. First, M 50 regularly distributed samples were chosen from the original PSD curve taken from the documentation of the Q1350-HP diode. Then the cubic-spline interpolation of these points was calculated to serve as a simulation of a real-world continuous spectrum passed to reference formula (14) [and used to generate Fig. 3(c)]. Next, a set of M 30 equidistant samples was picked from the interpolation curve, and the second cubic-spline interpolation was calculated from it for Eq. (14) [Fig. 3(b)]. This means that the two diagrams in Figs. 3(b) and 3(c) actually correspond to two separate light sources with slightly different PSDs Ψω, but in return we can be certain that in both cases (M 30 and M 50) the discrete functions Ψωm represent their continuous counterparts Ψω with comparable accuracy. Therefore, if proposed technique (21) turns out to be less accurate for M 30, it is not due to incomplete representation of the details in the PSD profile, but rather due to violated assumptions of the domain transformation procedure. As expected, we observe that the accuracy of the proposed simulation technique remains high as long as τ0L2 < T 0 ∕2. This confirms implicit presence of DTFT in formula (21). On the other hand, the slight discrepancies visible below the T 0 ∕2 line strongly depend on the initial number of the PSD samples: the change from M 30 to M 50 reduces
Fig. 3. Results: (a) degree of polarization of superluminescent diode light after passing a dispersive Lyot depolarizer calculated by the proposed technique, Eq. (21), and its discrepancy from the result of exact numerical integration, Eq. (14), for initial number of the spectrum samples, (b) M 30, and (c) M 50. The white arrows indicate the total length L1 L2 of the depolarizer that corresponds to the half-period T 0 ∕2 of the DTFT approximation of the autocorrelation Γ11 τL2 0 .
P. L. Makowski and A. W. Domanski
the maximal level of jΔDoPj (below the arrow), where DoP is degree of polarization, by a factor of 4, from less than 2 × 10−3 to less than 5 × 10−4 . These values reflect expected limitations of the approximate domain transformation procedure defined in Section 3 for the discretized CSD functions. The estimated absolute integration errors reported by the quad_qags algorithm did not exceed the level of 1 × 10−8 ; therefore the inaccuracies of the reference method can be ignored. In the more general case of the field described by the matrix Wω the domain transformation has to be performed on all four elements of the matrix. In the case of the complex-valued element W 12 ω W 21 ω it is performed for the real and the imaginary part separately.
6. CONCLUSION We have discussed conditions under which an uniformly polarized polychromatic field can be numerically propagated in dispersive birefringent systems by the use of the DTFT of the discretized frequency profile of the CSD matrix Wρ; ρ; ωn . The proposed numerical approach utilizes a form of the frequency-integrated Mueller–Stokes equation that involves spectral Stokes vectors. The dispersion dependence of birefringence of the (homogeneous) linearly birefringent elements in the system is included to all orders by means of a domain transformation of the frequency profiles of these vectors at the input plane. The sources of errors in the proposed approximation have been identified, and the expected accuracy range has been numerically confirmed by using parameters of a photonic crystal fiber Lyot depolarizer and wideband superluminescent diode source. The technique can be useful in modeling and design of some types of bulk optics devices as well as polarimetric and interferometric fiber-optic sensors.
REFERENCES B. H. Billings, “A monochromatic depolarizer,” J. Opt. Soc. Am. 41, 966–968 (1951). 2. A. Loeber, “Depolarization of white light by a birefringent crystal. II. The Lyot depolarizer,” J. Opt. Soc. Am. 72, 650–656 (1982). 3. W. K. Burns, “Degree of polarization in the Lyot depolarizer,” J. Lightwave Technol. 1, 475–479 (1983). 4. K. Mochizuki, “Degree of polarization in jointed fibers: the Lyot depolarizer,” Appl. Opt. 23, 3284–3288 (1984). 5. P. L. Makowski, M. Z. Szymanski, and A. W. Domanski, “Lyot depolarizer in terms of the theory of coherence-description for light of any spectrum,” Appl. Opt. 51, 626–634 (2012). 6. E. I. Alekseev and E. N. Bazarov, “Theoretical basis of the method for reducing drift of the zero level of the output signal of a fiber-optic gyroscope with the aid of a Lyot depolarizer,” Sov. J. Quantum Electron. 22, 834–839 (1992). 7. R. Barakat, “Theory of the coherency matrix for light of arbitrary spectral bandwidth,” J. Opt. Soc. Am. 53, 317–322 (1963).
1.
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
1831
8. P. L. Makowski and A. W. Domanski, “Approach for modeling influence of birefringence dispersion on polarization properties of multi-section systems,” Proc. SPIE 8697, 869703 (2012). 9. P. Réfrégier, T. Setälä, and A. T. Friberg, “Temporal and spectral degrees of polarization of light,” Proc. SPIE 8171, 817102 (2011). 10. T. R. Woliński, “Polarimetric optical fibers and sensors,” in Progress in Optics, E. Wolf, ed. (Elsevier, 2000), Vol. 40, pp. 1–75. 11. T. Nasilowski, F. Berghmans, T. Geernaert, K. Chah, J. Van Erps, G. Statkiewicz, M. Szpulak, J. Olszewski, G. Golojuch, T. Martynkien, W. Urbanczyk, P. Mergo, M. Makara, J. Wojcik, C. Chojetzki, and H. Thienpont, “Sensing with photonic crystal fibres,” in IEEE International Symposium on Intelligent Signal Processing, 2007. WISP 2007 (2007), pp. 1–6. 12. D. J. Liu and Z. X. Zhou, “Propagation of partially polarized, partially coherent beams in uniaxial crystals orthogonal to the optical axis,” Eur. Phys. J. D 54, 95–101 (2009). 13. R. C. Jones, “A new calculus for the treatment of optical systems. VII. Properties of the N-matrices,” J. Opt. Soc. Am. 38, 671–683 (1948). 14. N. Ortega-Quijano and J. L. Arce-Diego, “Generalized Jones matrices for anisotropic media,” Opt. Express 21, 6895–6900 (2013). 15. R. M. A. Azzam, “Propagation of partially polarized light through anisotropic media with or without depolarization: a differential 4×4 matrix calculus,” J. Opt. Soc. Am. 68, 1756–1767 (1978). 16. E. Wolf, “Unified theory of coherence and polarization of random electromagnetic beams,” Phys. Lett. A 312, 263–267 (2003). 17. F. Gori, M. Santarsiero, R. Simon, G. Piquero, R. Borghi, and G. Guattari, “Coherent-mode decomposition of partially polarized, partially coherent sources,” J. Opt. Soc. Am. A 20, 78–84 (2003). 18. J. Tervo, T. Setälä, and A. T. Friberg, “Theory of partially coherent electromagnetic fields in the space-frequency domain,” J. Opt. Soc. Am. A 21, 2205–2215 (2004). 19. M. A. Alonso and E. Wolf, “The cross-spectral density matrix of a planar, electromagnetic stochastic source as a correlation matrix,” Opt. Commun. 281, 2393–2396 (2008). 20. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995). 21. P. L. Makowski and A. W. Domanski, “Extended Mueller–Stokes description of polarization-mode transformation in linearly birefringent single mode optical fibers,” Opt. Lett. 38, 1107–1109 (2013). 22. O. Korotkova and E. Wolf, “Generalized Stokes parameters of random electromagnetic beams,” Opt. Lett. 30, 198–200 (2005). 23. A. Dogariu and G. Popescu, “Measuring the phase of spatially coherent polychromatic fields,” Phys. Rev. Lett. 89, 243902 (2002). 24. C. Brosseau, “Mueller matrix analysis of light depolarization by a linear optical medium,” Opt. Commun. 131, 229–235 (1996). 25. C. Tsao, Optical Fibre Waveguide Analysis (Oxford University, 1992). 26. K. Böhm, K. Petermann, and E. Weidel, “Performance of Lyot depolarizers with birefringent single-mode fibers,” J. Lightwave Technol. 1, 71–74 (1983). 27. L. Zhang, T. Luo, Y. Yue, and A. E. Willner, “High group birefringence in photonic crystal fibers with both positive and negative phase birefringences,” Pure Appl. Opt. 10, 035004 (2008).
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
1825
Approach for fast numerical propagation of uniformly polarized random electromagnetic fields in dispersive linearly birefringent systems Piotr L. Makowski* and Andrzej W. Domanski Faculty of Physics, Warsaw University of Technology, Koszykowa 75, 00-662 Warsaw, Poland *Corresponding author: [email protected] Received May 8, 2013; revised August 1, 2013; accepted August 1, 2013; posted August 2, 2013 (Doc. ID 190016); published August 22, 2013 An efficient simulation technique is proposed for computing propagation of uniformly polarized statistically stationary fields in linear nonimage-forming systems that includes dispersion of linear birefringence to all orders. The method is based on the discrete-time Fourier transformation of modified frequency profiles of the spectral Stokes parameters. It works under the condition that all (linearly) birefringent sections present in the system are described by the same phase birefringence dispersion curve, being a monotonic function of the optical frequency within the bandwidth of the light. We demonstrate the technique as a supplement for the Mueller–Stokes matrix formalism extended to any uniformly polarized polychromatic illumination. Accuracy of its numerical implementation has been verified by using parameters of a Lyot depolarizer made of a highly birefringent and dispersive monomode photonic crystal fiber. © 2013 Optical Society of America OCIS codes: (030.0030) Coherence and statistical optics; (260.1440) Birefringence; (260.2030) Dispersion; (260.5430) Polarization; (060.2420) Fibers, polarization-maintaining. http://dx.doi.org/10.1364/JOSAA.30.001825
1. INTRODUCTION In order to calculate the polarization state of a nonmonochromatic random electromagnetic beam passing through any system, one has to perform integration over spectral components of the field. In the case of one-dimensional propagation in systems with linear birefringence it can be done in several different ways. In the literature concerning description of the Lyot depolarizer we find the Mueller–Stokes matrix equations integrated over the power spectrum density (PSD) of the passing light [1,2] and depolarization formulas based on the complex degree of coherence being the Fourier transform of the PSD [3–5]. In [6] the problem of the optimal placement of a Lyot depolarizer within the fiber-optic gyroscope has been discussed with the use of the spectral coherency matrix formalism [7]. In our previous paper [8] we show how to approximate integrals of this kind by a summation [discretetime Fourier transform (DTFT)] in a way that accounts for an arbitrary dispersion dependence of linear birefringence in the system. That formulation, however, is restricted to fields originating from a single optical disturbance (as are those discussed in [1–5]) and thus always fully polarized in the spectral domain [9]. In the present paper this approach is extended for collimated beams or modal fields in single-mode optical fibers with any frequency distribution of spectral degree of polarization [9] as long as they are uniformly polarized. We also discuss practical aspects of numerical implementation of the concept. The presented technique, based on a frequencyintegrated variant of the Mueller–Stokes matrix formalism, can be useful in numerical modeling of a class of bulk or fiber-optic devices, such as polarimetric sensors [10,11]. It does not cover the general case of a partially coherent and 1084-7529/13/091825-07$15.00/0
partially polarized stationary beam in a uniaxial medium [12], for which the numerical integration of the (two-point) crossspectral density (CSD) matrix over frequency has to be performed in a different way. It is also not suitable for modeling light in inhomogeneous birefringent media other than optical systems formed by a series of homogeneous sections. In case of a well collimated field propagating in an anisotropic system with continuous or stochastic spatial variation of the material parameters the differential version of the Jones calculus or the Mueller–Stokes calculus should be used [13–15].
2. PHYSICAL PRINCIPLE In this section we derive the polarization (coherence) matrix of any uniformly polarized field transmitted through a single piece of uniaxial medium in a way to account for its phase birefringence dispersion curve. This will provide physical grounds for a numerical approach that can be extended for multisection systems with linear birefringence. The reasoning conducted here applies both to properly cut anisotropic crystals (for collimated beams in the paraxial region) and to weakly guiding nontwisted single-mode polarization-maintaining fibers. The second-order correlations of the electric vector of a statistically stationary (in the wide sense) beam-like field propagating in the paraxial region in the positive z direction is completely described by the CSD matrix Wρ1 ; ρ2 ; z; ω [16], where ρ x; y indicates a point in a given z constant plane. If the field is uniformly polarized (in the spacefrequency domain) throughout the region A in the z 0 plane and its transverse intensity profile is essentially kept unchanged along certain propagation distance z L, the CSD matrix, evaluated in coinciding points ρ1 ρ2 ρ, © 2013 Optical Society of America
1826
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
Wρ; ρ; ω hE i ρ; ωEj ρ; ωi;
i; j x; y;
P. L. Makowski and A. W. Domanski
(1)
is sufficient to describe propagation along any line parallel to the z axis bounded by the region A and the propagation distance z L. The angle brackets in Eq. (1) denote the ensemble average over all possible realizations of the complex amplitudes Ex ρ; ω, Ey ρ; ω of the random but strictly monochromatic waves that constitute the polychromatic beam, according to the coherent-mode (scalar-mode) expansion of the CSD matrix [17–19]. The matrix Wρ; ρ; ω Wω may be also called the spectral polarization (coherence) matrix of the field because of its strong resemblance to the spectral coherency matrix defined by Barakat [7]. In order to combine all the frequencies present in the beam at z 0 we switch to the complex mutual coherence (covariance) matrix using the generalized Wiener–Khintchine theorem [20]: Z Γτ
∞
−∞
Wω exp−iωτdω;
(2)
which is also defined by Γτ hE i tE j t τi;
i; j x; y:
(3)
Here, the correlated random functions Ex ρ; t, Ey ρ; t are the complex analytic signals associated with the Cartesian components of the fluctuating electric vector. For the time delay τ 0 the covariance matrix is reduced to the Wolf’s polarization (coherence) matrix: Γ0 J. We notice that transformation (2) can describe onedimensional propagation of an uniformly polarized field along a single section of a homogeneous lossless uniaxial medium with its optical axis perpendicular to the z axis. It follows directly from the physical model and from definition (3) that the field transmitted through such medium of length z is described by the polarization matrix of the form [9] Jz
Γxx 0 Γxy τz ; Γxy τz Γyy 0
(4)
with τz z∕cB being the physical time delay introduced between the birefringence axes of the medium; c is the speed of light in vacuum, and B is the phase birefringence. It follows from the conventions used in Eqs. (3) and (4) that one gets τz > 0 when y is the fast axis. This implies that the phase birefringence in τz z∕cB has to be defined eff as B neff x − ny . The result, Eq. (4), is valid only under the assumption that the delay τz is a constant for all spectral components of the field. In order to account for the birefringence dispersion we can simply acknowledge the frequency dependence of τ τz in Eq. (2), τz ω z∕cBω, which makes z the constant parameter of the integral: Z Γij z
∞
−∞
W ij ω exp−iωτz ωdω:
(5)
The resulting expression (5) is meant as a replacement for Eq. (2) in Eq. (4) and can raise doubts, as it is clearly no longer a Fourier integral and thus it is not supported by the Wiener– Khintchine theorem. Definition (3) also loses its meaning for
τ τω. There is a physical justification for Eq. (5), however, which we will provide in Section 4. For now it is of interest to notice that Eq. (5) can be reduced back to the Fourier integral form by applying a variable substitution: ω ~ Fω
Bω ω Bω0
(6)
with an arbitrarily chosen effective frequency ω0 . By applying Eq. (6) to Eq. (5) we define a new quantity, which can be called the effective complex mutual coherence function: Γij τz 0
Z
∞
−∞
~ ij ω ~ exp−iωτ ~ z ~ W 0 dω;
−1 dF ω ~ ~ ij ω ~ ~ W F −1 ω; W ~ ij dω
(7a)
(7b)
where τz 0 z∕cBω0 is the effective time delay between the polarization components in the medium at distance z. ~ indicates the inverse of the transformation function deF −1 ω ~ ~ fined by Eq. (6). Both τz 0 and the transformed profiles W ij ω depend on the particular choice of ω0 in a way that makes the overall result independent of ω0 . Transformation formula (7b) is independent of the propagation distance. Representations (7) make sense only if Fω is invertible in the range of optical frequencies where jW ij ωj > 0, i; j x; y. If this condition is met, relation (7a), being a Fourier integral again, carries the same information as Eq. (5). The physical meaning of this correspondence can be expressed as follows: In the case of one-dimensional propagation of a well collimated uniformly polarized field, the influence of the birefringence dispersion on the polarization state of light passing through a linearly birefringent medium can be reproduced by means of a frequency-domain transformation of the spectral polarization (coherence) matrix of the illuminating field. In other words, in place of a medium of length z and birefringence Bω and illuminated by a field described by the matrix Wω, one can take a medium of the same length z with constant birefringence Bω0 and illuminated by light ~ ω, ~ and the state of light at described by a different matrix W the output would be the same. This interpretation, however, has limited use. It follows from Eq. (6) that, when Bω switches its sign within the passband of the field, the domain ~ ij ω ~ will naturally include negof the transformed profiles W ~ ij ω ative frequencies. Therefore, in this case the profiles W ~ can no longer be regarded as standard CSD distributions, which implies that the discussed correspondence is no longer physically realizable. Consequently, the correlation function at the left-hand site of Eq. (7a) is not an analytic signal in general; thus it loses one of the basic properties of the ordinary complex mutual coherence function ([20], part 3.1.3). Nevertheless, it is still a valid expression of correlations in the vectorial field propagating through a medium of the considered type. This is ensured by the equivalence of the integrals of Eqs. (7a) and (5). In Section 5 we illustrate usefulness of the form (7a) in numerical simulations in cases where the negative frequencies occur. The effective complex mutual coherence functions (7a) substituted into Eq. (4) describe propagation in a single birefringent and dispersive element. It follows from the formulas
P. L. Makowski and A. W. Domanski
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
derived in [3–5] that in case of a two-section system, such as the Lyot depolarizer, we need mutual coherence functions of the form Γxy τz2 τz1 . Taking into account the birefringence dispersion, we can express them as Γxy τz2 ω τz1 ω Z ∞ z z W ij ω exp −iω 2 B2 ω 1 B1 ω dω: c c −∞
(8)
It is apparent from Eq. (8) that in order to enforce the linear frequency dependence of the term in square brackets by means of a variable substitution, as it was done for Eq. (5), it is necessary to have B1 ω B2 ω Bω for all ω. In conclusion, in general case of multisection systems the direct usage of the representation in Eqs. (6) and (7) of the effective mutual coherence functions is limited by the two conditions: (a) All the (linearly) birefringent sections have to be described by the same phase birefringence curve Bω and (b) The difference of effective propagation constants Δβω ω∕cBω being proportional to Fω defined by Eq. (6) must not possess local extrema within the bandwidth of the passing light, because Eq. (7b) requires the existence of the inverse of Fω. In cases where condition (b) is not met, an indirect use of representations (7) is still possible. It is done by splitting the initial integral of Eq. (5) into a sum of integrals with neighboring subdomains in which Fω is invertible separately. This sum can be transformed into a sum of integrals of the form ~ intervals (overlapping in general). of Eq. (7a) with finite ω Because the local extremal points of Fω will correspond to singularities in the transformed cross-spectral densities, ~ intervals), one can preEq. (7b) (at the related edges of the ω dict that every local extremum in Fω will cause an oscillatory component in the effective mutual coherence function Γij τz 0 ([20], part 2.4.2). It can be easily verified numerically by calculating the integral of Eq. (5) that a local extremum in Bω results in a correlation function whose absolute value does not approach zero for increasing medium length z.
3. NUMERICAL IMPLEMENTATION As the form of Eq. (7a) represents a Fourier integral again, we can easily approximate it for fast numerical computation. Since the real and imaginary parts of the elements Γij τ essential for determining the state of polarization of the beam often exhibit rapid oscillations within the range of the coherence time of the light wave, the DTFT is a convenient solution: ~ Γij τz 0 ≈ jΔωj
N X
~ ij ω ~ n exp−iω ~ n τz W 0 ;
(9)
n1
where N is the total number of equidistant sampling points ~ n g with the sampling step Δω ~ ω ~N −ω ~ 1 ∕N − 1. The fω ~ n g should be chosen in such a way to make range of fω ~ 1 − Δω∕2; ~ ~ FDω ω ω~ N Δω∕2, where Fω is defined by Eq. (6) and Dω fω∶jW ij ωj > 0g. This ensures full accuracy of expression (9) (exact equality) if the following conditions are strictly fulfilled: (a) we consider τz 0 ∈ −T 0 ∕2; T 0 ∕2, ~ is the period of the transform, and where T 0 2π∕jΔωj (b) the value T 0 is large enough (Δω ~ is set small enough)
1827
~ ij ω ~ n g to eliminate overlapping of the for given elements fW replicas of Γij τz . Because of the continuity in the transform 0 domain there is no periodicity in the signal domain, which implies that the sum of Eq. (9) is actually performed over the ~ but only infinite set of samples covering the entire domain ω, ~ n g contribute to the sum. This fact elimthose within the set fω inates spectral leakage errors. It is formally required that one ~ 0, which restricts the of the samples occupy the point ω choice of ω1 , ωN , and N, but violation of this rule does not affect the shape of the transform within the range of interest: −T 0 ∕2; T 0 ∕2. In order to utilize Eq. (9) in calculation of Jz one has to ~ ij ω ~ n first. obtain the transformed complex distributions W Using the DTFT approximation directly on Eq. (5) in the presence of the birefringence dispersion leads to errors discussed in our previous work [8] for the special case of light linearly polarized at input. These errors originate from the variable distance between the subsequent phases inside the term, exp−iωn τz ωn in Eq. (5), which is inconsistent with the definition of the DTFT. This is why transformation [Eq. (7b)] of the CSD profiles is needed before discretization on a regular grid. ~ ij ω ~ n through direct implementation of Obtaining W Eqs. (6) and (7b) is highly inefficient. However, both evalu~ −1 ω ~ and numerical inversion ation of the derivative d∕dωF of Fω itself can be avoided with little loss of precision if the following transformation procedure is applied. Suppose that Ωω denotes the real or the imaginary part of the original CSD distribution W ij ω. We define a set of samples ωm , where m 1; 2; …; M. If the points fωm g are equidistant, the values fΩωm g can be treated as “frequencies” of a histogram approximating the continuous distribution Ωω, with boundaries of the bins positioned at ωm − Δω∕2 and ωm Δω∕2 for each sample m 1; 2; …; M. Next, we ~ domain and associate center map these boundaries into the ω ~ 0m g, with new bins of a points of the obtained intervals, say fω transformed histogram whose frequencies are of the form ~ 0 ω ~ 0m κm Ωωm g. The coefficients κm Δω∕jFωm fΩ Δω∕2 − Fωm − Δω∕2j ensure an unchanged total area ~ 0 ω ~ 0m g are therefore covered by the bins. The new values fΩ (irregularly distributed) samples of the transformed density ~ ω. ~ They should be resampled on a regular distribution Ω grid with the desired number of points N. The resulting ~ ω ~ n g, n 1; 2; …; N, obtained for the real two data sets fΩ and the imaginary part of W ij ω determine the discrete com~ ij ω ~ n for formula (9). plex distribution W Although the described procedure does not seem to require Fω to be invertible, this restriction still holds. A local ~ ω extremum of Fω implies a singularity in Ω ~ (visible as κm → ∞), which is impossible to represent accurately by a finite number of samples (in the sense of the Nyquist–Shannon sampling theorem). The algorithm would also generate over~ 0 ω ~ 0m g. lapping bins in the new histogram fΩ The procedure is not exact due to the assumption that each sampling point ω ~ 0m corresponds to the central point between the boundaries ωm − Δω∕2 and ωm Δω∕2 mapped to the ~ domain. This is strictly correct only if Fω can be approxiω mated by a linear function within every region covered by every single bin of the histogram fΩωm g. Equation (9) supplemented by the above transformation procedure constitutes the basis of the proposed computational
1828
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
P. L. Makowski and A. W. Domanski
technique. It includes the linear birefringence dispersion dependence into the (effective) mutual coherence functions of the uniformly polarized field passing through a system with one or many linearly birefringent elements meeting the conditions listed in the previous section. These effective mutual coherence functions can be employed directly as replacements for the corresponding terms in some known descriptions based on the complex degrees of coherence [3–5,21]. Those solutions, however, are not in principle dedicated to numerical modeling. It is then of interest to ask how the presented technique relates to the Mueller–Stokes matrix formalism in polychromatic illumination [1,2,5,8], which is equivalent to the spectral coherency matrix formalism [6,7].
4. APPLICATION IN MUELLER–STOKES CALCULUS Every polychromatic field characterized by the frequency profile of the spectral polarization matrix Wω can be equivalently described by the spectral Stokes vector sω [22]. The ordinary (temporal) Stokes vector S of the field is obtained through the integral Z S
∞ −∞
sωdω;
3 3 2 J xx J yy Z ∞ W xx ω W yy ω 6 W xx ω − W yy ω 7 6 J xx − J yy 7 7 7 6 6 4 2ReJ xy 5 −∞ 4 2ReW xy ω 5dω; 2ImJ xy 2ImW xy ω
(10a)
2
(10b)
which reflects the generalized Wiener–Khintchine theorem (5) R evaluated for τ 0: J Wωdω. On the other hand, it follows from representation (1) of the matrix Wω that each spectral Stokes vector sω can be regarded as a second-order statistic of an ensemble of strictly monochromatic waves: fEρ; ω E x ρ; ω; Ey ρ; ωg. It has been confirmed experimentally [23] that the complex amplitudes of these waves are measurable and that they propagate in accordance with Maxwell’s equations. Therefore, it is justified to model propagation of each realization, 2 3 E x ωEx ω E y ωEy ω 6 E x ωEx ω − E y ωEy ω 7 7; sω 6 (11) 4 5 2ReE x ωE y ω 2ImEx ωEy ω such that hsωi sω by use of the Mueller–Stokes matrix formalism: s
out
ω Mωs ω: in
(12)
If we now take the ensemble average of both sides in Eq. (12), we get sout ω
hMωsin ωi
Mωsin ω;
(13)
where Mω is a Mueller matrix of a static medium or an averaged Mueller matrix of a stochastic medium that fluctuates independently from the field. Relation (13) combined with (10a) defines the Mueller–Stokes matrix formalism for any polychromatic illumination uniformly polarized in the space-frequency domain:
Z Sout
∞
−∞
Meff ωsin ωdω;
(14)
where Meff is the averaged effective Mueller matrix of a linear system. Because the resulting vector Sout is meant to describe the whole cross section of the field, we require that it be uniformly polarized in the space–time domain as well. The spectral Stokes vector sin ω represents monochromatic components whose spectral degrees of polarization Pω [9] can take any value between zero and unity. It can be shown that the Mueller–Stokes equations integrated over the optical frequency that we find in works by Billings [1] and Loeber [2] in the context of the Lyot depolarizer are special cases of Eq. (14). This correspondence takes place under the assumption that sin represents the Stokes vector common for all spectral components, which are also spectrally fully polarized and thus expressible by deterministic plane waves. If Meff is additionally derivable from a Jones matrix (a Mueller–Jones matrix [24]), Eq. (14) becomes equivalent to the coherency matrix formalism defined by Barakat [7]. In this technique, the polarization state of the field is represented directly by the coherency matrix Φω, which is essentially the same as Wω defined by the inverse form of Eq. (2), apart from the symmetry of the frequency profiles of the elements; the optical system is represented by the effective Jones matrix, which can express a single train of optical elements or a sum of two matrices corresponding to the two round-trip paths in an interferometer [6]. In order to switch from this technique to the formalism defined by Eq. (14) one has to replace the matrix Win ω by the equivalent spectral Stokes vector sin ω and transform the effective Jones matrix of the system into the corresponding Mueller– Jones matrix [24]. Because of the use of real-valued quantities only and because of the simple algebra of the polarization state transformation, the method of Eq. (14) is well suited for numerical modeling. As for the theoretical analysis, the formalism by Barakat [7] and the Jones formalism are more adequate [25]. Let us now consider propagation through a linearly birefringent element, as discussed in Section 1. We will prove that the output Stokes vector obtained by Eq. (14) is also expressed by elements of polarization matrix (4) for τz τz ω. This will provide physical justification for Eqs. (5)–(7) and allow us to approximate Eq. (14) on the same terms that apply to the DTFT form [Eq. (9)] of the coherence functions Γij τ0z . The Mueller matrix of an uniaxial medium with its optical axis perpendicular to the propagation z axis has the form 2 6 Mδ 6 4
1
3 1
7 7; cosδ sinδ 5 − sinδ cosδ
(15)
where δω Δβωz ωτz ω is the retardation for frequency ω. The placement of the minus sign in Eq. (15) is determined by conventions used in definitions of τz , Γτ, and sω and by the proportionality E ω t ∝ exp−iωt. The magnitude of δω is unlimited, as the spectral degree of polarization [16] of the monochromatic wave sω is unaffected by time shifts introduced between its Cartesian components [9]. We also notice that matrix (15) includes birefringence dispersion to all orders in a natural way, but
P. L. Makowski and A. W. Domanski
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
for the moment let us assume that the birefringence in the medium is nondispersive, meaning τz constant for all ω. If we substitute Meff ω Mδω into Eq. (14), then according to definitions (10) the two nontrivial output Stokes parameters take the forms Z S out 2 2
∞
−∞ ∞
Z
−2
−∞
Z S out 3 2
∞
ImW xy ω sin−ωτz dω;
(16a)
ReW xy ω sin−ωτz dω
−∞ ∞
Z
2
ReW xy ω cos−ωτz dω
−∞
ImW xy ω cos−ωτz dω:
(16b)
On the other hand, it follows from Eq. (2) that Γxy τz
Z
∞
−∞
ReW xy ω iImW xy ω
× cos−ωτz i sin−ωτz dω:
(17)
Comparing Eqs. (16) with Eq. (17) leads to the immediate result 3 2 Γxx 0 Γyy 0 6 Γxx 0 − Γyy 0 7 7 (18) Sout z 6 4 2ReΓxy τz 5; z 2ImΓxy τ which is equivalent to Eq. (4), as expected. In the case of a out multisection system we would obtain a vector P S composed of linear combinations of functions Γxy k ak τk with ak 0; 1 and τk τzk , where k numbers the birefringent sections. It is important to note that the established correspondence is satisfied both for the nondispersive τz constant and for the dispersive case, τz τz ω, where the latter case entails replacement of Γxy τz by Γxy z defined through Eq. (5). This is because for τz ω z∕cBω the integrals in Eqs. (16) can be collected to the form of the real and imaginary parts of the right-hand site of Eq. (5). Two essential conclusions follow from this fact. First, definition (5) of the modified coherence function Γxy z is given clear physical sense, as it leads to the same results as the physically well justified Mueller–Stokes equation (14). Second, because the spectral vector sω is composed of the functions W ij ω, the numerical approach for propagation of polychromatic light in an uniaxial medium expressed by Eq. (9) translates into ~ Sout ≈ jΔωj
N X
~ n τz ~ n ; Mω sin ω 0 ~
(19)
n1
which, unlike Eq. (9), can be easily extended for the case of multisection systems composed of linearly birefringent parts and possibly other nonbirefringent elements: ~ Sout ≈ jΔωj
N X n1
z2 zK in ~ n ; τz1 Meff ω ~s ω~ n ; 0 ; τ0 ; …; τ 0
(20)
1829
where k 1; …; K numbers the linearly birefringent elements in the optical system. Despite the fact that Eq. (20) does not resemble any kind of discrete Fourier transformation, it is subject to the same conditions as Eqs. (9) and (7). We keep ~ defined through Eqs. (7), will not always in mind that s~ in ω, represent the standard spectral Stokes vector, as formula (7b) can introduce negative frequencies to the transformed profiles. In the following section we examine the accuracy of expression (20) supplemented by the transformation algorithm described in Section 3. The simulation errors will be evaluated by comparing the output (temporal) degree of polarization predicted by Eq. (20) with the one obtained through direct numerical integration using Eq. (14).
5. NUMERICAL EXAMPLE For the sake of clarity we choose a simple two-section system (the fiber-optic Lyot depolarizer [26]) and uniformly coherent light [9] with unit spectral degree of polarization Pω for all ω, in particular a linearly polarized light wave described completely by its power spectral density distribution Ψω W xx ω W yy ω and the initial azimuth angle α. In this case Eq. (20) takes the form Sout τ
0
~ ≈ jΔωj
X N
~ ~ n ; τ0 ; 2τ0 Ψω ~ n sˆ in ; Meff ω
(21)
n1
~ ω where Ψ ~ is the original PSD Ψω after the domain transformation defined by Eqs. (6) and (7), sˆ in sin α∕sin 0 α is the normalized spectral Stokes vector of the input field linearly ~ n ; τ0 ; 2τ0 polarized with the azimuth angle α, and Meff ω ~ n τ0 is the effective Mueller matrix of M2ω~ n τ0 R45°Mω the depolarizer. In order to obtain the result in the coordinate system associated with the first section, the second matrix has to be surrounded by two rotation matrices: ~ n τ0 R45°. The set of sampling points fω~ n g R−45°M2ω meets analogical conditions as discussed for Eq. (9). The simulations were performed using the PSD of the superluminescent Broadlighter diode Q1350-HP from Superlum, Ltd. with spectral width of about 200 nm. The two sections of the depolarizer are made of a highly birefringent holey photonic crystal fiber [27], single-mode within the spectrum width of the Q1350-HP diode. Its phase birefringence, depicted in Fig. 1, switches sign within the spectrum range of the diode. Figure 2 depicts the result of the transformation defined by Eq. (7b) and the procedure described in Section 3. The original power spectrum samples fΨωm g, m 1; …; M (black points), are replaced by samples of the transformed profile ~ ω ~ n g, n 1; …; N (gray points), which encodes both the fΨ diode spectrum and the phase birefringence curve of the fiber. At the same time, the original curve Bω (black line) is replaced by a constant effective birefringence Bω0 (gray line). The effective frequency has been chosen arbitrarily at the right border of the Q1350-HP diode spectrum, ω0 ωM ; hence Bω0 takes the highest value of Bω. The number of the transformed samples N has been chosen in such a way to make Γ11 τz 0 fade almost completely for τ 0 < T 0 ∕2, where T 0 is the period of DTFT transforms (10). Because the phase birefringence Bω switches sign within the passband of Ψω, ~ ω the transformed samples fΨ ~ n g occupy both positive and
1830
J. Opt. Soc. Am. A / Vol. 30, No. 9 / September 2013
Fig. 1. Input data: phase and group birefringence of the highly birefringent holey photonic crystal fiber plotted in the region covered by the Q1350-HP diode spectrum. The phase birefringence Bλ is represented by a fifth-degree polynomial fitting to the samples of the original phase birefringence curve [27].
Fig. 2. Input data: PSD samples of the superluminescent diode source along with associated phase birefringence curve of the highly birefringent photonic crystal fiber (black) and the corresponding functions after domain transformation (7b) (gray).
negative frequencies. This is an expected property of Eqs. (6) ~ simply switches the phase and (7b); the negative frequency ω sign inside the exponential term in (7a), just like the negative sign of the phase birefringence Bω switches the phase sign inside integral (5). The two integrals have to be equivalent, hence the necessity for the negative frequencies in the domain ~ ω ~ in this example. of Ψ Figure 3(a) presents the residual degree of polarization P 1∕S 0 S 21 S 22 S 23 1∕2 [obtained from Eq. (21)] of light emerging from the fiber-optic depolarizer of total length L L1 L2 , where L2 2L1 . The horizontal axis α represents the azimuth angle of the input linear polarization. The white arrows indicate the total length L L1 L2 such that τL2 T 0 ∕2; therefore they determines the end of the 0
P. L. Makowski and A. W. Domanski
expected range of accurate results (the analogue of the Nyquist frequency). In the upper region of the diagram in Fig. 3(a) we observe a nonphysical compensation of depolarization, which corresponds to τL2 T 0 ; it is then an expected 0 manifestation of the periodicity of DTFT (9), which is implicitly present in expression (21), according to the correspondence established in Section 4. Figures 3(b) and 3(c) express the inaccuracy of formula (21) for M 30 and M 50, while the number of transformed samples N 30 remains unchanged (along with the value of the Nyquist delay). The reference degree of polarization values were obtained by direct numerical evaluation of integral (14) corresponding to the sum of Eq. (21) by use of the quad_qags procedure in Maxima CAS (v5.24.0, GCL 2.6.7). ~ ω ~ n g for The discrete spectra fΨωm g transformed to fΨ proposed solution (21) and the corresponding curves Ψω for reference method (14) were prepared in such way to allow evaluation of errors originating from the approximate character of the domain transformation procedure described in Section 3. First, M 50 regularly distributed samples were chosen from the original PSD curve taken from the documentation of the Q1350-HP diode. Then the cubic-spline interpolation of these points was calculated to serve as a simulation of a real-world continuous spectrum passed to reference formula (14) [and used to generate Fig. 3(c)]. Next, a set of M 30 equidistant samples was picked from the interpolation curve, and the second cubic-spline interpolation was calculated from it for Eq. (14) [Fig. 3(b)]. This means that the two diagrams in Figs. 3(b) and 3(c) actually correspond to two separate light sources with slightly different PSDs Ψω, but in return we can be certain that in both cases (M 30 and M 50) the discrete functions Ψωm represent their continuous counterparts Ψω with comparable accuracy. Therefore, if proposed technique (21) turns out to be less accurate for M 30, it is not due to incomplete representation of the details in the PSD profile, but rather due to violated assumptions of the domain transformation procedure. As expected, we observe that the accuracy of the proposed simulation technique remains high as long as τ0L2 < T 0 ∕2. This confirms implicit presence of DTFT in formula (21). On the other hand, the slight discrepancies visible below the T 0 ∕2 line strongly depend on the initial number of the PSD samples: the change from M 30 to M 50 reduces
Fig. 3. Results: (a) degree of polarization of superluminescent diode light after passing a dispersive Lyot depolarizer calculated by the proposed technique, Eq. (21), and its discrepancy from the result of exact numerical integration, Eq. (14), for initial number of the spectrum samples, (b) M 30, and (c) M 50. The white arrows indicate the total length L1 L2 of the depolarizer that corresponds to the half-period T 0 ∕2 of the DTFT approximation of the autocorrelation Γ11 τL2 0 .
P. L. Makowski and A. W. Domanski
the maximal level of jΔDoPj (below the arrow), where DoP is degree of polarization, by a factor of 4, from less than 2 × 10−3 to less than 5 × 10−4 . These values reflect expected limitations of the approximate domain transformation procedure defined in Section 3 for the discretized CSD functions. The estimated absolute integration errors reported by the quad_qags algorithm did not exceed the level of 1 × 10−8 ; therefore the inaccuracies of the reference method can be ignored. In the more general case of the field described by the matrix Wω the domain transformation has to be performed on all four elements of the matrix. In the case of the complex-valued element W 12 ω W 21 ω it is performed for the real and the imaginary part separately.
6. CONCLUSION We have discussed conditions under which an uniformly polarized polychromatic field can be numerically propagated in dispersive birefringent systems by the use of the DTFT of the discretized frequency profile of the CSD matrix Wρ; ρ; ωn . The proposed numerical approach utilizes a form of the frequency-integrated Mueller–Stokes equation that involves spectral Stokes vectors. The dispersion dependence of birefringence of the (homogeneous) linearly birefringent elements in the system is included to all orders by means of a domain transformation of the frequency profiles of these vectors at the input plane. The sources of errors in the proposed approximation have been identified, and the expected accuracy range has been numerically confirmed by using parameters of a photonic crystal fiber Lyot depolarizer and wideband superluminescent diode source. The technique can be useful in modeling and design of some types of bulk optics devices as well as polarimetric and interferometric fiber-optic sensors.
REFERENCES B. H. Billings, “A monochromatic depolarizer,” J. Opt. Soc. Am. 41, 966–968 (1951). 2. A. Loeber, “Depolarization of white light by a birefringent crystal. II. The Lyot depolarizer,” J. Opt. Soc. Am. 72, 650–656 (1982). 3. W. K. Burns, “Degree of polarization in the Lyot depolarizer,” J. Lightwave Technol. 1, 475–479 (1983). 4. K. Mochizuki, “Degree of polarization in jointed fibers: the Lyot depolarizer,” Appl. Opt. 23, 3284–3288 (1984). 5. P. L. Makowski, M. Z. Szymanski, and A. W. Domanski, “Lyot depolarizer in terms of the theory of coherence-description for light of any spectrum,” Appl. Opt. 51, 626–634 (2012). 6. E. I. Alekseev and E. N. Bazarov, “Theoretical basis of the method for reducing drift of the zero level of the output signal of a fiber-optic gyroscope with the aid of a Lyot depolarizer,” Sov. J. Quantum Electron. 22, 834–839 (1992). 7. R. Barakat, “Theory of the coherency matrix for light of arbitrary spectral bandwidth,” J. Opt. Soc. Am. 53, 317–322 (1963).
1.
Vol. 30, No. 9 / September 2013 / J. Opt. Soc. Am. A
1831
8. P. L. Makowski and A. W. Domanski, “Approach for modeling influence of birefringence dispersion on polarization properties of multi-section systems,” Proc. SPIE 8697, 869703 (2012). 9. P. Réfrégier, T. Setälä, and A. T. Friberg, “Temporal and spectral degrees of polarization of light,” Proc. SPIE 8171, 817102 (2011). 10. T. R. Woliński, “Polarimetric optical fibers and sensors,” in Progress in Optics, E. Wolf, ed. (Elsevier, 2000), Vol. 40, pp. 1–75. 11. T. Nasilowski, F. Berghmans, T. Geernaert, K. Chah, J. Van Erps, G. Statkiewicz, M. Szpulak, J. Olszewski, G. Golojuch, T. Martynkien, W. Urbanczyk, P. Mergo, M. Makara, J. Wojcik, C. Chojetzki, and H. Thienpont, “Sensing with photonic crystal fibres,” in IEEE International Symposium on Intelligent Signal Processing, 2007. WISP 2007 (2007), pp. 1–6. 12. D. J. Liu and Z. X. Zhou, “Propagation of partially polarized, partially coherent beams in uniaxial crystals orthogonal to the optical axis,” Eur. Phys. J. D 54, 95–101 (2009). 13. R. C. Jones, “A new calculus for the treatment of optical systems. VII. Properties of the N-matrices,” J. Opt. Soc. Am. 38, 671–683 (1948). 14. N. Ortega-Quijano and J. L. Arce-Diego, “Generalized Jones matrices for anisotropic media,” Opt. Express 21, 6895–6900 (2013). 15. R. M. A. Azzam, “Propagation of partially polarized light through anisotropic media with or without depolarization: a differential 4×4 matrix calculus,” J. Opt. Soc. Am. 68, 1756–1767 (1978). 16. E. Wolf, “Unified theory of coherence and polarization of random electromagnetic beams,” Phys. Lett. A 312, 263–267 (2003). 17. F. Gori, M. Santarsiero, R. Simon, G. Piquero, R. Borghi, and G. Guattari, “Coherent-mode decomposition of partially polarized, partially coherent sources,” J. Opt. Soc. Am. A 20, 78–84 (2003). 18. J. Tervo, T. Setälä, and A. T. Friberg, “Theory of partially coherent electromagnetic fields in the space-frequency domain,” J. Opt. Soc. Am. A 21, 2205–2215 (2004). 19. M. A. Alonso and E. Wolf, “The cross-spectral density matrix of a planar, electromagnetic stochastic source as a correlation matrix,” Opt. Commun. 281, 2393–2396 (2008). 20. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995). 21. P. L. Makowski and A. W. Domanski, “Extended Mueller–Stokes description of polarization-mode transformation in linearly birefringent single mode optical fibers,” Opt. Lett. 38, 1107–1109 (2013). 22. O. Korotkova and E. Wolf, “Generalized Stokes parameters of random electromagnetic beams,” Opt. Lett. 30, 198–200 (2005). 23. A. Dogariu and G. Popescu, “Measuring the phase of spatially coherent polychromatic fields,” Phys. Rev. Lett. 89, 243902 (2002). 24. C. Brosseau, “Mueller matrix analysis of light depolarization by a linear optical medium,” Opt. Commun. 131, 229–235 (1996). 25. C. Tsao, Optical Fibre Waveguide Analysis (Oxford University, 1992). 26. K. Böhm, K. Petermann, and E. Weidel, “Performance of Lyot depolarizers with birefringent single-mode fibers,” J. Lightwave Technol. 1, 71–74 (1983). 27. L. Zhang, T. Luo, Y. Yue, and A. E. Willner, “High group birefringence in photonic crystal fibers with both positive and negative phase birefringences,” Pure Appl. Opt. 10, 035004 (2008).
