Statistical studies on optical vortices in dynamic speckle fields

New parameters to statistically describe and differentiate between different decorrelation behaviors in dynamic speckle fields are described. These decorrelation behaviors are surrogate descriptors of the dynamics of the underlying processes in object space being observed. The statistical parameters are based on the temporal variations in the location of optical vortices in the speckle fields. The length and number of optical vortex trails, motion of the vortices in the plane of observation and the distance between the mean locations of the positive and negative vortices are investigated. The implementation of the statistical analysis presents new methods to quantify and describe biophysical dynamics. © 2018 Journal of Biomedical Photonics & Engineering.


Introduction
Random optical interference gives rise to granular patterns with bright and dark regions called speckle.Light scattering techniques, including those involving laser speckle analysis have been used to study dynamic biological systems over the last several decades [1][2][3][4].Speckle based imaging techniques have seen tremendous advancements since the 1990s and are now well-established in the field of biomedical optics [5][6][7].Within speckle fields are locations of phase discontinuity where the intensity of the field is zero and the phase is undefined.Such phase singularities were first described by Nye and Berry [8].A few years later they reported an analysis of phase singularity lines within fields of multiple beam interference [9].Detailed studies regarding the statistics of the amplitude distribution in such interference fields were conducted by Baranova et al [10].Although the existence of these phase singularities have been known since the 1970s, their in-depth study in speckle fields [11] has garnered interest much more recently.Speckle fields, ( , ) t r , like other electromagnetic fields, can be written as a general solution of the Helmholtz equation as , , Α ,t , i t t e r r r (1) where Α(, t) is the amplitude and (, ) is the phase of the field, both at location r and time t.The real and 27 May 2018 © J-BPE 020301-3 imaginary parts of this speckle field can be separated, with phase singularities occurring in locations where both the real and imaginary parts are zero (Fig. 1).Fig. 1 Zero-contours of real and imaginary field components.The red (solid) lines represent where Re ψ ( ) = 0 and the black (broken) lines represent where Im ψ ( ) = 0.The points of their intersection are the locations of phase singularities.
In the immediate vicinity of and surrounding these singular points, the optical phase rotates in either a clockwise or counter-clockwise direction over at least a full 2π radians and the structures are subsequently referred to as either negatively charged or positively charged optical vortices, respectively [12].The number of 2π rotations of the phase is known as the topological charge of the vortex, expressed mathematically as where the line integral is along a closed loop, l, around the vortex.Optical vortices always occur in pairs of opposite topological charge, and as a consequence of conservation of charge, there must always be an equal number of positively and negatively charged vortices in the scattered field It should be noted here that the generation of optical vortices does not require special circumstances.It is a phenomenon that accompanies any interference of randomly scattered plane waves [13].
However, experimentally, it has been observed that helically phased beams play a vital role in the formation of phase singularities [14].These beams, particularly the Laguerre-Gaussian mode (LG 01 ) have been shown to be produced from a combination of Hermite-Gaussian modes HG 01 and HG 10 [15].An alternative method of generation was demonstrated in 1990 where a diffractive element such as a grating with a dislocation was used to obtain a helically phased mode from the Gaussian beam of a laser [16].Subsequently, the generation of optical vortex beams has also been achieved using spiral phase plates [17] and spatial light modulators [18].
While the angular divergence of the wave-front was already described by Baranova et al [10], it was noted by Allen et al. [19] that helically phased beams carried orbital angular momentum with them.This was in addition to the already commonly known momenta associated with the polarization state and photon spin.A helical mode whose azimuthal component has the profile e iφ with φ being the phase, has been shown to have an angular momentum of lℏ per photon.This fact that the helical beams contained angular momentum, has been the basis of multiple applications of optical vortices involving trapping and rotating micro-particles [20][21][22].Further advancements were done when particles were made to revolve in an annular ring by using a circularly polarized LG 01 beam with a large azimuthal component [23].Developments in remote sensing and wireless communication have also focused on these beams containing orbital angular momentum [24].Additionally, analyzing the spatial coherence in a one-dimensional projection of the optical vortex field has been used as a method to delve into the information content of an optical vortex [25][26].
Thus, we note that the field of singular optics is a rapidly expanding one, both as a science for understanding the properties of beams creating these fields, as well as in engineering to apply this understanding.It has been previously shown [27] that the optical vortices in an interference field translate in both time and space due to an interplay of three factors, the amplitude gradient, phase gradient and intensity gradient.However, these motions have not yet been correlated to the dynamics of the physical systems which caused the interference.The work presented in this paper is targeted towards gaining a better understanding of the behavior of these optical vortices and their potential use in understanding dynamical biological systems.

Methods
Sequences of 50 speckle patterns with user-defined correlation behaviors were simulated using the concept of a copula [28] in the MATLAB environment.Decorrelation rates are defined in terms of the number of frames it takes for the overall speckle pattern autocorrelation coefficient, Γ I , to reduce to a value of 1/e.Three different rates of decorrelation behavior between the frames were studied, with their autocorrelation coefficients falling to 1/ e in 11, 6 and 3 frames, during the 50 frame simulation.Herein, they shall be referred to as the slow, medium and fast rates of decorrelation, respectively.For each of the rates, three different decorrelation behaviors (or modes) were studied.One mode was of constant sequential autocorrelation coefficient ( Γ I = 0.9984 for slow  decorrelation rate; Γ I = 0.99 for medium decorrelation rate and Γ I = 0.96 for fast decorrelation rate) between frames.The second behavior was one that displayed a Gaussian decorrelation line shape.This decorrelation behavior has generally been understood to be associated with ordered dynamics [29].The third behavior examined, Lorentzian decorrelation, is linked to Brownian dynamics.In general, biological systems likely exhibit a combination of ordered and Brownian (or disordered) dynamics.
Once the 50 speckle patterns were generated for each decorrelation behavior studied, we created a pseudo-phase representation for each of the patterns using a two-dimensional Hilbert filter (Fig. 2(b)) [11].The location of optical vortices in the speckle were identified and tracked on a frame-by-frame basis through these sequences (Fig. 2(c)).
This was done using a series of convolution operations on the phase distribution over the speckle fields as shown in Eq. (3) and Eq. ( 4) [11].This series of convolution operations gives non-zero values only at the locations of phase singularities.These values are the topological charge n t at those locations.
where φ(x, y) represents the pseudo-phase of the speckle field, ⊗ is the convolution operator and It is to be noted that Eq. ( 3) represents a discrete evaluation of the topological charge that is expressed in Eq. ( 2).The location of a positive (negative) topological charge is called a positive (negative) vortex.Physically, the two different charges represent the two opposing senses of rotation, clockwise and anti-clockwise, of the phase of the speckle field in the complex plane.Tracking any particular vortex, as it traverses through the frames, results in what has been termed a vortex trail [11].
Physically, coherent light scattering from dynamic biological systems over a period of time leads to temporal decorrelation of the speckle pattern in the observation plane.In our simulations, the individual frames can be considered as snapshots of the speckle field as it decorrelates over time.Thus, traversing through frames is equivalent to the temporal movement of a vortex.In this context, frame numbers can be considered as discrete points in time.
Identifying the x and y coordinates of each optical vortex as it passes through the frames, gives rise to a trail from its point of generation to its point of annihilation.Occasionally, a trail also ends due the vortex leaving the field of observation.A trail is called either a positive or a negative trail, according to the sign of the topological charge contained in the vortex that created it.When a negative and a positive vortex intersect at a single spatial and temporal location, they annihilate each other and both trails that they were forming end [11].
A challenge in the simulations discussed herein was to draw the boundary between a continuous trail and a new trail starting close to where a trail ended in the previous frame.As the locations of the vortices were defined in the terms of discrete pixels, it was essential to decide how close two vortices in adjacent frames must be to be considered the part of the same trail.Based on empirical observations of multiple simulation results, it was determined that a change of four or more pixels in any direction would be considered to be the beginning of a new trail.
To create the vortex trail displays like the ones shown in Fig. 3, the vortices in each frame were located as described above, using the series of convolution operators.The location of each vortex, found in this fashion, was then stored on a frame-by-frame basis.
27 May 2018 © J-BPE 020301-5 While this enabled us to obtain the visual representation of the vortex trails, the following algorithm was followed to investigate the statistics of the trails.The coordinates of the vortices found in each frame were stored in an individual list.This was followed by locating the first vortex in the first frame (or equivalently, the first list) and searching for a vortex in the second frame such that it was within four pixels in any direction of the coordinates of the previous vortex found.If such a vortex did not exist, then the trail ended and a new trail began with the next vortex in the first frame.If a vortex was found within the pre-defined four pixel radius, we moved to the third frame and searched for a vortex within the same radius with the coordinates of the second frame as the center location.In this manner, each trail continued until no vortex was found within the set radius of 4 pixels, or all the frames were exhausted.After all the trails starting from the first frame were tracked, the algorithm searched for trails starting from second frame using vortices that were not already considered to constitute any of the previously formed trails.This process continued through the final frame.A pictorial description of this algorithm is shown as a flowchart in Fig. 4. Once the individual trails were tracked, various parameters of the motion of the vortices were investigated.One parameter investigated was the average length of the positive and negative trails, as well as the overall average length of all trails regardless of charge.This was defined in terms of the number of frames the average trail survived.Another parameter studied was the average displacement of the vortices, as they traversed through the frames (equivalently, through time).For this, the positive and negative centers of the frames were determined.A positive (negative) center ( X c ,Y c ) of a frame is defined as the mean location of all the positive (negative) vortices in the frame: where (x i , y i ) is the coordinates of the i th vortex in the frame and n is the total number of positive (negative) vortices in the frame.The overall center of the frame ( X co ,Y co ) is the weighed mean location of the positive and negative centers:  and n n represent the number of positive and negative vortices in the frame, respectively.We defined the mobility, M , for each type of vortex as the mean distance moved by their center, per frame (Eq.( 7)): where M cq represents the mobility of the q th type of vortex and i X cq , i Y cq ( ) is the location of the center of q th type of vortex in the i th frame.As indicated, q can represent positive (p), negative (n) or overall (o).N is the total number of frames in the simulation.Additionally, the distance between the positive and negative centers of each the frame was studied.We have termed this distance as the charge separation, G(i).
Poincaré plots [30] were used to analyze the charge separation.Poincaré plots are used to display consecutive measurements against each other, i.e., the i th element of a series X, X(i) is plotted against the (i-1) th element, X(i-1).This is continued for all set of consecutive data points on X.
Once all the data points are plotted as above, standardized metrics SD1 and SD2 for the data set are measured [31] with the following definitions where Var(X) stands for the variance of the data series X.
An ellipse fitting process is then employed over the Poincaré plots.This ellipse has semi-major axis length of SD2 and semi-minor axis length of SD1 with the major axis inclined at 45 o to the horizontal axis of the Poincaré plot (Fig. 7).The metric SD1 represents the short-term variations of the quantity measured (such as standard deviation in successive differences), while SD2 is what remains when contribution of SD1 is removed from the variance (Eq.7) and thus indicates the longterm variations.The shape of the obtained ellipse indicates the comparative effects of each type of variation.The more dominant long-term variations are, the higher is the magnitude of the semi-major axis (given by SD2), resulting in a more elongated ellipse.

Results and Discussion
The average trail lengths and number of trails for each type of vortex is noted in Table 1.
The effects of decorrelation line shapes and rates were investigated (see above).By weighting the trail length of each type (positive and negative) by the corresponding number of trails, the overall average trail length of each type and rate of decorrelation was obtained.These results are displayed in Fig. 5.  difference of means was used for all statistical significance testing in this study.Additionally, it was also observed that as the rate of decorrelation increases, the average trail length of each mode became shorter.
For each type of decorrelation, the number of vortices in each frame roughly remained independent of all other simulation conditions.The number of vortices identified was on the order of half of the number of coherence areas in each speckle field [32].Thus, longer trails directly resulted in fewer number of trails.In addition to elongated trail lengths, slower decorrelation was also found to be associated with lesser number of individual trails.The third mode, constant rate of decorrelation ( Γ I = 0.9984 for slow, Γ I = 0.99 for medium and Γ I = 0.96 for fast decorrelation, respectively), formed the shortest, and by extension, the most number of trails at the medium and fast decorrelation modes.At the slowest decorrelation rate, the constant decorrelation mode formed trails longer than the Gaussian mode, but shorter than the Lorentzian mode.
Fig. 6 Average vortex displacement (or mobility) for each mode, under three different rates of decorrelation (same as Fig. 5).The statistics are from 20 measurements, each measurement comprising of all 50 frames initiated from a randomly generated speckle pattern; Error bars represent one standard deviation The next parameter that was investigated was the mobility, M , of the vortices.Mobility was defined in terms of the average displacement incurred by the center of each type of vortex field, while moving between consecutive frames.From Table (2), we can observe that the vortices experienced higher mobility under Gaussian decorrelation, as compared to Lorentzian decorrelation (p < 0.01).Additionally, increasing the rate of decorrelation resulted in higher motion of the vortices (p < 0.01).Combining the results from both the tables above, it can be seen that the longer trails were associated smaller displacements.
To study the variations in charge separation, Poincaré plots [7] were employed.As noted above, these plots are used to investigate the long-term and short-term variations in a variable.In this case, the variable in question was the charge separation over the 50 frames.A typical Poincaré plot from this set of measurements is shown in Fig. 7.It was observed that slow decorrelation resulted in elongated ellipses on the Poincaré plot, which is a result of higher SD2/SD1 ratio.As mentioned earlier, this is the ratio of the long-term variability to the short-term variability.Thus, we notice that as the decorrelation rate increases, the short-term variations tend to become more prominent as compared to the long-term variations.This ratio tends towards a value of 1.0 at higher decorrelation rates.The SD2/SD1 ratios can be visualized by placing the obtained ellipses in overlap with each other, as shown in Fig. 8.
It can be seen that the use of Poincaré plots of charge separation gives a stronger indication of the rate of decorrelation, rather than the lineshape of the decorrelation.Thus, this can be used as an indicator of the rate of dynamics in a system as is discussed below

Conclusion
In this paper, a new approach to studying dynamic systems by tracking phase singularities (or optical vortices) has been presented.Speckle patterns for biological systems inherently contain phase singularities, or optical vortices, points around which the phase of the field rotates.In this work we have related the observed dynamics of the optical vortices in a simulated speckle field, to the type and rate of decorrelation in the field.This decorrelation is a representative of the dynamic behavior of scatterers that, under coherent illumination, produced the speckle field under observation [29].
Analyzing the location of vortices at discrete time points in a speckle field simulated using a pre-defined rate of decorrelation, we note the expected result that as the rate increased, the lengths of the vortex trails decreased.In this study, the rate of speckle decorrelation was assumed to correspond to the rate of activity in the scattering object.Thus, in a scattering medium with fast motion among particles, the individual phase singularities tend to survive for shorter periods of time.Faster activity was also shown to result in a higher degree of vortex displacement per frame in the location of the singularities.Additionally, we compared different modes of activity; a Lorentzian decorrelation relationship between frames, corresponding to Brownian or unordered form of activity among the scatterers and a Gaussian decorrelation relationship, corresponding to an ordered mode of activity among the scatterers.The ordered motion resulted in more mobile phase singularities as they traversed through time.This also corresponds to shorter durations for which the trails existed.We remind the readers that the mobility, M, is defined as the average displacement per frame of the mean location (center) of the vortices.A possible reason for low vortex mobility for Brownian motion could be the homogeneous line profile of the scatterer distribution in a Lorentzian flow model [33].This homogeneity results in no preferential direction for the vortices to move.This results in the average summed motion to be less than the Gaussian flow model, which has an inhomogeneous line profile of the scatterers [33].We also examined a technique of analyzing sequential vortex data using Poincaré plots, by separating the long and short-term variations and calculating their relative prominence with regards to each other.When applied to the average charge separation per frame, this technique did not indicate any significant difference between the different modes of decorrelation.However, reducing the rate of decorrelation did result in the Poincaré ellipses elongating, indicating a higher prominence of the longterm variations.A possible reason could be that faster decorrelation did not allow sufficient time for the longterm trends in the data to dominate over the short-term trends.We note here that at no point at all do the shortterm variations SD1 dominate over the long-term variations SD2 .In the case in which SD1 > SD2 the ellipse would be more "broad" than "long', i.e., the major and minor axes would have reversed.Thus, the use of Poincaré plots was more useful in comparing the different rates of activity, rather than the different modes of activity.
Note that actual phase information is not available in the intensity speckle patterns, and thus we took the approach of calculating pseudo-phase estimate maps 27 May 2018 © J-BPE 020301-9 using convolution operators [11].While this approach will not necessarily re-create the true phase of an imaged speckle pattern, the behavior of the vortices is not affected.Also, the simulation only takes into account a limited spatial window of the speckle field.Thus, there are instances of minor topological charge imbalance (about 1 extra vortex in 4000 total vortices) in some of the simulations.This usually occurred in situations where a vortex was close to the edge of the observation window, and the opposite topological charge corresponding to this lone vortex was just outside the window.This oppositely charged vortex typically appeared within the window in the very next frame, thus balancing out the charges.This can be thought of as a case where a pair production happens right at the edge of the window, but only half the pair is within the window at that instant.It is also to be noted we are currently not in a position to speculate on the physical significance of some of the parameters noted here (such as the charge separation, G(i)).As such, they are introduced only as a measureable feature with some sensitivity to the simulation dynamics (in the case of G(i), sensitivity to the rate of activity).
We conclude by noting that the analysis in this paper can potentially be used to add to our current understanding of light scattering and propagation, the physics of coherent wave fields, and the field of singular optics.This can, in turn be used to improve existing microscopy and imaging techniques, particularly those aimed at quantifying biophysical dynamics.Simulating additional biological features such as tumors and other complex systems containing dynamics scattering elements, and identifying them using the dynamics of optical vortices as done in this paper, is among the possible directions in which this research can progress.

Disclosures
The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.

Fig. 2
Fig. 2 Locating optical vortices in a speckle field.(a) A speckle image.(b) Its pseudo-phase representation.(c) Location of positive (blue circles) and negative (red points) vortices.The colorbar in (a) indicates relative intensity while those in (b) and (c) indicate phase.The axes represent spatial coordinates.

Fig. 3
Fig. 3 Typical displays of vortex trails, as the vortices (red points: negative, blue circles: positive) are tracked through the 50 frames (a).Slow decorrelation with autocorrelation coefficient falling to 1/ e in 11 frames (b).A faster decorrelation with autocorrelation coefficient falling to 1/ e in 6 frames.

Fig. 4
Fig. 4 Flow chart representing the algorithm followed to track individual optical vortex trails.

Fig. 5
Fig. 5 Average trail length for each decorrelation type, under three different rates of decorrelation.(Slow: (1/e) = 11 frames; Medium: (1/e) = 6 frames; Fast: (1/e) =3 frames; Total frames in simulation = 50.The statistics are from 20 measurements, each measurement comprising of all 50 frames initiated from a randomly generated speckle pattern; Error bars represent 1 standard deviationThe Lorentzian decorrelation behavior typically resulted in longer and fewer trails compared to Gaussian decorrelation, given the same rate of decorrelation (p < 0.01) for both the slow and fast rates of decorrelation shown in Fig.5.The Student's t-test for

Fig. 7 A
Fig. 7 A typical Poincaré plot obtained for charge separation G i ( ) .The magnitude of descriptor SD1 indicates short-term variability in the value of the separation, while SD2 indicates the long-term variability.In this plot, SD1 = 1.01 and SD2 = 1.85.The lengths are in terms of coordinate points.This simulation was done for a 64 × 64 pixel speckle field.

Fig. 8
Fig. 8 Ellipses from the Poincaré plots of charge separation G(i) for (a) Gaussian (b) Constant and (c) Lorentzian decorrelation behavior.(O : Slow decorrelation; X : Medium decorrelation; -:Fast decorrelation).It can be noted that the ellipses are prominently elongated for the slowest rate of decorrelation, indicating a high SD2/SD1 ratio.This elongation of the ellipse indicates a greater effect of long-term variations on the conditions under which the measurement has been made.

Table 1
Average trail lengths and number of trails of each type (C: Constant; G: Gaussian; L: Lorentzian).The trail lengths are given as the number of frames the trail lasted.Total number of frames in the simulation = 50.

Table 2
The average displacement of each type of vortex center as they travel through the frames (C: Constant; G: Gaussian; L: Lorentzian).All distances are given in terms of coordinate points.Each frame was a square of 64 × 64 coordinate points.