1. Introduction

The spatial coherence properties of light are important—not least in several industrial applications, such as microlithography, laser ablation and annealing—where non-coherent high power lasers are used. Conventionally, the spatial coherence properties of an optical field are determined by carefully conducted interference measurements [1–4]. Unfortunately the classic interference measurement methods are often both cumbersome and time consuming. Therefore, increasing attention is paid to the problem of how to retrieve spatial coherence properties solely from information about the intensity distribution of the optical field. This information should be readily available, e.g., by measuring the intensity distribution in a finite number of planes transverse to the direction of beam propagation [5–9]. The algorithm introduced in this work operates much in the same manner as phase retrieval algorithms, such as the Gerchberg-Saxton algorithm [10] or the Yang-Gu algorithm [11], for the recovery of the phase distribution from the intensity information in a coherent beam. It should also be noted that the term intensity will be used in the common and convenient, though somewhat inadequate, sense. To be sure, for a partially coherent field the instantaneous intensity is very rapidly fluctuating. In this work, the term intensity refers to the time average of the instantaneous intensity. This is also the entity that is measured experimentally, since for most practical cases the integration time of the measurement equipment is far larger than the timescale of the fluctuations of the instantaneous intensity.

This paper is organized as follows: section 2 contains some preliminary material and a brief discussion on the uniqueness of the spatial coherence properties for a given intensity distribution. In section 3 the proposed algorithm is presented and in section 4 the numerical simulations carried out to evaluate the proposed algorithm are described. Finally, conclusions are given in section 5.

2. Theory

Spatial coherence properties describe the randomness of a field over spatial coordinates, usually in directions perpendicular to the main propagation direction of the optical field. In this work we will limit our scope of interest to propagating fields in the scalar paraxial approximation. Further, we classify the retrieval problem as one-dimensional for fields that are functions of only one transverse coordinate, and as two-dimensional for fields that are functions of both transverse coordinates. The two-point mutual coherence function of the field in time and space can be written as

where the brackets denote a time (ensemble) average. Here u denotes the complex amplitude of the quasimonochromatic field in the slowly-varying envelope formulation, u ^* is its complex conjugate, and ρ = (x,y) denotes spatial coordinates in a plane perpendicular to the direction of the propagating field, the latter which we take to be along the z-axis. For a temporally stationary field, the amount of mutual temporal coherence is solely dependent on the time difference τ = t ₂ - t ₁. Recalling the Wiener-Khintchine theorem, the cross-spectral density function for this field is written as [12],

The cross-spectral density function is useful since it expresses the spatial coherence as a function of oscillation frequency and gives us a tool to formally separate temporal and spatial coherence. We will assume that the detailed temporal coherence properties are not important, so that S is essentially independent of v within the narrow-band frequency range of the quasimono-chromatic field. In this case the function S(ρ₁,ρ₂) is often called the (spatial) mutual coherence function. This is the function we want to calculate with our retrieval algorithm.

2.1. Uniqueness of the spatial coherence properties for a certain intensity distribution

Since our ultimate goal is to recreate the mutual coherence function of the field in a certain transverse plane from only intensity distribution data, it is of interest to first elaborate on the theoretical problem of whether the intensity distribution of a field corresponds to a unique mutual coherence function. If this is not the case, we would hope that an efficient numerical algorithm would still provide us with a solution, although it would be one out of many possible solutions.

The problem of uniqueness of fields specified on continuous spatial coordinates has been treated by several authors. We are using sampled fields specified only in discrete positions so their results are not strictly valid. We can anticipate, though, that situations where uniqueness is established for continuous fields are beneficial also for sampled fields, and vice versa. This will also be shown in one of the numerical simulations. Returning to the case of continuous fields, in [13] it is stated that for the one-dimensional case the intensity distribution in the entire space does indeed uniquely determine the mutual coherence function. On the other hand, for a two-dimensional field the intensity in the entire space is not sufficient to uniquely determine the mutual coherence function. As an example, it is possible for a conventional Gaussian-Schell beam [12] to have the same intensity distribution in the entire space as a so-called twisted Gaussian-Schell beam [14], although their spatial coherence properties are significantly different in character.

Fortunately, as pointed out by Dragoman [15], the two-dimensional case can in practice be made unique by inserting an optical element lacking rotational symmetry, such as a cylinder lens, in the beam path and using the intensity distribution data both from the space in front of the inserted element and behind, where the insertion of the element has changed the intensity distribution. An example of this technique will be shown in Section 4.2 in which the retrieval algorithm is used for a two-dimensional problem.

3. The proposed retrieval algorithm

As discussed in [16], a pulsed, partially coherent field can be represented as an ensemble of a limited number of one- or two-dimensional functions (depending on the dimensionality of the optical fields in the beam cross section). For our purposes, each function in the ensemble can be interpreted as a fully coherent optical field, which is entirely incoherent with respect to any of the other fields. The longer the duration of the pulse, the larger the number of coherent fields needed to accurately represent the partially coherent pulse. In this work we are dealing with non-pulsed, continuous-wave fields and we would thus need an infinite number of coherent fields in the representation of the partially coherent field. Since, obviously, we can only use a finite number of coherent fields, the accuracy of the recovered spatial coherence properties will, as we shall see, depend on the number of coherent fields that we use in the ensemble. Thus, in a transverse plane located at longitudinal coordinate z, the partially coherent field can be represented by an ensemble of coherent, but mutually incoherent, fields

Evidently, the partially coherent field is represented by N coherent fields, U_i(ρ, z), i = 1,…,N. Further, ρ is the transverse coordinates in the plane. Since the fields are mutually incoherent, there is no correlation between different fields so that the intensity distribution of the representation becomes

and the spatial mutual coherence function in the transverse plane located at z

It is now straightforward to apply a projection based coherence retrieval algorithm [17] to this representation. In our numerical examples we use a numerical technique similar to that of the Gerchberg-Saxton algorithm [10]. Other phase retrieval algorithms for coherent fields, such as the Yang-Gu algorithm which can account for energy loss, might also be possible to use [11]. Moreover, in our method the Gerchberg-Saxton-like enforcement of similarity between calculated and measured (reference) fields is performed in an avalanche approach going from one plane directly to its neighbor. A different stepping scheme between the planes would also be possible. In such a case a phase retrieval method which can efficiently account for inserted optical elements between the planes of propagation may be particularly beneficial [18], since this situation is encountered when we break the degeneracy of two-dimensional problems with, e.g., an inserted cylinder lens.

The proposed algorithm requires that the intensity distribution, I_ref (ρ, z), is known in a number, N_t, of transverse planes located at z = z_n, n = 1,…,N_t. The algorithm is based on repeated projection of the representation of the partially coherent field, in one plane after another, onto a new representation having an intensity distribution that is equal to the known intensity distribution in that plane. The amplitude of every coherent field in the representation is scaled by the same real positive multiplicative function. After the projection in one plane, each coherent field is individually propagated to the next transverse plane where the intensity distribution is known, and a new projection is performed. To be more specific, starting from the first operation after having propagated all coherent fields to a new plane, the proposed algorithm can be expressed as the following sequence of basic operations,

where P is the propagation operator for a coherent field from the plane located in z = z_n to a plane located in z = z _n+1. Either the two-step method, presented in [16], or the conventional angular spectrum propagation method, presented for instance in [19], can be used to conveniently perform the propagation. The index ⁺ indicates that the coherent fields are corrected, i.e. projected, so that the intensity distribution of the field representation in the plane equals the known reference intensity distribution in each lateral position ρ. Further,U_i(ρ, z_n) is the distribution in the plane at z = z_n of the i:th coherent field, and I(ρ, z_n) is the intensity distribution of the entire field representation, as given by Eq. 4. The expression in Eq. 6 defines the amplitude correction factor that is required in the current transverse plane in order for the intensity of the representation to equal the known intensity in every lateral position. Equation 7 shows how this factor is multiplied to each individual coherent field in the representation to produce the corrected coherent field. By doing so the intensity of the field representation containing the modified coherent fields equals the known intensity in the plane, as desired. Finally, Eq. 8 describes how each coherent field is then individually propagated to the next transverse plane, after which the field operations continue from Eq. 6 for the fields in the plane at z = z _n+1. The operations of the algorithm are illustrated in Fig. 1. Since the algorithm is iterative, there must be a starting distribution for all coherent fields in the plane from which the retrieval algorithm begins. In the examples in section 4, each of the coherent fields was generated from a unique random field as described in [16], so that each coherent field had a rapid spatial variation associated with a field with a much smaller coherence length, i.e. much less spatially coherent, than the actual, and therefore also the finally retrieved, fields. All the coherent starting fields were thus different in their details, but corresponding to a spatially partially coherent field with the same statistical properties.

When the algorithm has converged, one finally obtains the spatial mutual coherence function in a desired plane at z = z_n by applying Eq. 5, using the ensemble of coherent fields that was obtained in this plane, i.e., {U ₁(ρ, z_n),U ₂(ρ, z_n),…,U_N(ρ, z_n)}.

There are distinct computational advantages by using an ensemble of coherent fields rather than directly trying to use the mutual coherence function. First, the number of dimensions of the mutual coherence function is twice that of the optical fields, which tends to exhaust the computer memory if the coherence function is to be sampled with an acceptable accuracy. Second, the implementation of the ensemble of coherent fields is intuitive and straightforward and lends itself to parallelization, if one would wish to speed up the retrieval by using multiple processors or computers.

 

Fig. 1. Illustration of the operations in the coherence retrieval algorithm. In this example the intensity of the beam is measured in four transverse planes at longitudinal positions z = {z ₁,…,z ₄}. The ensemble of fields, representing the partially coherent field, is propagated one coherent field, U_i, at a time from one plane to the next, which is represented by the operator P. Each coherent field is then adjusted so that the intensity distribution of the ensemble at the new plane corresponds to the measured intensity distribution. This operation is denoted by⁺.

Download Full Size | PDF

4. Evaluation of the retrieval algorithm

To verify the efficiency of the proposed algorithm we conducted some numerical experiments. From simulated intensity measurements we strive to retrieve the coherence function. The simulated intensity distribution in a number of transverse planes was obtained analytically since analytical forms for the intensity distribution of the propagated field exist for all the examples of partially coherent fields in this work. Alternatively, one could calculate the simulated intensity distribution numerically by using the method presented in [16].

4.1. Retrieval for a one-dimensional field: the Gaussian-Schell beam

In the one-dimensional case the spatial coherence properties are uniquely determined by the intensity distribution. Thus, knowledge of the intensity distribution in multiple transverse planes (the intensity along one coordinate, say x, in the planes would be sufficient, since the field is one-dimensional) in combination with an efficient retrieval algorithm should enable us to calculate the spatial coherence properties quite accurately.

In the demonstration, we take as our example the intensity distribution of a partially coherent so-called Gaussian-Schell beam, with the goal to retrieve the spatial coherence function. The Gaussian-Schell beam has a Gaussian intensity profile and a parabolic average phase distribution of the corresponding instantaneous field. The coherence properties are quasistationary, which means that the degree of coherence between any two lateral positions is primarily dependent only on their separation, ∣x ² - x ¹∣. For the Gaussian-Schell beam this dependency is also described by a Gaussian function. The intensity distribution of the Gaussian-Schell beam at any point in space can be calculated directly from an analytical expression [12]. Moreover, the mutual coherence function in any transverse plane can be described analytically as

where the first two exponential factors relate to the Gaussian intensity of the field. The third factor is the Gaussian dependence of the coherence between two lateral positions, and the last factor is caused by the average parabolic phase curvature. In this expression, w is the 1/e ²-radius of the intensity profile, 2σ_g is the 1/e ²-radius of the two-point correlation function, and f is the radius of the phase curvature. Of course, since the problem under study is one-dimensional, the fields and the mutual coherence function are invariant along the y-axis. This means, naturally, that all fields are sampled only along the x-direction, and x ₁ and x ₂ denote the x-coordinates for two positions between which the mutual coherence is calculated. However, by convention we use the term radius for w and σ_g although the term distance would perhaps be less confusing in the one-dimensional case. In the simulations conducted the following numerical values were used for the parameters: w = 2.4∙10³λ, σ_g = 2.8∙10³λ, and f = 1.6∙10⁷λ, where λ is the center wavelength of the quasimonochromatic radiation.

To appreciate the convergence of the proposed algorithm we further introduce an error measure, E(S_calc,S_ref) that indicates the normalized distance between the mutual coherence function that is calculated with the retrieval algorithm, S_calc(x ₁,x ₂), and the known solution, S_ref (x ₁,x ₂) from Eq. 9 with the mentioned numerical values of the parameters. The error measure is defined as

A small value of E(S_calc,S_ref) indicates that the two functions, S_calc and S_ref, are similar and thus that the retrieval was successful. The error measure is normalized in the sense that E(S_calc,S_ref) ∈ [0,1] and equals zero only if S_calc and S_ref are identical.

The analysis starts by analytically calculating the intensity distribution in four planes, z = {z ₁, …,Z ₄} = {0, ⅓, ⅔, 1} ∙ 8 ∙ 10⁶λ. Thereafter the analytically calculated intensity distributions are used to retrieve the mutual coherence function S(x ₁,x ₂) in one of the planes (the one with z = z ₁ = 0 in our example) with the described retrieval algorithm. In this investigation the retrieval algorithm was executed for a sufficient number of iterations for it to converge. It was found that the accuracy of the retrieved spatial coherence function depends on several factors. First an analysis of the accuracy as a function of the number of fields, N, included in the ensemble was performed. Figure 2 shows the error measure, E(S_calc,S_ref), for the retrieved coherence function, S_calc(x ₁,x ₂), in the plane z = z ₁ as a function of TV. It is evident that the accuracy increases with the number of fields in the ensemble, although a fair agreement can be obtained with quite a low number of fields (< 100). Second, an analysis of the convergence properties of the retrieval algorithm was conducted. Figure 3 shows the error measure, E(S_calc,S_ref), of the retrieved coherence function as a function of the number of iterations, where one iteration constitutes a full cycle of going through all transverse planes once and a backpropagation of each coherent field from the final plane at z = z_Nt to the starting plane at z = z ₁. In this numerical simulation the number of fields, N, included in the ensemble was large enough not to affect the results of the simulations; see Fig. 2 for information on the minimum number of fields needed. Three different cases were analyzed for the same beam, but where the transverse planes were located differently before and after the waist of the focused beam. The separation of the planes was uniform for each of the three cases. In Fig. 3, the case shown in the top inset is for the planes at the previously mentioned locations z = {z ₁,…,z ₄}. The simulations indicate that the convergence rate may be considerably slower for an unfortunate choice of the positions of the transverse planes. Intuitively, these results should not be too surprising, as the retrieval is based on how the coherence properties influence the propagation behavior. For this reason the planes should be separated far enough for the propagation to have considerably changed the intensity distribution. The focal region of a beam displays a particularly profound change of the intensity distribution over a rather limited propagation distance. Inserting transverse planes before and after the beam waist should thus facilitate the extraction of coherence information and reduce the required number of iterations in the retrieval.

 

Fig. 2. The residual error in the retrieval of the mutual coherence function as a function of the number of fields, N, in the ensemble representing the partially coherent field.

Download Full Size | PDF

 

Fig. 3. The residual error in the retrieval of the mutual coherence function as a function of the number of iterations in the retrieval algorithm. The three different cases are simulated with the same beam, but the locations of the four transverse planes from which the intensity information is collected are chosen differently, as indicated in the inserts.

Download Full Size | PDF

4.2. Retrieval for a two-dimensional field: the twisted Gaussian-Schell beam

As mentioned, the two-dimensional case fundamentally differs from the one-dimensional case since in two dimensions the spatial coherence properties are not uniquely determined by the intensity distribution. Hence it cannot be expected that any coherence retrieval algorithm could be efficient without further measures. However, as discussed in section 2.1, if the beam is manipulated by inserting a non-rotationally symmetric optical element in its path, the uniqueness can be effectively restored by measuring intensity distributions both in front of and behind the inserted element. In the retrieval algorithm, the only thing that is changed by the insertion of the element is the propagation of each of the coherent fields from transverse plane n to n+1, where these two planes are the ones between which the element is inserted. For this propagation, the propagation operator in Eq. 8 is modified to a sequence of three operations so that

where z_element is the position of the inserted element and P_zn→z_element and P_zelement→z _n+1 are the same free space coherent propagation operators as the one in Eq. 8, only here the propagation is from the plane at z = z_n to the plane at z = z_element, and from that plane to the plane at z = z _n+1, respectively. The P_element operator is the spatial complex transmission function of the inserted element—for a lens it is thus simply a parabolically varying phase function.

To demonstrate the capabilities of the retrieval algorithm for a two-dimensional intensity distribution a so-called twisted Gaussian-Schell beam is used. A twisted Gaussian-Schell beam has a Gaussian intensity profile and a parabolic average phase curvature just like a (non-twisted) Gaussian-Schell beam. The difference is that there is an extra, non-stationary, phase factor in the mutual coherence function. In a transverse plane, the twisted Gaussian-Schell beam has a spatial mutual coherence function given by [14],

where the last exponential factor is the extra, non-stationary, phase factor that distinguishes the twisted Gaussian-Schell beam from the quasi-stationary Gaussian-Schell beam. The beam parameters that appeared in Eq. 9 have the same meaning in Eq. 12, only now the geometry under consideration is two-dimensional, so that ρ₁ = (x ₁,y ₁) and ρ₂ = (x ₂,y ₂) are used for the coordinates of two points in the transverse plane. The parameter u is the twist parameter, for which the value u = 2.0∙10^-8λ^-1 is assigned in our simulations. The other beam parameters have the same value as for the one-dimensional beam.

The effect of the non-stationary phase factor of the twisted Gaussian-Schell beam is to rotate the intensity distribution about the optical axis upon propagation. Naturally, if the beam intensity profile is circularly symmetric, as for our beam whose coherence function is given by Eq. 12, this effect is not visible. For such a beam there exists a conventional, non-twisted, Gaussian-Schell beam that has exactly the same intensity distribution in the entire space. Directly measuring the intensity distribution for either of these two beams and applying our retrieval algorithm, one could expect the algorithm to converge to a field that can be any combination of these two beam types and, possibly, any other beam type with the same intensity distribution. As mentioned, insertion of a cylinder lens removes the ambiguity of the problem, since the lens produces an elliptic cross-sectional intensity distribution in the transverse planes after the lens. In these planes the beam rotation of the twisted beam is easily detected and therefore the retrieval algorithm can now distinguish between the twisted and the conventional Gaussian-Schell beam. We illustrate this fact by performing simulations for a twisted Gaussian-Schell beam, both without any symmetry-breaking elements inserted in the beam path and with an inserted cylinder lens.

 

Fig. 4. The residual error in the retrieval of the mutual coherence function as a function of the number of iterations for the case of a two-dimensional beam. The same twisted Gaussian-Schell beam is considered in two different set ups, shown in the two inserts, with the beam propagating from left to right. The top insert shows the intensity distribution generated by the twisted beam upon propagation through a rotationally symmetric lens. The retrieval algorithm converges but not to the correct mutual coherence function because of the multitude of possible solutions. The bottom insert illustrates how the beam propagates after the insertion of a cylinder lens with different focal lengths along the x- and y-axes. Above each beam-path diagram is indicated the intensity distribution in the different planes. The non-rotationally symmetric cylinder lens produces an elliptical intensity distribution so that the rotation also becomes visible, which restores the uniqueness of the problem. As the diagram shows, the retrieved mutual coherence function now tends toward the correct solution.

Download Full Size | PDF

To evaluate the algorithm convergence the error measure, E(S_calc,S_ref), is again calculated to compare the retrieved spatial coherence function, S_calc, with S_ref. The error measure as a function of the number of iterations is shown in Fig. 4. As expected, when we use an ordinary rotationally symmetric lens the retrieved spatial coherence function does not converge toward our reference S_ref. Although the algorithm itself converges it is unable to find the desired solution because the problem is not unique. Clearly, as explained, inserting the non-rotationally symmetric cylinder lens can restore the uniqueness which enables the retrieved spatial coherence function to converge toward the actual spatial coherence function.

4.3. Evaluation of the sensitivity to noise in intensity data and measurement setup inaccuracy

In a practical retrieval of the coherence properties of a field using the proposed algorithm, the first step is the measurement of the intensity in a number of transverse planes. The measurement process inevitably produces some degree of noise contamination. It is also possible that there exists an uncertainty with regard to the position of the planes in which the intensity is measured. Both these error sources force the problem toward a state where no physical optical field can realize the measured intensity distribution. To study the behavior of the proposed retrieval algorithm under such conditions, a number of numerical simulations were performed.

 

Fig. 5. Performance of the retrieval algorithm with noisy measured intensity data. The residual error in the retrieval of the mutual coherence function is shown as function of the intensity noise root-mean-square level expressed in percent of the noise-free intensity. The plot shows the average behavior of several simulations, each generated with different distortions of the intensity data.

Download Full Size | PDF

The sensitivity to noise was simulated by again considering the Gaussian-Schell beam from section 4.1. Before the intensity data was fed into the algorithm it was intentionally distorted by a varying degree of noise. In Fig. 5 the error measure is calculated for the spatial coherence function retrieved from noisy intensity data, as a function of the root-mean-square of the noise in the intensity data. From the numerical simulations it is clear that the algorithm can handle quite noisy intensity data, and that the error in the retrieved coherence function gracefully decreases with lower noise levels.

 

Fig. 6. Performance of the retrieval algorithm when the planes where the transverse intensity data is measured are displaced compared to their locations assumed in the algorithm. The residual error in the retrieved mutual coherence function is shown as a function of the root-mean-square error, expressed in percent of the nominal distance between plane z ₁ and z ₄, of the longitudinal position of the transverse intensity planes. The plot shows the average behavior of several simulations, each generated with a different error for the location of each plane.

Download Full Size | PDF

To evaluate the impact of an error in the location of the transverse planes where the intensity is measured, a numerical simulation was conducted of the Gaussian-Schell beam from section 4.1. In Fig. 6 the error measure for the retrieved spatial coherence function is displayed as a function of the root-mean-square error of the position of the planes, along the direction of propagation. Once again the numerical simulation indicates that the algorithm is relatively stable and can handle input conditions that do not precisely represent any physically realizable optical field.

5. Conclusions

In this paper we have presented an algorithm for the retrieval of the spatial mutual coherence function of the optical field in a beam from intensity measurements in a moderate number of transverse planes. The proposed numerical algorithm is iterative and straightforward to implement. We show with numerical simulations that the proposed algorithm can robustly handle one-dimensional fields. For these fields there is a one-to-one correspondence between the intensity distribution and the coherence function. Furthermore, we show that the proposed method can also accurately retrieve the spatial mutual coherence function for two-dimensional fields, if a one-to-one correspondence is created between the intensity distribution and the coherence function. This can be achieved by inserting into the beam an optical element that lacks rotational symmetry and measuring the intensity distribution in transverse planes located both before and after the inserted element. Simple numerical illustrations of the ability of the algorithm to handle noisy or deterministic erroneous data were also included.