Multiple interfacing between classical ray-tracing and wave-optical simulation approaches : a study on applicability and accuracy

In this study the applicability of an interface procedure for combined ray-tracing and finite difference time domain (FDTD) simulations of optical systems which contain two diffractive gratings is discussed. The simulation of suchlike systems requires multiple FDTD↔RT steps. In order to minimize the error due to the loss of the phase information in an FDTD→RT step, we derive an equation for a maximal coherence correlation function (MCCF) which describes the maximum degree of impact of phase effects between these two different diffraction gratings and which depends on: the spatial distance between the gratings, the degree of spatial coherence of the light source and the diffraction angle of the first grating for the wavelength of light used. This MCCF builds an envelope of the oscillations caused by the distance dependent coupling effects between the two diffractive optical elements. Furthermore, by comparing the far field projections of pure FDTD simulations with the results of an RT→FDTD→RT→FDTD→RT interface procedure simulation we show that this function strongly correlates with the error caused by the interface procedure. ©2014 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (050.1950) Diffraction gratings; (080.1753) Computation methods; (230.0230) Optical devices; (350.4600) Optical engineering; References and links 1. M. Lang and T. D. Milster, “Investigation of optics in the 10 – 200 μm size regime,” Opt. Rev. 14(4), 189–193 (2007). 2. T. Bocksrocker, J. B. Preinfalk, J. Asche-Tauscher, A. Pargner, C. Eschenbaum, F. Maier-Flaig, and U. Lemmer, “White organic light emitting diodes with enhanced internal and external outcoupling for ultra-efficient light extraction and Lambertian emission,” Opt. Express 20(S6), A932–A940 (2012). 3. J. Hauss, T. Bocksrocker, B. Riedel, U. Geyer, U. Lemmer, and M. Gerken, “Metallic Bragg-gratings for light management in organic light-emitting devices,” Appl. Phys. Lett. 99(10), 103303 (2011). 4. C.-H. Lin, C.-Y. Chen, D.-M. Yeh, and C.-C. Yang, “Light extraction enhancement of a GaN-based lightemitting diode through grating-patterned photoelectrochemical surface etching with phase mask interferometry,” IEEE Photon. Technol. Lett. 22(9), 640–642 (2010). 5. D. N. Weiss, H.-C. Yuan, B. G. Lee, H. M. Branz, S. T. Meyers, A. Grenville, and D. A. Keszler, “Nanoimprinting for diffractive light trapping in solar cells,” J. Vac. Sci. Technol. B-Microelectron. and Nanom. Struct. 28, C6M98 (2010). 6. R. Dewan and D. Knipp, “Light trapping in thin-film silicon solar cells with integrated diffraction grating,” J. Appl. Phys. 106(7), 074901 (2009). 7. Y. M. Song, J. S. Yu, and Y. T. Lee, “Antireflective submicrometer gratings on thin-film silicon solar cells for light-absorption enhancement,” Opt. Lett. 35(3), 276–278 (2010). 8. S. Feng, X. Zhang, H. Wang, M. Xin, and Z. Lu, “Fiber coupled waveguide grating structures,” Appl. Phys. Lett. 96(13), 133101 (2010). #209272 $15.00 USD Received 31 Mar 2014; revised 13 Jun 2014; accepted 17 Jun 2014; published 23 Jun 2014 (C) 2014 OSA 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.016048 | OPTICS EXPRESS 16048 9. G. Roelkens, D. Van Thourhout, and R. Baets, “High efficiency grating coupler between silicon-on-insulator waveguides and perfectly vertical optical fibers,” Opt. Lett. 32(11), 1495–1497 (2007). 10. T. Pustelny, I. Zielonka, P. Karasinski, and J. Jurusik, “Bragg ’ s grating coupler in planar optical sol-gel waveguides,” Opt. Appl. XXXIV, Vol. 4 (2004). 11. K. Schmitt and C. Hoffmann, High-Refractive-Index Waveguide Platforms for Chemical and Biosensing, Optical Gu, Springer Series on Chemical Sensors and Biosensors (Springer Berlin Heidelberg, 2009), Vol. 7, pp. 21–55. 12. S. Grego, S. Naskar, A. M. Patel, A. Huffman, C. A. Bower, and B. R. Stoner, “Novel optical-waveguide sensing platform based on input grating coupler,” Proc. SPIE 6123, 61230D (2006). 13. M. Prasciolu, D. Cojoc, S. Cabrini, L. Businaro, P. Candeloro, M. Tormen, R. Kumar, C. Liberale, V. Degiorgio, A. Gerardino, G. Gigli, D. Pisignano, E. Di Fabrizio, and R. Cingolani, “Design and fabrication of on-fiber diffractive elements for fiber-waveguide coupling by means of e-beam lithography,” Microelectron. Eng. 68, 169–174 (2003). 14. P. Struk and T. Pustelny, “Design and numerical analyses of the planar grating coupler,” Bull. Polish Acad. Sci. Tech. Sci. 58, 509–512 (2010). 15. D. Taillaert, F. Van Laere, M. Ayre, W. Bogaerts, D. Van Thourhout, P. Bienstman, and R. Baets, “Grating couplers for coupling between optical fibers and nanophotonic waveguides,” Jpn. J. Appl. Phys. 45(8A), 6071– 6077 (2006). 16. G. Roelkens, D. Van Thourhout, and R. Baets, “High efficiency Silicon-on-Insulator grating coupler based on a poly-Silicon overlay,” Opt. Express 14(24), 11622–11630 (2006). 17. S. S. Trieu and X. Jin, “Study of top and bottom photonic gratings on GaN LED with error grating models,” IEEE J. Quantum Electron. 46(10), 1456–1463 (2010). 18. I. J. Buss, G. R. Nash, J. G. Rarity, and M. J. Cryan, “Finite-difference time-domain modeling of periodic and disordered surface gratings in AlInSb light emitting diodes with metallic back-reflectors,” J. Lightwave Technol. 28(8), 1190–1200 (2010). 19. Y. Jin, F. Yang, Q. Li, Z. Zhu, J. Zhu, and S. Fan, “Enhanced light extraction from a GaN-based green lightemitting diode with hemicylindrical linear grating structure,” Opt. Express 20(14), 15818–15825 (2012). 20. S. S. Trieu, X. Jin, B. Zhang, T. Dai, K. Bao, X.-N. Kang, and G.-Y. Zhang, “Light Extraction Improvement of GaN-based Light Emitting Diodes using Patterned Undoped GaN Bottom Reflection Gratings,” Proc. SPIE 7216, 72162Q (2009). 21. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58(5-6), 449–466 (2011). 22. LightTrans,” http://www.lighttrans.com/. 23. A. Rohani, S. K. Chaudhuri, and S. Safavi-Naeini, “Gaussian Beam-Based Hybrid Method for Quasi-Optical Systems,” IEEE Trans. Antenn. Propag. 59(12), 4679–4690 (2011). 24. Breault Research Organization, ASAP Getting started guide,” http://www.breault.com/resources/kbasePDF/broman0108_getstart.pdf. 25. Lumerical Solutions, FDTD Solutions Knowledge Base,” http://docs.lumerical.com/en/fdtd/knowledge_base.html. 26. A. W. Greynolds, “Vector Formulation Of The Ray-Equivalent Method For General Gaussian Beam Propagation,” Proc. of SPIE 0679 (1986). 27. Breault Research Organization, Wave Optics in ASAP,” http://www.breault.com/resources/kbasePDF/brotg0919_wave.pdf. 28. Breault Research Organization, Modeling Incoherent, Extended Sources,” http://www.breault.com/resources/kbasePDF/bropn1153_sources.pdf. 29. Y. Wang, S. Safavi-Naeini, and S. K. Chaudhuri, “A hybrid technique based on combining ray tracing and FDTD methods for site-specific modeling of indoor radio wave propagation,” IEEE Trans. Antenn. Propag. 48(5), 743–754 (2000). 30. Y. Wang, S. K. Chaudhuri, and S. Safavi-Naeini, “An FDTD/ray-tracing analysis method for wave penetration through inhomogeneous walls,” IEEE Trans. Antenn. Propag. 50(11), 1598–1604 (2002). 31. Y. Kawaguchi, K. Nishizono, J.-S. Lee, and H. Katsuda, “Light extraction simulation of surface-textured lightemitting diodes by finite-difference time-domain method and ray-tracing method,” Jpn. J. Appl. Phys. 46(1), 31– 34 (2007). 32. L.-W. Cheng, Y. Sheng, C.-S. Xia, W. Lu, M. Lestrade, and Z.-M. Li, “Simulation for light power distribution of 3D InGaN/GaN MQW LED with textured surface,” Numer Simul Optoelectron Devices 1-2 (2010). 33. Y. Sheng, C. S. Xia, Z. M. Simon Li, and L. W. Cheng, “Simulation of InGaN/GaN light-emitting diodes with patterned sapphire substrate,” 2012 12th Int. Conf. Numer. Simul. Optoelectron. Devices 23–24 (2012). 34. C. Leiner, S. Schweitzer, F. P. Wenzl, P. Hartmann, U. Hohenester, and C. Sommer, “A simulation procedure interfacing ray-tracing and finite-difference time-domain methods for a combined simulation of diffractive and refractive optical elements,” J. Lightwave Technol. 32(6), 1054–1062 (2014). 35. C. Leiner, S. Schweitzer, V. Schmidt, M. Belegratis, F. P. Wenzel, P. Hartmann, U. Hohenester, and C. Sommer, “Multi-scale simulation of an optical device using a novel approach for combining ray-tracing and FDTD,” Proc. SPIE 8781, 87810Z (2013). 36. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995). 37. M. Born and E. Wolf, Principles of Optics (1999), Vol. 1. 38. A. Taflove and S. C. Hagness, “Computational Electrodynamics,” The Finite-Difference Time-Domain Method 2010, 1006 (2005). #209272 $15.00 USD Received 31 Mar 2014; revised 13 Jun 2014; accepted 17 Jun 2014; published 23 Jun 2014 (C) 2014 OSA 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.016048 | OPTICS EXPRESS 16049 39. K. S. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antenn. Propag. 14(3), 302–307 (1966).


Introduction
The last few decades have witnessed remarkable progress of processing speed and storage capacity of computers.This facilitates the use of numerical simulations to address more and more problems in a lot of academic and industrial research and development areas.By this, in comparison with conventional trial and error experimental techniques, the research and innovation speed can be enhanced and costs can be saved.
A remarkable example in this regard is the design and optimization of optical devices that comprise different optical elements.Currently, there is a broad range of optical simulation techniques available to treat suchlike problems.The application ranges of the respective simulation tools mainly depend on the size ranges of the optical elements in the devices in relation to the wavelength of light used.In the "macro" regime, where structure sizes are more than 100 times larger than the wavelengths, the simulation methods are based on the theory of refractive optics (RO), such as classical ray-tracing (RT), Gaussian beam tracing (GBT) or the beam propagation method (BPM) [1].Contrarily, in the so called "nano" regime, for which the ratio between structure size and wavelength becomes smaller than a factor of ten, simulation methods based on ray optics lose their accuracy.Therefore, for an appropriate simulation of devices containing diffractive optical elements (DOEs), numerical methods that are capable of solving problems in the context of wave optics (WO), in particular the Maxwell's equations, are required.Well known examples in this regard are methods like the finite difference time domain method (FDTD), the finite element method (FEM) for more general applications, as well as semi analytic methods like the rigorous coupled-wave analysis method (RCWA), which is restricted to more specific problems [1].However, numerical WO methods require computers with a high memory capacity and often prevent complete simulations of macroscopically optical devices even with the capacity of today's workstation technology.Therefore, an all-embracing simulation of optical devices consisting of both DOEs and refractive optical elements constitute a problem for the abovementioned simulation methods.
On the other hand, as state of the art, DOEs are used in more and more macroscopic devices such as light emitting diodes (LEDs) to enhance, e.g., light extraction or to tailor the emission patterns [2][3][4], in solar cells for light trapping [5][6][7], in fiber-waveguides or planar waveguides as coupling structures [8][9][10], as well as in the field of optical sensing or integrated optics [11,12].As a result of the size restrictions of the different simulation approaches mentioned above and the limited memory capacity of the computers, optical simulations of these devices are limited in many cases to the simulation of individual diffractive or refractive components [13][14][15][16][17][18][19][20].In particular for devices that consist of several diffractive and refractive components these limitations restrict a combined simulation and optimization as well as the consideration of interaction phenomena.Therefore, simulation techniques which are capable of overcoming the abovementioned drawbacks and that allow to handle suchlike complex devices as a whole are highly recommended.
In this regard, several authors already suggested different approaches to interface between RO and WO based simulation tools, using either self-written or commercial RT, GBT, BPM, FDTD or FEM based simulation tools.Such interface approaches allow to address and to consider the individual ROEs and DOEs of a device by an appropriate data transfer from one simulation tool to the other.Still, the biggest challenge in this regard is the set-up of the interface and an appropriate and correct data transfer between WO and RO based simulation tools.
One possibility in this regard is the use of RO-based approaches which allow the consideration of both field and phase data, as shown, e.g., by Wyrowski and Kuhn [21], [22] or Rohani et al. [23].In addition, also suppliers of well-known commercial simulation programs provide tools addressing this problem.For example, Breault Research [24] (RT/GBT/BPM-Program: ASAP) and Lumerical [25] (FDTD Solutions) set up an interface for the exchange of field data between the GBT-mode of ASAP [26,27] and FDTD and vice versa (GBT↔FDTD).This interface allows to handle a wide scope of optical problems, for example, the simulation of an LCD camera or the optical response of microstructures on a DVD surface [24,25].However, this interface is limited to simulation tasks for which in the GBT simulation part the light is focused on the DOE.
Another option to realize an interface between WO and RO based simulation approaches is to combine classical RT with WO approaches, which is a very simple method to address this problem without further restrictions for the RO part of the simulation.Anyway, it has to be considered that in classical RT phase relations between different rays are neglected.This restricts transitions from RT to FDTD to applications for which these relations can be assumed to be random.Nonetheless, this condition is fulfilled for broadband and spatially incoherent light sources like incandescence bulbs, halogen lamps, and light emitting diodes.Another advantage of this approach is that commercial RO simulation programs like ASAP offer a wide range of possibilities to design the emission patterns of such light sources [28].
An example for such a kind of interface was presented already in the year 2000 by Wang et al. who used a combination of classical RT and FDTD for the investigation of radio wave propagation in buildings [29,30].The authors used RT for the simulation of large areas of the building and switched to FDTD for the simulation of areas for which RT was not sufficiently accurate.Further approaches for a combination of RT and FDTD were discussed in [31][32][33].The authors used a combination of FDTD and classical RT for the simulation of the output coupling efficiency of textured surface structures on LED chips.Thereby, FDTD was applied for determining a "scattering" function of these structures mimicking the physical behavior of a DOE in the RT simulation.A further example can be found in the knowledge database of FDTD Solutions [25] where the far fields of different dipole sources of an organic LED were simulated and incoherently combined with FDTD to create an RT source for the ASAP simulation (FDTD→RT).
In addition to these approaches, we recently introduced a step by step simulation procedure for interfacing between the two commercial programs ASAP and FDTD by using the Poynting vector representation of either rays or wave propagation directions to connect the RO and WO approaches [34].Additionally, we demonstrated the accuracy and physical correctness of this interface procedure by comparing the simulation results of an RT→FDTD→RT simulation with real world measurements of a DOE with a goniometer setup [34].In another study [35], we used this RT→FDTD→RT simulation approach for the optimization of an FDTD model setup describing a diffraction grating with zero-order suppression.These studies proved that such an approach allows to get physically accurate results for simulation settings comprising a single interface steps from RT to FDTD and backwards.
In the RT parts of this simulation procedure, DOEs are represented by arbitrary planes.At these planes, the rays are stopped and their positions, directions and fluxes are recorded.Subsequently, FDTD is used to simulate the optical behavior of the DOEs.However, since phase effects are not considered in classical RT, any phase information of the electromagnetic fields is lost during the subsequent FDTD→RT step.For this reason, in our earlier studies the interface procedure was restricted to one FDTD→RT step in order to avoid inaccuracies due to disregarding phase effects between different DOEs.
Nevertheless, in order to be able to apply the simulation procedure for more complex simulation settings -in which e.g. the rays pass the first DOE and hit another DOE or the same DOE again -arbitrary bidirectional stepping between ASAP and FDTD (RT↔FDTD) is indispensable.For such simulation settings the coherence correlation between different DOEs gives reason for additional interference effects which, due to the loss of the phase relations of the different wave fronts upon switching back from FDTD to classical RT, cannot be considered.Therefore, it is necessary to gain a better knowledge on the impact of disregarding these phase relations on the final result.In this regard, the derivation of an error function for estimating the maximum error for the case that these interference effects are not taken into consideration is of crucial importance, not only from the viewpoint of the present interface procedure, but also for all classical RT-WO interface methods mentioned above.
Therefore, in the following we will discuss this on the basis of an optical system which consists of two diffraction gratings and therefore requires an additional FDTD→RT step.In order to minimize the error due to the loss of the phase information in the FDTD→RT step, we will derive an equation for a maximal coherence correlation function (MCCF) which describes the degree of the maximum impact of phase effects between these two different diffraction gratings and which depends on: the spatial distance between the gratings, the degree of spatial coherence of the light source and the diffraction angle of the first grating for the wavelength of light used.This MCCF envelops the oscillations of the square of the time averaged Poynting vector |<S>| 2 caused by the distance dependent coupling effects between the two DOEs.Furthermore, by comparing the far field projections of pure FDTD simulations with the results of an RT→FDTD→RT→FDTD→RT interface procedure simulation we will show that this function strongly correlates with the error caused by the interface procedure.

Simulation -consideration
In order to derive a MCCF that describes the maximal degree of the impact of phase effects between two diffraction gratings one has to adapt the fundamental theory to the case under consideration.This starts with the mathematical formulation of the far-zone form of the Van Cittert-Zernike theorem [36]: ' in which j(r 1 , r 2 ) describes the equal-time degree of coherence at the two points P 1 (r 1 ), P 2 (r 2 ), which are located in the far field of a quasi-monochromatic circular source σ which has a radius a and which is centered at O (see Fig. 1(a)) [36].The vectors s 1 and s 2 are defined as unity vectors, pointing from origin O to the points P 1 (r 1 ) and P 2 (r 2 ) (see Fig. 1(a)).The Van Cittert-Zernike theorem assumes that every point of the source emits spatially completely incoherent light which is an acceptable approximation for the description of incoherent extended light sources like incandescent bulbs, halogen lamps, and light emitting diodes.
Under the assumptions of a uniform intensity distribution (I(r') = constant) and that all points P 1 (r 1 ), P 2 (r 2 ) have the same distance (r = |r 1 | = |r 2 |) from the origin as well as switching to spherical coordinates [r' = (ρcosθ, ρsinθ)] and setting the projections of the unit vectors s 1 and s 2 onto the source plane [s 2⊥ -s 1⊥ ≡ (wcosψ, wsinψ)], it is possible to express Eq. ( 1) written by the integral [36] 2 cos( ) .
Considering the theory of Fraunhofer diffraction for the case of a circular aperture [37], this integral can be solved as in which J 1 (n) is a Bessel function of first kind and order and k represents the wavenumber., where another grating is located.The main challenge for the interface procedure is to derive a coherence correlation function between these two gratings, which allows one to determine the minimal distance d BC between the gratings for which the maximal degree of coherence drops below a predefined threshold value and, thus, allows to disregard possible interference effects without the risk of too large errors in the simulation results.For this reason, Eq. ( 3) has to be converted into a coherence correlation function which depends on the distance d BC in between the two DOEs.
The two points P 1 and P 2 (separated by a distance d 12 ) are located on plane B at the same distance r from the origin O of the source σ.For many reasons it is better to replace the dependency of j(r 1 , r 2 ) on the radius a of the source and the distance d AB by the "divergence angle" α of the incident light on points on plane B, which similarly represents the opening angle of the source σ (see Fig. 1(b)).One of these reasons is that in many cases refractive optical elements, for example lenses, modify the virtual optical distance d AB .Therefore, with the substitution and using the small-angle approximation for α, Equation ( 3) can be written as Assuming a diffraction grating with three orders of diffraction (m = -1, 0, + 1) that are separated by the diffraction angle β, the intensity at every point on plane C is given by a superposition of either the zeroth or the first orders of two points on plane B, which are separated by the distance d 12 = d BC tanβ.For example, in Fig. 1(b) the zeroth order emitted from point P 1 and the first order emitted from point P 2 superpose in point P x .Using Eq. ( 6), one gets an expression for the MCCF  ( , , , tan which depends on the lateral distance between the two gratings d BC , the divergence angle α of the incident light on the first diffraction grating (at plane B) and its diffraction angle β.In Fig. 1(c) the function |j(α, β, d BC )| is plotted for a diffraction angle of β = 43.43°(details later on in the simulation part) for a grating constant of 800 nm and a wavelength of 550 nm.From this image it becomes evident that with decreasing divergence angle α the spatial coherence of the light at plane B increases.Therefore, larger distances d BC are needed to fall below a defined threshold for the MCCF and hence for the maximum error caused by the phase effects.A "divergence angle" of α = 0° would represent a point source.In this case, for which the light is spatially completely coherent, the MCCF is independent from the distance d BC and remains at maximum.

Simulation -simulation setup
In order to verify the assumption that the derived MCCF can be applied for estimating the maximum error caused by the interference effects between the two gratings, different simulations using the program FDTD Solutions from Lumerical [25] were performed and compared with results of simulations using the interface procedure.The FDTD method solves the Maxwell's equations by approximating the derivatives of the Maxwell's curl equations with a central difference approach [38].The simulation area is divided into a fine mesh consisting of cells, so called Yee cells [39], at which the electrical and magnetic field components get discretized and solved in a temporal progressive leapfrog algorithm.The sizes of these unit cells affect the correctness of the simulation results.In case that the size goes to zero, the result becomes exact, on the other hand a coarse mesh is leading to "numerical dispersion" in the simulation area [38] and yields incorrect simulation results.For this reason, the FDTD simulation method is restricted to smaller areas whose sizes depend on the amount of grid points the random access memory of the computer can handle.The simulation area is defined by different boundary conditions, in our case Bloch boundary conditions (B-BC) and perfectly matched layer boundary conditions (PML-BC).While the PML-BC are absorbing any incident field, the B-BC are "re-injecting" incident fields to the other side, enabling the simulation of a periodic grating by simulation of a unit cell of the grating (compare the simulation settings of Fig. 2(b), schematic depicted in Fig. 2(a)).A simple plane wave with a polarization of 45° with respect to the z-axis is used as source for the simulations, averaging polarization effects due to s and p polarization.The source injects a time dependent field signal, which represents a certain range of frequencies in the frequency domain.Additional information about the source, like the mathematical formulation etc., can be found in the literature [38] and the knowledge base of Lumerical [25].Since FDTD is a time domain method, a standard Fourier transformation (SFT) is needed to record the frequency dependent electrical and magnetic field data at a certain position in the simulation area.With further calculations, the averaged pointing vector |<S>| can be determined from the electrical and magnetic field data.Furthermore, a near-to-far field transformation gives the angle dependent intensity distribution in the far field [38].In Fig. 2(a) a sketch of the two dimensional FDTD simulation setting is shown.For simplicity in regard to the near and far field distributions, two identical gratings with identical rectangular shapes and grating constants are used at plane B and plane C (Fig. 1(b)).Figure 2(b) shows an illustration of the FDTD simulation area which covers a unit cell of the grating at plane B as well as a unit cell of the grating at plane C, which are both separated by the distance d BC .Normal to the y-direction, which is the propagation direction of the plane wave, PML-BC are used to absorb all transmitted and reflected fields.Normal to the x-direction B-BC are used to consider the periodicity of the two gratings (see direction of the periodicity in Fig. 2 The derivations discussed in the literature, which lead to Eq. ( 3) and to Eq. ( 7), assume a light source with a circular shape [36].For the two dimensional FDTD simulations of the present case however, a plane wave source of linear shape is used.Therefore, starting again with Eq. ( 1) and using the same assumptions but for a source with a linear shape [r' = (x, 0)] gives the integral ( ) which can be solved similarly to the theory of Fraunhofer diffraction for a rectangular aperture [37] as sin( ) ( , , ) , tan .
With the FDTD simulation setting mentioned above, we tried to validate the distance dependency of d BC in Eq. ( 9) by calculating the MCCF, which is shown by the blue line in Fig. 3(c).These calculations were compared with FDTD simulation results of two (identical) gratings using a partially spatial incoherent light source.Basically, a plane wave source in FDTD is a spatially completely coherent source.However, it is possible to get a spatially incoherent result by superimposing the SFT field data of many simulations using different injection angles from -1.5° to + 1.5° in 0.1° steps to mimic a "divergence angle" α = 1.5° of the source.This value of the divergence angle is close to the divergence angle of the light source of α = 1.4° discussed in our recent study [34].The lattice constants of both gratings were chosen to be 800 nm, which gives reason for an angle of diffraction β = 43.43°for a wavelength of 550 nm (see Fig. 1(c)).), this simulation setup is equivalent to a point source with a radius of zero and therefore represents a spatially completely coherent source.In this case, a distance dependent modulation of the values of |<S>| 2 is caused by destructive and constructive coherence effects between the two gratings.However, the strength of the modulation is independent of the distance d BC , in accordance with Eq. ( 7) and Eq. ( 9) for the case α = 0°, which is why the MCCF remains at maximum for every distance d BC .

Results
Figure 3(b) shows the superposition of many simulation results with different angles of incidence (-1.5° up to + 1.5° in 0.1° steps) representing the case of a partially spatial coherent light source with a divergence angle of α = 1.5°.Contrary to Fig. 3(a), the coherence correlation between the gratings is now decreasing with increasing distance, leading to a damping of this modulation.For large distances, the SFT-Monitor data for the square of the time average Poynting vector |<S>| 2 is approaching a distribution which is independent of d BC representing the result of incoherent coupling of the two gratings.For comparison of this modulation of the incoherent result with the MCCF, the modulation was extracted by subtracting the incoherent part (output for d BC → ∞) of the results as determined by the SFT monitor for every pixel.These absolute values were averaged over the respective pixels.The resulting graph is plotted in Fig. 3(c) (green line), representing the averaged deviation from the incoherent part of the SFT monitor results.As evident, the damping behavior of the modulation fits pretty well the calculated values j(α = 1.5°, β = 43.43°,d BC ) of Eq. ( 9).In this regard, it can be concluded that the calculated MCCF is suitable to estimate the maximal deviation from the incoherent result and furthermore to use the calculated values of j(α, β, d BC ) as a threshold condition for the RT↔FDTD interface.
To highlight this as well as another point of interest which is the impact of the interference effects on the far fields of the gratings, near to far field projections of the data in Fig. 3(b) were performed and compared with an interface simulation of the grating using a RT→FDTD→RT→FDTD→RT interface simulation (Fig. 4(a), (c), (e)).The interface simulation represents the case of a completely incoherent coupling of the two gratings since any phase information is lost due to the RT step in between the respective FDTD calculations of the two gratings.In Fig. 4(a) the thick black line represent the result from the simulations using the interface procedure, while the colored lines represent the calculated far fields obtained from sole FDTD simulations for 50 different distances d BC , ranging from 0 µm to 5 µm in 0.1 µm steps.The FDTD results show a pronounced distance dependency and notable deviations from the results using the interface simulation procedure because of the strong coherence correlation of the two gratings.Figure 4(b) shows a zoom of the zeroth order of the far field.The far field result of the interface simulation procedure is apparently an average of the FDTD far field results, with a deviation of up to ΔI = 18.5% for the maximum of the zeroth order.The use of the variation of the maximum fulfills the requirements of the present study which aims at highlighting the general coherences, or to be more exactly, the area enclosed by the far field intensity distribution should be used as a measure.
In Fig. 4(c), the same RT→FDTD→RT→FDTD→RT far field results are compared with the far field projections of 50 FDTD simulations with the distances d BC ranging from 5 µm to 10 µm in 0.1 µm steps.Similar to Fig. 4(b) the zeroth order as obtained by the interface procedure is in accordance with the average of the FDTD far fields.However, the deviation (ΔI = 8.8%) is much smaller.This is in good accordance with our assumption that the differences between a pure FDTD simulations and the interface simulation are depending on the phase effects between the different orders and furthermore with Eq. ( 9) (the MCCF value becomes smaller with increasing distances d BC ).
In the case of Fig. 4(e), the results of the interface procedure are compared with 50 FDTD far field simulations for large distances d BC, ranging from 20 µm to 25 µm in 0.1 µm steps.In this case the deviations for the zeroth order between FDTD and the interface procedure even drop below 5% (3.7%, see Fig. 4(f)), which could be a reasonable threshold value for a given simulation task.Since this value is comparably small, the impact of disregarded coherence effects on the error of the simulation results will be quite low.Therefore, in this exemplary simulation setting, distances in the range > 25 µm should be suitable for the use of the interface procedure without the risk of a distinct error by disregarding interference effects.Anyhow, the exact value for this distance will also depend on the required accuracy of the given problem, for which the threshold can be set accordingly.
The procedure is not limited to two diffraction gratings, each of them having structures with a rectangular shape.In an additional simulation setting a triangular sawtooth structure was used for the diffraction grating at plane B. In this case the obtained results are quite similar to the ones shown in Fig. 3(a).This can be explained by the fact that in this case only the shape of the structures was changed; however, the angle of diffraction remained the same.On the other hand, changing the grating constant of the diffraction gratings has a much stronger impact on the correlation between the two diffraction gratings; in particular, since the number of orders is increased the respective coherence correlations between these additional orders become important.
Therefore, in another simulation setting (two diffraction gratings with structures having a rectangular shape) we have increased the grating constant from 800 nm to 1400 nm.This gives reason for a far field distribution with 5 different diffraction orders located at −51.78°, −23.13°, 0°, + 23.13° and + 51.78° for an angle of incidence of 0°.In this case, it is necessary to consider all possible combinations of respective coherence correlations between the different orders which can occur in this example (e.g. the −2nd with the 0th or the −2nd with the −1st etc.), to obtain more precise information about the correlation between these two diffraction gratings.As evident from Eq. ( 9) reducing β from 43.43° to 23.12° will lead to a MCCF which is decreasing slower than in the previous example and FDTD simulations up to a distance d BC < 50 µm would be necessary for getting meaningful results.For this reason a divergence angle of α = 3° was used to counteract larger FDTD simulation areas, since this manuscript basically focuses on the verification of the accuracy and applicability of the interface simulation method in general.
For the image shown in Fig. 5(a) the modulation of |<S>| 2 is extracted and compared with the calculated MCCF j(α = 3°, β = 23.13°,d BC ), which represents the coherence correlation between the 0th and the ± 1st diffraction orders.As one can see in this case the modulation of |<S>| 2 is decreasing faster than predicted by the MCCF.This faster decrease of the modulation of |<S>| 2 can be explained by the fact that coherence correlations between the other orders have a significant additional impact.These coherence correlations have to be calculated separately in a way similarly to the calculation of the MCCF and subsequently averaged with different weighting factors to gain an averaged coherence correlation (ACCF) between the two diffraction gratings.We assume that the weighting factors will depend on the intensity of the different orders, because orders with insignificant intensity should not have a large impact on the final distribution even in case that the associated coherence correlation is high.The ACCF plotted in Fig. 5(a) is calculated with equal weighting factors of 1 (arithmetic mean of all coherence correlations) and shows a somewhat better envelopment of the modulation of |<S>| 2 than the MCCF does.However, the values of the ACCF never exceed the values of the MCCF because the diffraction angle between the 1st and the 0th orders is smaller than for every other possible combination of orders, of which the ACCF is composed of.Since it is much more complex and time-consuming to calculate the ACCF than the MCCF (especially in cases with increasing numbers of diffraction orders and therefore of possible combinations between these), the MCCF is still more suitable and ascertainable as a threshold condition.Furthermore, in order to show the applicability and accuracy of the interface procedure for more complex structures, the far field distributions are calculated and depicted in Fig. 5(b),  (c).The corresponding far fields are obtained either by pure FDTD calculations or by using the interface procedure for the two diffraction gratings as a function of different distances d BC and a divergence angle of α = 3°.Similar to the results in Fig. 4(a), smaller distances d BC are leading to higher deviations between FDTD and interface simulation results in Fig. 5(b), whereas in Fig. 5(c) the deviations are much smaller due to higher distances d BC which cause smaller coherence correlations between the two diffraction gratings.These results again demonstrate the possibility of using the interface procedure without the risk of a distinct error by neglecting interference effects for distances of d BC > 25 µm for this combination of diffraction gratings and divergence angle of α = 3° of the light source.However, for such more complex examples und when using a light source with a smaller divergence angle like in the previous example (α = 1.5°), the minimal distance d BC would slightly increase according to the Eq. ( 9).

Summary and conclusions
As shown in our recent work [34,35], the Poynting vector representation of either rays or wave propagation directions can be applied as an interface between two different simulation approaches, ray-tracing and the FDTD method.However, interfacing from FDTD to RT results in a loss of information, namely the phase relations between waves with different propagation directions.Generally, this information loss and the unknown impact of phase effects between different diffractive optical elements prohibit a multiple use of such an interface because of its unknown impact on the accuracy of the simulations and the risk of a massive error summation.
In order to quantify the impact of such interference effects between diffractive optical elements, such as, diffraction gratings, we derived a maximum coherence correlation function j(α, β, d BC ) from the van Cittert-Zernike theorem, which allows to express the coherence correlation between two diffraction gratings in dependence of their distance and the degree of spatial coherence of the light source.For spatially coherent sources, like e.g. point sources or laser sources, the coherence correlation is independent from the distance between the diffractive elements, which restricts the multiple use of such an interface for such light sources.Nevertheless, for more extended light sources, like e.g.blackbody radiators or LEDs, which are only partially spatial coherent the coherence correlation becomes damped for higher distances between the diffractive optical elements.Once the divergence angle of the light source and e.g., the grating constant of the diffraction gratings are known one can calculate the minimum distance between the DOEs for which the impact of interference effects drops below a defined value.In particular, as demonstrated for two diffraction gratings with β = 43.43°separated by a distance d BC > 25 µm and a source with α = 1.5°, results comparable to a sole FDTD simulation of the same setting can be obtained.In another example we verified the extensibility of the basic concept to diffraction gratings having a higher number of diffraction orders, where the correlation is much more complex because of the increasing importance of coherence correlations different from those that are considered by the MCCF.The example shows that the influence of these orders results in an even faster decrease of the oscillations of |<S>| 2 than predicted by the MCCF and which can be described by an ACCF more exactly.However due to the easier calculability of the MCCF especially for more complicated simulation setups, the MCCF is more suitable to serve as threshold condition for the appliance of the interface procedure.
The far field results, plotted in Fig. 4(b),(c) and in Fig. 5(b),(c), indicate that the interface procedure is suitable for a simulation of a wide range of optical devices which contain partially spatial coherent sources where the distance between the diffractive optical elements is large enough so that j(α, β, d BC ) falls below a pre-defined threshold value, which can be set in accordance with the desired accuracy of the simulation task.For simulation tasks for which the distances are smaller it might be possible to combine the two diffractive elements into a single simulation setting in FDTD and use this setting for the interface as a whole.Still, such a combined consideration of the DOEs again will need a compromise between the increasing need of random access memory capacity of the workstation with increasing distance of the DOEs and the decreasing impact of an error, which can be qualified with the help of the derived MCCF.

Fig. 1 .
Fig. 1.(a) Sketch to illustrate the notation for the far-zone form of the van Cittert-Zernike theorem, adopted from [36], (b) Illustration of the parameters needed for the derivation of the MCCF, (c) Plot of the MCCF for a diffraction angle of the first grating β = 43.43°.

Figure 1 (
Figure 1(b) illustrates the case under consideration which consists of two neighboring diffraction gratings.A source σ (Fig. 1(a)) located at plane A (with origin at point O) emits light towards a first diffraction grating located at plane B, which splits the light into different orders towards plane C, where another grating is located.The main challenge for the interface procedure is to derive a coherence correlation function between these two gratings, which allows one to determine the minimal distance d BC between the gratings for which the maximal degree of coherence drops below a predefined threshold value and, thus, allows to disregard possible interference effects without the risk of too large errors in the simulation results.For this reason, Eq. (3) has to be converted into a coherence correlation function which depends on the distance d BC in between the two DOEs.The two points P 1 and P 2 (separated by a distance d 12 ) are located on plane B at the same distance r from the origin O of the source σ.For many reasons it is better to replace the dependency of j(r 1 , r 2 ) on the radius a of the source and the distance d AB by the "divergence angle" α of the incident light on points on plane B, which similarly represents the opening angle of the source σ (see Fig.1(b)).One of these reasons is that in many cases refractive optical elements, for example lenses, modify the virtual optical distance d AB .Therefore, with the substitution

Fig. 2 .
Fig. 2. (a) Scheme of the orientation the FDTD simulation area between the two periodic gratings.(b) Illustration of the FDTD simulation area and its components: The orange lines define the boundary conditions of the simulation area; the green line shows the position of the injected plane wave; the yellow line indicates the position of the standard Fourier transformation (SFT) monitor.The simulation area contains unit cells of grating B and grating C which are separated by the distance d BC.They are replicated in x-direction by the Bloch boundary conditions to simulate two gratings with an orientation as given in Fig. 2(a).
(a)).The yellow line in Fig. 2(b) represents the SFT-Monitor where the electrical and magnetic field data are recorded and subsequently transformed into the |<S>| 2 data, which are plotted in Figs.3(a) and 3(b).

Fig. 3 .
Fig. 3. (a), (b) Simulation results for the simulation model of Fig. 2, recorded with the SFTmonitor for different distances d BC using two identical diffraction gratings with β = 43.43°and (a) for an angle of incident 0°, (b) for a superposition of different angles of incidence from −1.5° up to + 1.5° in 0.1° steps.(c) Modulation of |<S>| 2 extracted from the data of (b) compared with the MCCF j(α = 1.5°, β = 43.43°,d BC ).

Figure 3 (
Figure 3(a) and 3(b) show the square of the time average pointing vector |<S>| 2 for the different pixels of the SFT-Monitor plotted versus different distances d BC .In Fig. 3(a) the simulation results of a plane wave with an injection angle of 0° is shown.According to the applied model (Fig. 1(b)), this simulation setup is equivalent to a point source with a radius of zero and therefore represents a spatially completely coherent source.In this case, a distance dependent modulation of the values of |<S>| 2 is caused by destructive and constructive coherence effects between the two gratings.However, the strength of the modulation is independent of the distance d BC , in accordance with Eq. (7) and Eq.(9) for the case α = 0°, which is why the MCCF remains at maximum for every distance d BC .Figure3(b) shows the superposition of many simulation results with different angles of incidence (-1.5° up to + 1.5° in 0.1° steps) representing the case of a partially spatial coherent light source with a divergence angle of α = 1.5°.Contrary to Fig.3(a), the coherence correlation between the gratings is now decreasing with increasing distance, leading to a damping of this modulation.For large distances, the SFT-Monitor data for the square of the time average Poynting vector |<S>| 2 is approaching a distribution which is independent of d BC representing the result of incoherent coupling of the two gratings.For comparison of this modulation of the incoherent result with the MCCF, the modulation was extracted by subtracting the incoherent part (output for d BC → ∞) of the results as determined by the SFT monitor for every pixel.These absolute values were averaged over the respective pixels.The resulting graph is plotted in Fig.3(c) (green line), representing the averaged deviation from the incoherent part of the SFT monitor results.As evident, the damping behavior of the modulation fits pretty well the calculated values j(α = 1.5°, β = 43.43°,d BC ) of Eq. (9).In this regard, it can be concluded that the calculated MCCF is suitable to estimate the maximal deviation from the incoherent result and furthermore to use the calculated values of j(α, β, d BC ) as a threshold condition for the RT↔FDTD interface.To highlight this as well as another point of interest which is the impact of the interference effects on the far fields of the gratings, near to far field projections of the data in Fig.3(b) were performed and compared with an interface simulation of the grating using a RT→FDTD→RT→FDTD→RT interface simulation (Fig.4(a), (c), (e)).

Fig. 4 .
Fig. 4. (a),(c),(e) Far field transformations of the SFT-monitor data of Fig. 3(b) (colored lines) compared with the simulation results using the interface procedure (thick black line).(a) SFT data for 0 µm < d BC < 5 µm, (c) SFT data for 5 µm < d BC < 10 µm, (e) SFT data for 20 µm < d BC < 25 µm, (b) (d) (f) show zooms of the zeroth order of the far field distributions and the deviations of the FDTD far field results from the interface result.

Fig. 5 .
Fig. 5. (a) Modulation of |<(S)>| 2 compared with the MCCF and an average of the different possible combinations of coherence correlations between the orders (b),(c) Far field transformations of the SFT-monitor data for two rectangular gratings with a grating constant of 1400 nm (colored lines) compared with the simulation results using the interface procedure (thick black line).(b) SFT data for 0 µm < d BC < 15 µm, (c) SFT data for 20 µm < d BC < 25 µm.