A comprehensive theoretical framework for the inverse design of correlation induced effects with optical beams is introduced. Correlation induced effects are able to modify the intensity distribution of an optical beam drastically via effects such as correlation induced splitting, focusing, and shifting. The inverse design steps are given analytically, which allows the analysis of several related experiments. Finally, an algorithm for more complex numerical inverse design is overviewed and demonstrated.

Topics
MATLAB, Weight function, Functions and mappings, Fresnel diffraction, Fraunhofer diffraction, Geometrical optics, Coherence theory, Laser physics, Wigner functions

Research on the coherence of light has traditionally focused on common (naturally occurring) low coherence light sources or highly coherent laboratory equipment. This has undoubtedly led to several important advances. For example, it was shown by Wolf that the spectrum of the light radiated by stars should not change on propagation,1 which is a key enabling factor for astronomical spectroscopy. Additionally, the considerations of Townes and Schawlow paved the way for the laser,2 triggering a revolution in optical sciences.

Within the past few decades, the area of partial coherence has gained increased importance due to the unique properties of partially coherent light. The first investigations of the unintuitive properties of partially coherent light were made in 1987 with the discovery of the Wolf effect3,4 and Bessel correlated fields.5 The Wolf effect was quickly confirmed experimentally,6 which sparked some interest in partially coherent fields. On the other hand, the inverse problem of determining the structure of a source from measurements of the field’s radiation pattern was found, in general, not to admit a unique solution.7 Given a priori knowledge, attempts to identify the source have been made with the delta correlated fields,8 the quasi-homogeneous sources,9,10 the Gaussian Schell-model beams,11,12 the coherent-mode representation,13 as well as through scattering media.14 

Interest in partial coherence was renewed with the discovery of correlation induced focusing15–17 and splitting18,19 of light beams, among many other correlation induced effects.20–24 Notably, novel experiments have validated several of the theoretical predictions for the properties of nonconventional correlations.18,23,24 Fields with novel correlation functions usually exhibit unique propagation properties; for example, nonuniformly correlated light beams have a shifted focus on propagation.15 It is important to note that although these are reminiscent of nonlinear effects, no material properties are involved. That is all correlation induced effects also occur upon propagation through a vacuum.

However, the process of discovering new correlation induced effects has been somewhat cumbersome so far, as it normally involves the introduction of a genuine representation25 for a correlation function, its analytical propagation, and the investigation of whether there are new effects due to the particular correlation function used. Although there is nothing wrong with this approach per se, it is quite inefficient. There are infinitely many correlation functions, and going through them one by one in hopes of discovering something new is not feasible. Still, it is desirable to investigate beams that produce correlation induced effects since they can often be superior with respect to completely coherent beams. For example, they are resistant to signal degradation when propagating through turbulent media.26–30 Such a forward design approach is, therefore, fundamentally limited.

In the present work, we introduce a comprehensive framework for the inverse design of correlation functions for experimentally realizable, on-demand effects. In other words, we start with a desired correlation induced effect and examine which type of coherence function can produce such an effect. We find the analytical forms for a variety of effects and compare these to some recent experimental efforts. Furthermore, we show the potential of correlation induced effects in problems such as inverse design.

The present work is structured as follows: we begin by reviewing the propagation of partially coherent beams with the angular spectrum approach at the beginning of Sec. II. We then move on to the Fresnel and Fraunhofer approximations for partially coherent light and derive propagation equations based on linear transformations of the corresponding Wigner function. In Sec. III, we briefly consider the genuine representation in the context of our propagation function, after which we show some particularly important analytical expressions for inverse design in Sec. IV. Finally, we showcase the wider capabilities of our approach with a numerical design algorithm in Sec. V before presenting conclusions in Sec. VI.

To begin with, we need convenient propagation formulas for the correlation functions. The cross spectral density (CSD) at the initial plane z = 0 is defined as31 
W(ρ1′,ρ2′)=⟨E*(ρ1′)E(ρ2′)⟩,
(1)
where E(ρ) is a member of a statistical ensemble of monochromatic realizations, asterisk denotes complex conjugation, ρj′=(xj′,yj′)⁠, with j = 1, 2, contains the transverse spatial coordinates at the initial plane, and the frequency ω is left implicit for brevity. The spatial signal forms a Fourier transform pair with the angular domain signal,
E(κ)=1(2π)2∫−∞∞E(ρ′)exp(−iκ⋅ρ′)d2ρ′,
(2)
which is also known as the angular spectrum. Here, κ = (kx, ky) contains the transverse wave vector components. Therefore, the CSD can be transformed to the angular correlation function (ACF) with the generalized Wiener–Khintchine theorem,
A(κ1,κ2)=⟨E*(κ1)E(κ2)⟩=1(2π)4∬−∞∞W(ρ1′,ρ2′)expiκ1⋅ρ1′−κ2⋅ρ2′×d2ρ1′d2ρ2′.
(3)
The ACF can be used to propagate the field to any plane z by imposing a propagation phase of the form
A(κ1,κ2;z)=A(κ1,κ2)expikz2−kz1*z,
(4)
where kzj=kj2−κj2 is the wave vector component toward the propagation direction, kj = |kj|, kj being the wave vector in free space at a frequency ω. In general, kzj can be complex, so evanescent waves are implicitly accounted for. The spatial domain correlation function can be retrieved with
W(ρ1,ρ2,z)=∬−∞∞A(κ1,κ2)expikz2−kz1*z×exp−iκ1⋅ρ1−κ2⋅ρ2d2κ1d2κ2,
(5)
where ρj = (xj, yj) contains the transverse coordinates at the plane z > 0. Note that the initial plane coordinates were defined under Eq. (1) with the primed symbols.
Combining Eqs. (3)–(5), we can write the propagated CSD in the explicit form
W(ρ1,ρ2,z)=1(2π)4∬−∞∞∬−∞∞W(ρ1′,ρ2′)×expiκ1⋅ρ1′−κ2⋅ρ2′d2ρ1′d2ρ2′×expikz2−kz1*z×exp−iκ1⋅ρ1−κ2⋅ρ2d2κ1d2κ2.
(6)
Furthermore, according to the Fubini–Tonelli theorem,32,33 we can rearrange the order of integrals if the function to be integrated fulfills the requirement,
∬|f(x,y)|dxdy<∞,
(7)
which is always fulfilled with genuine optical correlation functions due to the requirement of finite energy. Therefore, we can cast the propagated CSD into the simplified form
W(ρ1,ρ2,z)=1(2π)4⨌−∞∞W(ρ1′,ρ2′)×exp−iκ1⋅(ρ1−ρ1′)−κ2⋅(ρ2−ρ2′)×expikz2−kz1*zd2κ1d2κ2d2ρ1′d2ρ2′.
(8)
We emphasize that these propagation equations hold for any quasi-monochromatic statistically stationary field of any state of spatial coherence described by W(ρ1′,ρ2′) at the initial plane. The employed notation is shown in Fig. 1. However, evaluating four two-dimensional integrals is tedious. We can further simplify the equations by employing a few common approximations.
FIG. 1.
Basic geometry of a propagation problem. The intensity at the initial plane (left, primed coordinates) is transformed to the intensity at a propagated plane (right); the dotted line shows a constant intensity contour. The angle θ corresponds to the divergence angle when the propagation has been taken to the far-zone, while θ′ is significant in the near-zone.
View largeDownload slide

Basic geometry of a propagation problem. The intensity at the initial plane (left, primed coordinates) is transformed to the intensity at a propagated plane (right); the dotted line shows a constant intensity contour. The angle θ corresponds to the divergence angle when the propagation has been taken to the far-zone, while θ′ is significant in the near-zone.

FIG. 1.
Basic geometry of a propagation problem. The intensity at the initial plane (left, primed coordinates) is transformed to the intensity at a propagated plane (right); the dotted line shows a constant intensity contour. The angle θ corresponds to the divergence angle when the propagation has been taken to the far-zone, while θ′ is significant in the near-zone.
View largeDownload slide

Basic geometry of a propagation problem. The intensity at the initial plane (left, primed coordinates) is transformed to the intensity at a propagated plane (right); the dotted line shows a constant intensity contour. The angle θ corresponds to the divergence angle when the propagation has been taken to the far-zone, while θ′ is significant in the near-zone.

Close modal
As a first approximation, we consider the propagation of paraxial narrowband fields such that the first term of the Taylor expansion of
kzj≈k−κj22k,
(9)
is sufficient to model the propagation properties of such fields. Here k = |k| = ω0/c is the magnitude of the wavevector evaluated at the center frequency of the field, ω0, and c is the speed of light in vacuum. Inputting this into Eq. (8) and integrating with respect to κ1 and κ2 yields the usual Fresnel propagation integral for partially coherent light,
W(ρ1,ρ2,z)=k2πz2∬−∞∞W(ρ1′,ρ2′)×exp−ik(ρ1−ρ1′)2−(ρ2−ρ2′)22zd2ρ1′d2ρ2′.
(10)
If we now expand the squares in the exponential and introduce the average ρ = (ρ1 + ρ2)/2 and Δρ = ρ2 − ρ1, difference coordinates, we can rewrite the Fresnel propagation equation in the following, inspiring, form:
W(ρ,Δρ,z)=k2πz2⁡expikρΔρz∬−∞∞W(ρ′,Δρ′)×exp−ikρΔρ′−ρ′Δρ′+ρ′Δρzd2ρ′d2Δρ′.
(11)
where W(ρ′, Δρ′) = W(ρ′ − Δρ′/2, ρ′ + Δρ′/2) is the rotated CSD. This form is reminiscent of a Wiener–Khintchine relation, except that there is a term that couples the average and difference coordinates together.
The above result can be further simplified by collecting the Δρ′ factors in the exponent and introducing the Wigner representation of the CSD as34 
Wg(ρ′,κ)=12π2∫−∞∞W(ρ′,Δρ′)exp(−iκ⋅Δρ′)d2Δρ′,
(12)
so that we can write the Fresnel integral in the following, final form:
W(ρ,Δρ,z)=kz2⁡expikρΔρz×∫−∞∞Wgρ′,kρ−ρ′zexp−ikρ′Δρzd2ρ′.
(13)
This is the main result of our work. We have reduced the Fresnel propagation integral to a single two-dimensional Fourier transform over a deformed Wigner distribution. The deformation due to propagation is governed by a linear transformation of the form
ρ′k(ρ−ρ′)/z=Tρ′ρ,
(14)
where
T=10−k/zk/z,
(15)
is the transformation matrix due to propagation. The matrix T has two contributions of equal weight: the on-diagonal element corresponds to stretching, and the off-diagonal element causes shearing.
Let us consider the intensity distributions observed at different propagation lengths. The intensity is found along the diagonal of the CSD by setting ρ1 = ρ2 = ρ, which is equivalent to setting Δρ = 0 in the average and difference coordinates, W(ρ, 0, z) = S(ρ, z). This yields
S(ρ,z)=kz2∫−∞∞Wgρ′,kρ−ρ′zd2ρ′.
(16)
That is, we can calculate the intensity distribution of the propagated partially coherent field by simply integrating its deformed Wigner distribution. This is a particularly useful formula for the design of more complex correlation induced effects, which will be investigated in more detail in Sec. V.

The Fraunhofer approximation is a special case of the Fresnel formula, where the beam has propagated to the far-zone. This propagation length has customarily been evaluated as a distance after which the beam no longer changes shape, and mathematically, it has been found to be a length that is much greater than the Rayleigh length for a partially coherent beam, given by zR=kw02β/2⁠.35–37 Here, w0 is the initial width of the beam, and β ∈ [0, 1] is a parameter that is related to the overall degree of coherence (DOC) of the field. Since it is limited within the interval between 0 and 1, it is enough to consider the Rayleigh length of a completely coherent field when evaluating whether a beam is in the far-zone.

To move from the Fresnel to the Fraunhofer approximation, we assume that the term ρ′Δρ′ in Eq. (11) can be neglected, as for propagation distances z ≫ zR, it becomes insignificant with respect to the other length scales at play. With this approximation, the Fresnel integral can be evaluated directly in terms of the Fourier transform of the CSD, i.e.,
W(ρ,Δρ,z)=k2πz2⁡expikzρΔρAkzρ,kzΔρ,
(17)
where A(κ1, κ2) is the angular spectrum given by Eq. (3).

Therefore, the far-zone correlation function has the same functional form as the ACF. Note that the ACF is not modified upon propagation and, therefore, has the same form at the far-zone as in the initial plane (apart from a change of scale). The Fourier transform relationship between the initial plane CSD and the ACF, therefore, limits the possible correlation induced effects obtained in far-zone propagation. However, a lot of freedom is retained since changes to the phase of the initial plane CSD are not visible in its intensity distribution.

Moreover, in the Fraunhofer regime, we can write the intensity distribution in three equivalent ways, as in
S(ρ,z)=k2πz2∬−∞∞W(ρ′,Δρ′)exp−ikzρΔρ′d2ρ′d2Δρ′=kz2∫−∞∞Wgρ′,kρ/zd2ρ′=k2πz2Akρ/z,0.
(18)
Therefore, it is evident that the far-zone intensity distribution is entirely determined by the diagonal of the ACF, A(κ, 0), which is influenced by the form of the initial plane CSD. The second equality is again marginal over a Wigner distribution, with only stretching and no shear. In other words, Eq. (16) asymptotically approaches the functional form of Eq. (18) as the beam propagates toward the far-zone. The Wigner distribution is an object that connects the intensity distribution at the source plane to any other plane of propagation, and if we take the other marginal instead, we obtain the source plane intensity,
S(ρ′,0)∝∫−∞∞Wgρ′,kρ/zd2ρ.
(19)

The Wigner distribution allows for an alternative path for computing the dependence between the initial plane CSD and the far-zone intensity distributions. We emphasize that these propagation equations are quite general and hold for paraxial and quasi-monochromatic statistically stationary fields of any state of spatial coherence.

Apart from the propagation formulas, another important ingredient for the inverse design is the genuine representation. That is, for a correlation function to be physically realizable, it must admit a genuine representation, which can be written as25 
W(ρ1′,ρ2′)=∫−∞∞p(v)H*(ρ1′,v)H(ρ2′,v)d2v.
(20)
Here, p(⁠v⁠) is a non-negative weight function, and H(ρ′, v⁠) represents an arbitrary kernel, usually written as a Fourier kernel
H(ρ′,v)=τ(ρ′)exp−i2πf(ρ′)⋅v,
(21)
where f(ρ′) is a function that depends only on the initial plane position ρ′, and τ(ρ′) is in general a complex amplitude function.
With these choices, we can evaluate Eq. (20) to obtain W(ρ1′,ρ2′)=C(ρ1′,ρ2′)g(Δf)⁠, where
C(ρ1′,ρ2′)=τ*(ρ1′)τ(ρ2′),
(22)
Δf=f(ρ2′)−f(ρ1′)⁠, and
g(Δf)=∫−∞∞p(v)exp−i2πΔf⋅vd2v,
(23)
is proportional to the spectral DOC. Notice how all intensity information and constant phase factors are contained in C(ρ1′,ρ2′)⁠, while the correlations are shaped by g(Δf).

However, both the function g(Δf) itself and its argument f cannot be chosen arbitrarily but must fulfill at least the following three conditions:

  • g(Δf) is equal to the spectral degree of coherence up to a phase term, so Re{g(Δf)} ∈ [−1, 1];

  • due to (23), the function g(Δf) must be Fourier-transformable, i.e., g:L1(R2)→C and;

  • if f(ρ1′)=f(ρ2′)⁠, then it must follow that g(Δf) = g(0) = 1, to ensure normalization of p(⁠v⁠).

Within the constraints dictated by these conditions, we can freely choose the mathematical forms of g(Δf) and f(ρ′) while still retaining compliance with the genuine representation. For example, we can choose p(v) as a Gaussian function and f as a linear function to get the usual Gaussian Schell-model correlation. If we choose f as a quadratic function, then we get a nonuniformly correlated CSD. Essentially, p(v) and f can be any linear combination of elementary functions that fulfills the conditions outlined earlier.

Similarly, the complex amplitude τ(ρ) cannot be chosen arbitrarily, but must abide by the condition of being a finite-energy-carrying function, i.e., τ(ρ)∈L2(R2)⁠, and it must be continuous and smooth so that τ ∈ C∞ to ensure compliance with Maxwell’s equations.

It is important at this point in our analysis to remark on the physical meaning of the function g(Δf). Although it was noted earlier that the correlations in the field are shaped by the function g(Δf) itself, this does not necessarily mean that g(Δf) is the spectral DOC. The latter is, in fact, defined as
μ(ρ1′,ρ2′)=W(ρ1′,ρ2′)S(ρ1′)S(ρ2′).
(24)
If we now insert the definition of the CSD in terms of the function g(Δf), i.e., W(ρ1′,ρ′2)=C(ρ1′,ρ2′)g(Δf)⁠, above, we get
μ(ρ1′,ρ2′)=C(ρ1′,ρ2′)S(ρ1′)S(ρ2′)g(Δf).
(25)
From this result, we can clearly see that μ(ρ1′,ρ2′)=g(Δf) only if C(ρ1′,ρ2′)=S(ρ1′)S(ρ2′)⁠, which implies that τ(ρ′)∈R⁠.

This result allows us to complement our claim made in Sec. II A, that the propagation of a partially coherent field amounts to a deformation of the governing Wigner distribution. In fact, the above-mentioned results clearly show that any CSD admitting a genuine representation at the initial plane will also have a genuine representation at any other plane of propagation. This is mainly due to the fact that the linear map imposed by propagation only changes the coordinates and not the genuine representation itself.

To show this explicitly, let us substitute the expression of W(ρ1′,ρ2′) in terms of the genuine representation into the Fresnel equation of Eq. (11) to obtain
W(ρ1,ρ2,z)=k2πz2∬−∞∞C(ρ1′,ρ2′)g(Δf′)×exp−ik(ρ1−ρ1′)2−(ρ2−ρ2′)22zd2ρ1′d2ρ2′,
(26)
from which we can derive the corresponding intensity distribution as
S(ρ,z)=k2πz2∬−∞∞C(ρ1′,ρ2′)g(Δf′)×exp−ik(ρ−ρ1′)2−(ρ−ρ2′)22zd2ρ1′d2ρ2′,
(27)
which is, in fact, equivalent to Eq. (16). Unfortunately, these equations do not appear to admit an analytical solution for general expressions of the functions C(ρ1′,ρ2′) and g(Δf′). To gain more insight on the physics of the above-mentioned equation, let us then consider the quasihomogeneous case and consequently assume that C(ρ1′,ρ2′)=S(ρ1′)S(ρ2′)≈S(ρ′) and g(Δf′) = μ(Δf′) depends only on the distance between coordinates (i.e., it follows the Schell-model), Δf′=Δρ′=ρ2′−ρ1′⁠, while being narrow with respect to S(ρ′). With these approximations, we can perform the integration over Δρ′ in the above-mentioned expression, which gives the following result:
S(ρ,z)=k2πz2∫−∞∞S(ρ′)pρ′−ρλzd2ρ′,
(28)
where p is the weight function introduced in Eq. (20). That is, the propagated intensity is essentially a cross-correlation between S(ρ′) and the deformed weight function p.

This is the second result of our work. If a partially coherent field is described by a certain genuine representation at the initial plane z = 0, its representation at any other plane z > 0 will still be genuine, and the intensity distribution of the field at any plane z > 0 is given by the cross-correlation between the initial intensity distribution of the field and the genuine weight function p(x), suitably deformed according to the propagation matrix T introduced in Sec. II.

We may start our design process by specifying that we have an initial plane CSD W(ρ1′,ρ2′)⁠, from which we fix only the diagonal element, i.e., the initial plane intensity distribution S(ρ′). Then we say that we want this intensity distribution to evolve upon propagation to a different functional form, S(ρ, z). That is we choose what the initial plane and propagated intensity distributions look like and try to find the CSD that allows this to happen.

For a general field, one could, at least in principle – apply the general propagation formula of Eq. (8) between the initial and propagated planes and compute the required form of W(ρ1′,ρ2′)⁠. However, this problem is underdetermined and involves tedious integrals—it is not an analytically tractable problem. Therefore, when the field is strongly diverging or three dimensional, a numerical solver would be required. This could be in the form of an iterative Fourier transform algorithm (IFTA), where one would start from an initial guess for the initial plane CSD, then propagate it a small step, update, and repeat until the field behaves in the desired way. Due to the underdetermined nature of the problem, the solutions obtained with such an algorithm may not be unique.

If the field is paraxial, then the analytical expressions simplify greatly, and one merely needs to find the marginals of a linearly transformed Wigner distribution, as in Eq. (16), and from there, reconstruct the correct CSD for the desired correlation induced effect. While computationally less expensive, even this requires a numerical solver, which we briefly discuss in Sec. V. Finally, by approximating that the field is quasihomogeneous, we obtain equations that are relatively simple to solve analytically. If one is concerned with transient effects (such as focusing or object avoidance), then Fresnel domain expressions such as Eq. (28) need to be considered. However, a large variety of correlation induced effects can be manufactured with the far-zone (Fraunhofer) expressions, as we will show next.

Let us now take the inverse Fourier transform over both sides of the first equality of Eq. (18). We can recast the Fraunhofer propagation equation so that
∫−∞∞W(ρ′,Δρ′)d2ρ′=∫−∞∞Sρ,z≫zRexpikρΔρ′zd2ρ.
(29)
If we assume again a quasihomogeneous field, we can approximate W(ρ′, Δρ′) ≈ S(ρ′)μ(Δρ′), in which case we have
μ(Δρ′)=∫−∞∞Sρ,z≫zRexpikρΔρ′zd2ρ∫−∞∞S(ρ′)d2ρ′,
(30)
which allows us to build any quasihomogeneous CSD unambiguously since the right hand side contains only terms, which we are free to choose in any way we want (as long as the quasi-homogeneity is fulfilled). As long as the far-zone intensity distribution is not pathological, the above-mentioned equation fulfills the necessary conditions for admitting a genuine representation. Furthermore, the constraint of the spectral DOC, 0 ≤ |μ(Δρ′)| ≤ 1, needs to be met, and |μ(0)| = 1 holds, which is always fulfilled when energy is conserved.

In Secs. IV A–IV C, we will look at some well-known examples that have been extensively studied and experimentally confirmed. These examples validate our approach and highlight the possibilities we have for controlling the intensity distribution at the target plane.

As the first example, we consider a beam that has a Gaussian distribution at the initial plane and otherwise retains its functional form but shifts on propagation. The spectral intensity distribution at the source plane takes the form of
S(ρ′)=exp−ρ′2w02,
(31)
while the spectral intensity in the far zone is written as
S(ρ,z)=zR2z2exp−zR2(ρ±ρ0)2w02z2,
(32)
where ρ0 is the shift parameter and zR=kw02β/2 is the Rayleigh range introduced earlier, with β=1+w02/σ02−1/2 chosen to correspond to a Gaussian Schell-model beam, with w0 and σ0 being the rms width and coherence width of the source, respectively. In the case of a quasihomogeneous field with σ0 ≪ w0, the Rayleigh range is approximated as zR = kw0σ0/2.
The source spectral DOC can now be obtained with the help of Eq. (30) as
μ(Δρ′)=exp−Δρ′2σ02∓iρ0⋅Δρ′.
(33)
In other words, the shifting of a partially coherent beam on propagation is caused by a linear phase term in the coherence function, as one would expect from the basic properties of the Fourier transform. This is entirely consistent with the theory and experiments performed in Ref. 23.

Next, we consider beams, which attain a completely new type of intensity distribution as they propagate to the far-zone. In practical terms, there are infinitely many different intensity distributions one could choose from. For the sake of this example, we look at a family of beams that have the same initial plane spectral intensity distribution as in Eq. (31) but instead attain flat-topped intensity distributions after propagation to the far-zone.

One such spectral intensity at the far zone is given by the multi-Gaussian function38,
S(ρ,z)=1C0∑m=1MMm(−1)m−1zR2mz2exp−zR2ρ2w02z2,
(34)
where C0=∑m=1MMm(−1)m−1m is the normalized factor. Taking similar steps as earlier, the corresponding source complex DOC for a quasihomogeneous field is expressed as
μ(Δρ′)=1C0∑m=1MMm(−1)m−1mexp−Δρ′2mσ02,
(35)
which is consistent with the form of DOC in Ref. 38 as well as the experiments laid out in Ref. 24. In a similar fashion, we could design initial plane correlations that produce any other physically realizable intensity distribution in the far-zone.
As the last simple example, we take beams that split into several sub-beams upon propagation. Consider the Laguerre Gaussian Schell-model beams, whose spectral intensity at the source plane is still Gaussian as earlier but splits at the far-zone, where its spectral intensity is expressed as39,
S(ρ,z)=zR2z2Ln−zR2ρ2σ02z2exp−zR2ρ2w02z2.
(36)
Again, recalling that we consider the quasihomogeneous case, the source spectral DOC is derived as
μ(Δρ′)=LnΔρ′2σ02exp−Δρ′2σ02,
(37)
which takes a form similar to the one employed in Ref. 39 and is compatible with the experiments laid out in Ref. 18. These examples show that the behavior of the beam can be deduced from a few simple Fourier type relations applied to the initial plane correlation function. In Sec. V, we will introduce a more advanced example, which is most conveniently handled via numerical techniques.

Although the analytical approach allows us to design a large variety of different types of correlation induced effects, it is inherently limited to rather simple cases. More specifically, the analytical method outlined earlier works only for transformations between the source plane and the far-zone, for beams that follow the Schell-model and are quasi-homogeneous. These are rather strict limitations, and therefore, it is desirable to consider more general design cases where the beam may not be quasi-homogeneous, the correlations might not follow the Schell-model, and there are no direct transformations between the source plane and the far-zone.

In the following, we introduce the specific numerical implementation for more general design tasks. This will allow us to later design a beam with transient effects, i.e., a beam that shows unconventional properties upon propagation but only within a limited volume. Let us consider the Fresnel domain formula, which is of the form
S(ρ,z)=kz2∫−∞∞Wga,b−cd2ρ′.
(38)
This formula describes a beam that moves from the far-zone toward the source plane as the parameter c is increased. That is, a larger shear corresponds to shorter propagation, and when c → ∞, the expression is at the source plane. This behavior can be confirmed by considering Eq. (18), where c is zero. This causes two problems: first, numerical propagation cannot handle infinities. Second, a large shear also requires a large data window, which shrinks as the beam propagates toward the far-zone. Without rescaling the transverse axes, it would appear that the beam narrows as it propagates, and proper scaling must be introduced to recover the correct behavior.
A computationally simple method to overcome this is by changing variables in Eq. (16) so that κc = k(ρ − ρ′)/z and, therefore, ρ′ = ρ − zκc/k,
S(ρ,z)=∫−∞∞Wgρ−zκc/k,κcd2κc.
(39)
In this equation, we can keep κc constant, and the propagation is entirely handled by shearing along the ρ-axis. Note that the amplitude scaling factors are absorbed into the integral in this representation, further simplifying the numerical propagation.

By employing Eq. (39), it is rather straightforward to construct a numerical algorithm for propagating beams with any state of coherence. The algorithm can start from any physically realizable CSD, which is then propagated by shearing the corresponding Wigner function. The shearing can be accomplished with a two-dimensional affine transformation or by shifting each row of the Wigner function with commands like “circshift” in Matlab. The affine transformation corresponds directly to the shear map, but it resizes the data window as the shear increases. Shifting each row accomplishes the same effect without resizing the data window and may be the preferred method, especially when memory is an issue. However, one needs to be careful that the Wigner function does not wrap onto itself when using this method, which will happen as long as the beam is propagated far enough (depending on the zero padding around the Wigner distribution).

Next, we will consider how to numerically design correlation induced effects. We consider one dimensional correlations for simplicity and employ a type of coherent mode superposition,
W(x1,x2)=∑ncnen*(x1)en(x2),
(40)
where cn is a non-negative weight function and en(x) are shifted and tilted Gaussian coherent modes, with an explicit form of
e(x)=exp−(x−xn)2w2exp(ixan),
(41)
with xn and an being the shift and tilt parameters, respectively. The modes may not be orthonormal, and thus they are not defined in the same manner as laser modes. Still, any incoherent superposition of the form of Eq. (40) yields a genuine correlation function.

One can start the design process in any plane of propagation and move into Wigner space. In the Wigner space, the shift parameter causes the elementary contribution to move along the ρ′-axis, while the tilt parameter will move the contribution along the ρ-axis. This process is graphically outlined in Fig. 2, and it gives access to the whole Wigner space. By suitably distributing mutually uncorrelated elementary modes, it is possible to add energy where it is needed.

FIG. 2.
Schematic representation of the design process in Wigner space. Two Gaussian modes are shown, with the effect of varying the shift and tilt parameters. Solid red lines go through the centroids of the modes, whereas the dashed red lines denote the origin. The blue curves on the edges of the Wigner space show the corresponding marginals.
View largeDownload slide

Schematic representation of the design process in Wigner space. Two Gaussian modes are shown, with the effect of varying the shift and tilt parameters. Solid red lines go through the centroids of the modes, whereas the dashed red lines denote the origin. The blue curves on the edges of the Wigner space show the corresponding marginals.

FIG. 2.
Schematic representation of the design process in Wigner space. Two Gaussian modes are shown, with the effect of varying the shift and tilt parameters. Solid red lines go through the centroids of the modes, whereas the dashed red lines denote the origin. The blue curves on the edges of the Wigner space show the corresponding marginals.
View largeDownload slide

Schematic representation of the design process in Wigner space. Two Gaussian modes are shown, with the effect of varying the shift and tilt parameters. Solid red lines go through the centroids of the modes, whereas the dashed red lines denote the origin. The blue curves on the edges of the Wigner space show the corresponding marginals.

Close modal

In principle, it is also possible to remove energy from the Wigner function, and it is not uncommon to have areas where the Wigner function takes on negative values. However, causing negative dips in the Wigner space by hand does not seem like an effective design method since it very easily causes the intensity to take on negative values. Therefore, the safest method is to add energy in the desired region via the shifted and tilted coherent modes.

We can now take one of two paths: either (i) design a propagated Wigner function at the plane where we want to see the effect and back-propagate it to the initial plane, or (ii) design the initial plane Wigner function while keeping in mind that propagation will shear it. In both cases, one needs to consider the distribution in terms of the tilt and shift parameters visualized in Fig. 2.

Let us say that we want to produce a novel correlation function that combines two different effects: self-splitting and self-focusing. Moreover, we want this to happen inside a limited volume, such that it is a transient effect. For this purpose, it is convenient to consider option (ii). That is, we can arrange energy in the Wigner function such that the source plane and far-zone intensities are Gaussian, while a shear along one axis (propagation) will cause the intensity distribution to evolve in an unexpected manner. This can be achieved by arranging the energy in slanted columns: as the beam propagates, the column will first turn upright (self-focusing), after which it will slant in the opposite direction. By adding more columns side-by-side, we can introduce sub-beams (self-splitting).

In Fig. 3, we present a beam that has specifically designed correlation induced properties (see supplementary material for the corresponding Matlab code). It consists of 303 shifted and tilted Gaussian modes as described by Eqs. (40) and (41). The modes form three sub-beams, each containing 101 modes, weighted with a Gaussian function. The beam starts off with a shape close to a Gaussian beam, which then splits into three self-focusing sub-beams. After further propagation, the beams recombine and approach a near-Gaussian shape in the far-zone. This example shows novel propagation behavior for partially coherent beams by combining two well-known correlation induced effects (namely, self-splitting and self-focusing). This shows that the inverse design of such effects offers great flexibility.

FIG. 3.
Beam evolution upon propagation for the inverse designed beam. (a) intensity distribution, (b) the corresponding DOC at z = 0, (c) in the focal plane at ∼z = 50, and (d) in the far-zone at z = 120. In (b)–(d), the intensity distribution at each plane is overlaid with a black dashed line.
View largeDownload slide

Beam evolution upon propagation for the inverse designed beam. (a) intensity distribution, (b) the corresponding DOC at z = 0, (c) in the focal plane at ∼z = 50, and (d) in the far-zone at z = 120. In (b)–(d), the intensity distribution at each plane is overlaid with a black dashed line.

FIG. 3.
Beam evolution upon propagation for the inverse designed beam. (a) intensity distribution, (b) the corresponding DOC at z = 0, (c) in the focal plane at ∼z = 50, and (d) in the far-zone at z = 120. In (b)–(d), the intensity distribution at each plane is overlaid with a black dashed line.
View largeDownload slide

Beam evolution upon propagation for the inverse designed beam. (a) intensity distribution, (b) the corresponding DOC at z = 0, (c) in the focal plane at ∼z = 50, and (d) in the far-zone at z = 120. In (b)–(d), the intensity distribution at each plane is overlaid with a black dashed line.

Close modal

Note that the correlation functions shown in Figs. 3(b)–3(d) do not correspond to the Schell-model, although the beam has low coherence (the CSDs are narrow along Δx direction). This presents, therefore, a novel tool for designing correlation functions.

In the present work, we have developed a comprehensive framework for inverse designing partially coherent beams, which exhibit unusual effects upon propagation. The process consists of first defining what type of effect we wish to see and then finding the form of the correlation function that can achieve this. We analytically examine three well-known examples, which have been demonstrated experimentally as well. We then move on to demonstrate the numerical implementation with a novel correlation induced effect, combining self-focusing and self-splitting. The method is fairly flexible, and it can be used to design a wide variety of correlation functions for quasi-monochromatic beams. Note that we are not limited to Schell-model sources. The method is not applicable to surface waves, highly divergent fields, or pulses in the present form.

Although we consider only propagation through free space, the method is also applicable to beams propagating in linear media with a refractive index that is not unity. This will simply affect the magnitude of the wavevector and, therefore, the distance at which these effects occur. Turbulent and/or nonlinear media would affect the design process less predictably; in some cases, the method might be applicable, whereas in others, it does not work at all. The outcome would strongly depend on the properties of the medium and the field.

This study is entirely theoretical, and a few words on the experimental realizability are in order. First of all, the analytical examples in Sec. IV concern quasihomogeneous (Schell-model) beams, which may be generated with standard tools of coherence research. That is, by first taking a low coherence beam—generated from a single mode laser that is sent through a rotating ground glass, for example—and subsequently filtering the beam with a suitable phase and/or intensity mask. On the other hand, the example shown in Sec. V does not follow the Schell-model and, therefore, its synthesis is more complex as well.

The shifted and tilted coherent mode superposition that was used to generate the example can be realized experimentally as well. For example, one can have an intensity mask with pinholes at the desired locations dictated by the values of xn. When incoherent light is incident from behind the mask, the pinholes will filter the beam such that the light emerging from each individual pinhole is completely coherent, but there are no correlations between the two different pinholes. This gives us the desired shifts, and the mutual weights of the modes can be controlled with either a secondary intensity mask or by varying the pinhole sizes. Then, the tilt an can be realized by imposing a linear phase at each pinhole. This can be performed via a 4f imaging system together with a spatial light modulator or with a lens array on top of the intensity mask (which could be produced with 3D printing, for example).

Therefore, the concepts presented here are entirely possible to realize experimentally as well. Paired with the inverse design steps outlined here, this opens up many new possibilities for designing novel correlation induced phenomena, including on-demand effects.

See supplementary material for a Matlab script used to produce Fig. 3. The code contains three parts: first part constructs the correlation function with the shifted and tilted modes, second propagates the resulting correlation function by deforming the corresponding Wigner function, while the third part produces a figure similar to Fig. 3. The deformation of the Wigner function is plotted in real time.

M.K. acknowledges useful discussions with Professor Etienne Brasselet. M.L. acknowledges the financial support of the National Natural Science Foundation of China (NSFC) (Grant No. 61805080), the Hunan Provincial Natural Science Foundation of China (Grant No. 2019JJ50366), and the China Scholarship Council (CSC). The authors acknowledge financial support from the Academy of Finland Flagship Program, Photonics Research and Innovation (PREIN), and Decisions Grant Nos. (320165) and (320166).

The authors have no conflicts to disclose.

Meilan Luo: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal). Marco Ornigotti: Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Matias Koivurova: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (lead); Supervision (equal); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

1.
E.
Wolf
,
Phys. Rev. Lett.
56
,
1370
(
1986
).
https://doi.org/10.1103/physrevlett.56.1370
Crossref
Search ADS
PubMed
 
2.
A. L.
Schawlow
 and
C. H.
Townes
,
Phys. Rev.
112
,
1940
(
1958
).
https://doi.org/10.1103/physrev.112.1940
Crossref
Search ADS
 
3.
E.
Wolf
,
Nature
326
,
363
(
1987
).
https://doi.org/10.1038/326363a0
Crossref
Search ADS
 
5.
F.
Gori
,
G.
Guattari
, and
C.
Padovani
,
Opt. Commun.
64
,
311
(
1987
).
https://doi.org/10.1016/0030-4018(87)90242-2
Crossref
Search ADS
 
6.
M. F.
Bocko
,
D. H.
Douglass
, and
R. S.
Knox
,
Phys. Rev. Lett.
58
,
2649
(
1987
).
https://doi.org/10.1103/physrevlett.58.2649
Crossref
Search ADS
PubMed
 
7.
A. J.
Devaney
,
J. Math. Phys.
20
,
1687
(
1979
).
https://doi.org/10.1063/1.524277
Crossref
Search ADS
 
8.
H. P.
Baltes
,
M.
Bertero
, and
R.
Jost
,
Inverse Scattering Problems in Optics
(
Springer
,
1980
), Vol.
20
.
Google Scholar
Crossref
Search ADS
 
9.
W. H.
Garter
 and
E.
Wolf
,
J. Opt. Soc. Am. A
2
,
1994
(
1985
).
https://doi.org/10.1364/josaa.2.001994
Crossref
Search ADS
 
10.
G.
Gbur
,
Opt. Commun.
187
,
301
(
2001
).
https://doi.org/10.1016/s0030-4018(00)01115-9
Crossref
Search ADS
 
11.
A.
Beckus
,
A.
Tamasan
,
A.
Dogariu
,
A. F.
Abouraddy
, and
G. K.
Atia
,
J. Opt. Soc. Am. A
35
,
959
(
2018
).
https://doi.org/10.1364/josaa.35.000959
Crossref
Search ADS
 
12.
D.
Kuebel
 and
T. D.
Visser
,
Opt. Lett.
45
,
1375
(
2020
).
https://doi.org/10.1364/ol.385644
Crossref
Search ADS
PubMed
 
13.
T.
Habashy
,
A. T.
Friberg
, and
E.
Wolf
,
Inverse Probl.
13
,
47
(
1997
).
https://doi.org/10.1088/0266-5611/13/1/005
Crossref
Search ADS
 
14.
M.
Koirala
,
R.
Sarma
,
H.
Cao
, and
A.
Yamilov
,
Phys. Rev. B
100
,
064203
(
2019
).
https://doi.org/10.1103/physrevb.100.064203
Crossref
Search ADS
 
15.
H.
Lajunen
 and
T.
Saastamoinen
,
Opt. Lett.
36
,
4104
(
2011
).
https://doi.org/10.1364/ol.36.004104
Crossref
Search ADS
PubMed
 
17.
Y.
Chen
 and
Y.
Cai
,
High Power Laser Sci. Eng.
4
,
e20
(
2016
).
https://doi.org/10.1017/hpl.2016.19
Crossref
Search ADS
 
18.
Y.
Chen
,
J.
Gu
,
F.
Wang
, and
Y.
Cai
,
Phys. Rev. A
91
,
013823
(
2015
).
https://doi.org/10.1103/physreva.91.013823
Crossref
Search ADS
 
19.
C.
Ding
,
M.
Koivurova
,
J.
Turunen
, and
L.
Pan
,
Phys. Rev. A
97
,
053838
(
2018
).
https://doi.org/10.1103/physreva.97.053838
Crossref
Search ADS
 
20.
J.
Turunen
 and
A. T.
Friberg
,
Prog. Opt.
54
,
1
–
88
(
2010
).
https://doi.org/10.1016/S0079-6638(10)05406-5
Crossref
Search ADS
 
21.
L.
Ma
 and
S. A.
Ponomarenko
,
Opt. Lett.
39
,
6656
(
2014
).
https://doi.org/10.1364/ol.39.006656
Crossref
Search ADS
PubMed
 
22.
O.
Korotkova
,
Opt. Lett.
39
,
64
(
2014
).
https://doi.org/10.1364/ol.39.000064
Crossref
Search ADS
PubMed
 
23.
Y.
Chen
,
S. A.
Ponomarenko
, and
Y.
Cai
,
Sci. Rep.
7
,
039957
(
2017
).
https://doi.org/10.1038/srep39957
Crossref
Search ADS
 
24.
Y.
Cai
,
Y.
Chen
,
J.
Yu
,
X.
Liu
, and
L.
Liu
,
Prog. Opt.
62
,
157
–
223
(
2017
).
https://doi.org/10.1016/bs.po.2016.11.001
Crossref
Search ADS
 
25.
F.
Gori
 and
M.
Santarsiero
,
Opt. Lett.
32
,
3531
(
2007
).
https://doi.org/10.1364/ol.32.003531
Crossref
Search ADS
PubMed
 
26.
G.
Gbur
 and
E.
Wolf
,
J. Opt. Soc. Am. A
19
,
1592
(
2002
).
https://doi.org/10.1364/josaa.19.001592
Crossref
Search ADS
 
27.
L.
Zhao
,
Y.
Xu
, and
S.
Yang
,
Optik
227
,
166115
(
2021
).
https://doi.org/10.1016/j.ijleo.2020.166115
Crossref
Search ADS
 
28.
A.
Dogariu
 and
S.
Amarande
,
Opt. Lett.
28
,
10
(
2003
).
https://doi.org/10.1364/ol.28.000010
Crossref
Search ADS
PubMed
 
29.
Y.
Gu
 and
G.
Gbur
,
Opt. Lett.
38
,
1395
(
2013
).
https://doi.org/10.1364/ol.38.001395
Crossref
Search ADS
PubMed
 
30.
M.
Luo
,
M.
Koivurova
,
M.
Ornigotti
, and
C.
Ding
,
New J. Phys.
24
,
093036
(
2022
).
https://doi.org/10.1088/1367-2630/ac90e0
Crossref
Search ADS
 
31.
L.
Mandel
 and
E.
Wolf
,
Optical Coherence and Quantum Optics
(
Cambridge University Press
,
1995
).
Google Scholar
Crossref
Search ADS
 
32.
G.
Fubini
,
Rend. Acc. Naz. Lincei
16
,
608
(
1907
).
33.
L.
Tonelli
,
Rend. Acc. Naz. Lincei
5
,
246
(
1909
).
34.
M. J.
Bastiaans
,
J. Opt. Soc. Am. A
3
,
1227
(
1986
).
https://doi.org/10.1364/josaa.3.001227
Crossref
Search ADS
 
35.
G.
Gbur
 and
E.
Wolf
,
Opt. Commun.
199
,
295
(
2001
).
https://doi.org/10.1016/s0030-4018(01)01467-5
Crossref
Search ADS
 
36.
F.
Gori
,
Opt. Lett.
30
,
2840
(
2005
).
https://doi.org/10.1364/ol.30.002840
Crossref
Search ADS
PubMed
 
37.
J.
Turunen
,
A.
Halder
,
M.
Koivurova
, and
T.
Setälä
,
J. Opt. Soc. Am. A
39
,
C214
(
2022
).
https://doi.org/10.1364/josaa.475374
Crossref
Search ADS
 
38.
S.
Sahin
 and
O.
Korotkova
,
Opt. Lett.
37
,
2970
(
2012
).
https://doi.org/10.1364/ol.37.002970
Crossref
Search ADS
PubMed
 
39.
Z.
Mei
 and
O.
Korotkova
,
Opt. Lett.
38
,
91
(
2013
).
https://doi.org/10.1364/ol.38.000091
Crossref
Search ADS
PubMed
 
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Supplementary Material

Supplementary materials- zip file