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.

## I. INTRODUCTION

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 effect^{3,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 focusing^{15–17} and splitting^{18,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 representation^{25} 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.

## II. PROPAGATION OF PARTIALLY COHERENT BEAMS

*z*= 0 is defined as

^{31}

*E*(

**) is a member of a statistical ensemble of monochromatic realizations, asterisk denotes complex conjugation, $\rho j\u2032=(xj\u2032,yj\u2032)$, 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,

**= (**

*κ**k*

_{x},

*k*

_{y}) contains the transverse wave vector components. Therefore, the CSD can be transformed to the angular correlation function (ACF) with the generalized Wiener–Khintchine theorem,

*z*by imposing a propagation phase of the form

*k*

_{j}= |

**k**

_{j}|,

**k**

_{j}being the wave vector in free space at a frequency

*ω*. In general,

*k*

_{zj}can be complex, so evanescent waves are implicitly accounted for. The spatial domain correlation function can be retrieved with

*ρ*_{j}= (

*x*

_{j},

*y*

_{j}) contains the transverse coordinates at the plane

*z*> 0. Note that the initial plane coordinates were defined under Eq. (1) with the primed symbols.

^{32,33}we can rearrange the order of integrals if the function to be integrated fulfills the requirement,

### A. Fresnel propagation

*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,

**= (**

*ρ*

*ρ*_{1}+

*ρ*_{2})/2 and Δ

**=**

*ρ*

*ρ*_{2}−

*ρ*_{1}, difference coordinates, we can rewrite the Fresnel propagation equation in the following, inspiring, form:

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

*ρ***′ factors in the exponent and introducing the Wigner representation of the CSD as**

*ρ*^{34}

**has two contributions of equal weight: the on-diagonal element corresponds to stretching, and the off-diagonal element causes shearing.**

*T*

*ρ*_{1}=

*ρ*_{2}=

**, which is equivalent to setting Δ**

*ρ***= 0 in the average and difference coordinates,**

*ρ**W*(

**, 0,**

*ρ**z*) =

*S*(

**,**

*ρ**z*). This yields

### B. Fraunhofer propagation

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\beta /2$.^{35–37} Here, *w*_{0} 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.

**′Δ**

*ρ***′ in Eq. (11) can be neglected, as for propagation distances**

*ρ**z*≫

*z*

_{R}, 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.,

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

*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,**

*κ*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.

## III. GENUINE REPRESENTATION

^{25}

*p*($v$) is a non-negative weight function, and

*H*(

**′, $v$) represents an arbitrary kernel, usually written as a Fourier kernel**

*ρ***(**

*f***′) is a function that depends only on the initial plane position**

*ρ***′, and**

*ρ**τ*(

**′) is in general a complex amplitude function.**

*ρ**g*(Δ

**).**

*f*However, both the function *g*(Δ** f**) itself and its argument

**cannot be chosen arbitrarily but must fulfill at least the following three conditions:**

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

*g*(Δ) must be Fourier-transformable, i.e., $g:L1(R2)\u2192C$ and;*f*if $f(\rho 1\u2032)=f(\rho 2\u2032)$, 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

**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,**

*f**p*(

**v**) and

**can be any linear combination of elementary functions that fulfills the conditions outlined earlier.**

*f*Similarly, the complex amplitude *τ*(** ρ**) cannot be chosen arbitrarily, but must abide by the condition of being a finite-energy-carrying function, i.e., $\tau (\rho )\u2208L2(R2)$, and it must be continuous and smooth so that

*τ*∈

*C*

^{∞}to ensure compliance with Maxwell’s equations.

*g*(Δ

**). Although it was noted earlier that the correlations in the field are shaped by the function**

*f**g*(Δ

**) itself, this does not necessarily mean that**

*f**g*(Δ

**) is the spectral DOC. The latter is, in fact, defined as**

*f**g*(Δ

**), i.e., $W(\rho 1\u2032,\rho \u20322)=C(\rho 1\u2032,\rho 2\u2032)g(\Delta f)$, above, we get**

*f*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.

*g*(Δ

**′). To gain more insight on the physics of the above-mentioned equation, let us then consider the quasihomogeneous case and consequently assume that $C(\rho 1\u2032,\rho 2\u2032)=S(\rho 1\u2032)S(\rho 2\u2032)\u2248S(\rho \u2032)$ and**

*f**g*(Δ

**′) =**

*f**μ*(Δ

**′) depends only on the distance between coordinates (i.e., it follows the Schell-model), $\Delta f\u2032=\Delta \rho \u2032=\rho 2\u2032\u2212\rho 1\u2032$, while being narrow with respect to**

*f**S*(

**′). With these approximations, we can perform the integration over Δ**

*ρ***′ in the above-mentioned expression, which gives the following result:**

*ρ**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.

## IV. ANALYTICAL INVERSE DESIGN

We may start our design process by specifying that we have an initial plane CSD $W(\rho 1\u2032,\rho 2\u2032)$, 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(\rho 1\u2032,\rho 2\u2032)$. 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.

*W*(

**′, Δ**

*ρ***′) ≈**

*ρ**S*(

**′)**

*ρ**μ*(Δ

**′), in which case we have**

*ρ**μ*(Δ

**′)| ≤ 1, needs to be met, and |**

*ρ**μ*(0)| = 1 holds, which is always fulfilled when energy is conserved.

### A. Shifting beams

*ρ*_{0}is the shift parameter and $zR=kw02\beta /2$ is the Rayleigh range introduced earlier, with $\beta =1+w02/\sigma 02\u22121/2$ chosen to correspond to a Gaussian Schell-model beam, with

*w*

_{0}and

*σ*

_{0}being the rms width and coherence width of the source, respectively. In the case of a quasihomogeneous field with

*σ*

_{0}≪

*w*

_{0}, the Rayleigh range is approximated as

*z*

_{R}=

*kw*

_{0}

*σ*

_{0}/2.

### B. Intensity redistribution

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.

^{38},

### C. Self-splitting beams

^{39},

## V. NUMERICAL INVERSE DESIGN

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.

### A. Implementation

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

*κ*_{c}=

*k*(

**−**

*ρ***′)/**

*ρ**z*and, therefore,

**′ =**

*ρ***−**

*ρ**z*

*κ*_{c}/

*k*,

*κ*_{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).

*c*

_{n}is a non-negative weight function and

*e*

_{n}(

*x*) are shifted and tilted Gaussian coherent modes, with an explicit form of

*x*

_{n}and

*a*

_{n}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.

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.

### B. Transient effects

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.

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.

## VI. DISCUSSION AND CONCLUSIONS

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 *x*_{n}. 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 *a*_{n} can be realized by imposing a linear phase at each pinhole. This can be performed via a 4*f* 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Inverse Scattering Problems in Optics*

*Optical Coherence and Quantum Optics*