Improved AWG Fourier optics model

In this paper we present an improved Fourier Optics model to calculate the transmission characteristics between any arbitrary pair of input/output ports (IOPs) of an Arrayed Waveguide Grating (AWG). In this model the input and output sections of the AWG are modeled using the same approximations, thus removing some reciprocity-related inconsistencies present in previously existing models. The expressions which summarize the model are compact and easily interpretable. Simple quasi-analytical expressions are also derived under the Gaussian approximation of the mode field profiles. ©2004 Optical Society of America OCIS codes: (230.7390) Waveguides, planar; (060.4230) Multiplexing; (070.2580) Fourier optics; (070.6110) Spatial Filtering; (350.2460) Filters, interference; (230.1150) All-optical devices; (230.1950) Diffraction gratings References and links 1. M.K. Smit and C. van Dam, “PHASAR-based WDM-devices: Principles, design and applications,” IEEE J. Select. Topics Quantum Electron. 2, 236–250 (1996) 2. H. Takenouchi, H. Tsuda, and T. Kurokawa, “Analysis of optical-signal processing using an arrayedwaveguide grating,” Opt. Express 6, 124-135 (2000) 3. P. Muñoz, D. Pastor, and J. Capmany, “Modeling and design of arrayed waveguide gratings,” J. Lightwave Technol. 20, 661–674 (2002) 4. Y. Chu, X. Zheng, H. Zhang, X. Liu and Y. Guo, “The impact of phase errors on arrayed waveguide gratings,” IEEE J. Select. Topics Quantum Electron. 8, 1122–1129 (2002) 5. P. Muñoz, D. Pastor, J. Capmany and S. Sales, “Analytical and numerical analysis of phase and amplitude errors in the performance of arrayed waveguide gratings,” IEEE J. Sel. Top. Quantum Electron. 8, 1130– 1141 (2002) 6. S. Vallon, P. Chevallier, L. Guiziou, G. Alibert, L.S. How Kee Chun and N. Boos, “40-band integrated static gain-flattening filter,” IEEE Photon. Technol. Lett. 15, 554-556 (2003) 7. J. Zhou, N. Q. Ngo, K. Pita, C.H. Kam, P.V. Ramana and M.K. Iyer, “Determining the minimum number of arrayed waveguides and the optimal orientation angle of slab for the design of arrayed waveguide gratings,” Opt. Commun. 226, 181-189 (2003) 8. A.A. Bernussi, L Grave de Peralta and H. Temkin, “Electric field distribution in a grating of a folded Arrayed-Waveguide Multiplexer,” IEEE Photon. Technol. Lett. 16, 448-490 (2004)


Introduction
Arrayed waveguide gratings (AWGs) are integrated optical devices which make it possible to spatially separate the different wavelengths comprising an optical signal in a very efficient way.They have become one of the key components of wavelength division multiplexing (WDM) optical networks.During the last 15 years a wide variety of applications involving AWGs have been reported, mainly as multiplexers/demultiplexers and wavelength routers in WDM based networks [1], but also in other important applications such as dispersion compensators, antenna beamforming networks, etc. AWGs are constructed as Planar Lightwave Circuits (PLCs), mainly based on Silica-on-Silicon or InP technologies, and can be manufactured and packaged using standard low-cost techniques, thus enabling to cover mass market applications.
AWGs are complex devices, comprising a great number of components and exhibiting a great number of different effects, therefore designers usually rely on theoretical models which partially approximate the device behavior under specific conditions.Fourier Optics modeling of AWGs [2] [3], has become one of the preferred choices because: a) it provides a very simple yet accurate means of simulating the main device characteristics, and b) it can be efficiently simulated by means of the FFT (Fast Fourier Transform).This model has been used, for example, in [4] and [5] to study the impact of phase and amplitude errors on device performance, and in other applications such as [6] to design an integrated static gain-flattening filter, in [7] for determining the minimum number of arrayed waveguides or in [8] to obtain the electric field distribution in the grating of a folded AWG.Nevertheless, as will be shown in this work, such a model does not satisfy the reciprocity principle, that is to say, the forward and reverse transmission coefficients between two arbitraries input and output ports are not the same.
In this paper we present a new AWG Fourier optics model, based on the ones in [2] and [3], which makes it possible to calculate the transmission coefficient between any input/output waveguide pair located at arbitrary positions.The model is developed with a minimum number of approximations and removes some inconsistencies related to the non-reciprocity of the aforementioned AWG Fourier optics model, thus making it more precise and interpretable.
Although the models in [2] and [3] share many things in common, our approach will more closely follow the notation in [3].In Section 2 we develop the new model completely, first for a general situation and then obtaining simpler expressions under the Gaussian approximation of the waveguide mode profiles.Section 3 is devoted to elucidating the obtained results and comparing them with previously existing models.

Layout and notation
The layout of the device under consideration is plotted in Fig. 1.It is comprised of a single monomode input waveguide situated at an arbitrary position d i of the input plane x 0 , followed by a first free propagation region (FPR 1 ) with focal length L f .At the output plane of FPR 1 (x 1 ) the N arrayed waveguides (AW), which are separated by a distance d w from each other, are placed symmetrically with respect to the x 1 axis (N is assumed to be an odd number), to sample the outgoing light.The lengths of the AW differ in a constant amount ∆l which is equal to an integer number (m) of the design wavelength (λ o ) inside the AWs, where m is known as the grating order.The AW drive the light samples from the output plane of FPR 1 to the input plane (x 2 ) of the second FPR (FPR 2 ), which has the same geometrical characteristics as the first one.Light diffracted in FPR 2 illuminates the output plane (x 3 ) where it is collected by a single monomode output waveguide situated at an arbitrary position d o .The input, output and arrayed waveguides are characterized by their power normalized modal field profiles b i (x), b o (x) and b g (x) respectively, which are defined with respect to a local x-axis centered in the waveguides.
The aim of the model proposed in this work is to obtain a closed expression for t(d i ,d o ,ν) which is the transmission coefficient relating the complex amplitudes of the mode of the input waveguide located at position d i (input port IP), and the mode of the output waveguide located at position d o (output port OP), at frequency ν.Notice that, for any combination of the selected parameters d i and d o , the defined transfer function t(d i ,d o ,ν) is the frequency response of the system, that is, its squared amplitude gives the power fraction coupled from the fundamental mode of the input waveguide (b i (x)) to the fundamental mode of the output waveguide (b o (x)).Of course, transmission characteristics of real devices, containing a higher number of IOPs, can be directly analyzed with this model due to the device linearity.

Improved Fourier optics model
To obtain the transmission coefficient t(d i ,d o ,ν), we will start with the optical field in the input plane (x 0 ) and we will calculate the field in the different planes of the device (x 1 , x 2 , x 3 ).As there is a process of unguided light being injected into guided modes of the AW in plane x 1 , it is convenient to distinguish between what exists immediately before the plane (we will refer to this as plane x 1 -and to the corresponding field f 1 () ), and after the plane (we will refer to this as plane x 1 + and to the corresponding field f 1d ()) in the understanding that the field to be calculated in plane x 1 + only includes the transmitted field transported by the guided modes of the AW and thus it represents what would remain in the AW after radiation.A similar situation takes places in plane x 3 .
In this subsection the proposed model for the transmission coefficient between two arbitrarily located IOP's is developed.For the sake of clarity only the main steps of the process are given here, with a more detailed explanation being provided in the appendixes.Also for clarity, we will use upper cases only to refer to Fourier transforms of functions, as is usually done in signal theory.We will thus, for example, define as the spatial Fourier transform of the mode field profile of any waveguide in the array.

Field at x 0
We will consider an input excitation in the input plane f 0 (x 0 ,ν) corresponding to the fundamental mode of the input monomode waveguide at frequencyν.Thus, the input field can be written as where b ν i (x) is the power normalized mode field profile of the input waveguide at frequency ν assumed to be centered at the origin, which depends only on the width and refractive indexes of the input waveguide but not on its position.As we will be dealing only with single frequency excitation and in order to simplify the notation, we will hereafter eliminate the frequency ν and position d i dependence of the fields and will only explicitly use them when necessary.So we will write The field is radiated in the first FPR and, in the Fraunhofer region and under the paraxial approximation [3], the field illuminating the plane x - 1 can be expressed as where F 0 (f) is the Fourier Transform of the input f 0 (x 0 ), and α is the wavelength focal length product or spatial frequency factor, given by where L f is the focal length and n s is the refractive index of the FPR 1 .

Field at x + 1
The field illuminating the AWs is partially coupled to the guided modes of each of the waveguides comprising the array.If coupling between waveguides can be neglected, the guided field coupled to the AW can be written as a sum of each of the individual guided modes where b g (x) is the power normalized mode field profile of any waveguide in the array, assumed to be centered at the origin, and a r is the coupling coefficient to waveguide r, which can be calculated with the following overlap integral Now, instead of assuming that the illuminating field f 1 (x 1 ) is approximately constant throughout the whole AW mode profile b g (x 1 -rd w ) so that a r can be approximated by a sample of the illuminating field f 1 (rdw) multiplied by a constant (the area under b g ()), as is done in [2] [3] [4], we prefer not to make any further approximations yet.The reason is that, although this assumption usually holds true for the amplitude of the illuminating field, this is not the case for its phase.In fact, for non-centered input waveguides, due to the tilting of the wavefronts, the phase of the illuminating field varies linearly with x 1 with a slope proportional to the offset of the input waveguide, and therefore it cannot be considered constant along the AW mode profile.So we prefer to make use of the fact that, as detailed in Appendix A, Eq. ( 7) can be viewed as the sampling of the convolution of f 1 (x)* b g (-x) so that f 1d (x 1 ) can be written as ( ) with δ ω (x 1 ) being the sampling function and g 1 (x 1 ) defined by ( ) where ∏( ) is the rectangular window function.
At this point it is worth noting that under the Gaussian approximation, i.e. with b g (x) being a power normalized Gaussian function with mode field radius ω g and area under the curve 4 2 2 g πω , and assuming that its mode field radius is much smaller that the variation of the illumination f 1 (x 1 ), then b g (x) could have been approximated in Eq. (10) by a weighted delta function, obtaining Then, putting this value into Eq.( 8), we would have which then reduces to Eq. ( 14) of the model proposed in [3].However, we will not do such approximations to avoid reciprocity related issues.

Field at x 2
Neglecting the attenuation in the AWs, the field at the output of the AWs is a spatially phase modulated version of the field at its input.If the lengths of adjacent waveguides differ by a constant value ∆l, which is set to an integer that is a multiple of the operating wavelength in the array, i.e.
where m is the grating order, n c is the refraction index in the waveguides, and ν 0 is the design frequency, then propagation through the AW can be modeled [3] by a simple multiplication of the function The field at the end of the FPR 2 can now be calculated as where F 2 (f) is the Fourier transform of f 2 (x 2 ).Applying the properties of the Fourier transform and after some algebra, detailed in appendix A, we finally get the following for the illumination at the end of FPR 2 where the new function has been introduced.
As traditionally stated [1], we see that Eq. ( 17) describes the field distribution at the output plane x - 3 as a product of a periodic array function (the sum), and the radiation pattern of the singular waveguides of the array (the one proportional to B g (x 3 /α)), and clearly shows the effect of light being focused, with decreasing amplitude, to higher diffraction orders at the output of the AWG.Now it is convenient to point out that the function h(x 3 ), given by Eq. ( 18), depends only on the input section characteristics of the AWG and so it must take into account the two basic phenomena occurring there: windowing and coupling losses.In effect, the right term of the convolution in Eq. ( 18) is a weighted sinc function whose area is always one.This function takes into account the fact that part of the light diffracted in FPR1 is not collected by the finite size (N d w ) of the array aperture, which then acts as a diaphragm, spreading out the light.
When the number of waveguides in the array N tends to infinity, this term tends to δ(x 3 ), and the windowing effect disappears.The left term of the convolution in Eq. ( 18) takes into account the coupling mechanism of FPR 1 to the AW, which due to the reciprocity [1], must have its counterpart in the output section as light being focused to higher diffraction orders.Notice that, as expected, this term again includes the function B g (-x 3 /α), identified in the previous paragraph as the radiation pattern of each of the individual waveguides of the array, and which governed the amplitude of the diffraction orders of the field at the output of the AWG.
Although Eq. ( 17) and ( 18) allow for the calculation the optical field at the output plane x - 3 with a minimum number of approximations, in order to continue progressing, making the results clearer, it is convenient to realize that, for typical AWGs, the left hand term of Eq. ( 18) can be approximated as This approximation holds true as long as one can assume that the radiation pattern of each of the individual waveguides of the array (B g (x 3 /α)) is approximately constant along the width of the input mode profile.This is true in most practical devices in which the non-uniform losses (directly related to the variation of B g (x 3 /α)) are kept to a small value and thus the range of variation of B g (x 3 /α) is limited.Notice that although approximation (19) seems similar to the one we avoided making after Eq. ( 7) this is not actually the case, because the functions in (19) do not have rapidly changing phases that were present in the functions at Eq. ( 7).
Introducing Eq. ( 19) into (18) yields where, following a similar but not identical notation to [3], we have defined The field at the output plane (17) can now be written as Notice that we have recovered the complete notation f 3 (x 3 ,d i ,ν) to clearly highlight that (22) gives the illumination at plane x - 3 due to input excitation located at position d i , at frequency ν.Eq. ( 21) and ( 22) are equivalent to those obtained previously [2] [3] [4], except that these ones predict a different illumination level in the output plane, depending on the position of the input waveguide through the term (B g (d i /α)).This dependence of the illumination at plane x - 3 with the position of the input waveguide d i was lacking in previous models.Notice also that compared to [3] there are two differences in the used notation: a) firstly, the function f M () has been normalized by α 1 so that the weighted Sinc function appearing in Eq. ( 21) always has an area of one.This way, when the windowing effect due to the finite aperture of the array is neglected (making N ∞), f M () tends to the input waveguide mode profile b i (-x 3 ); b) and secondly, as stated before, B g () stands for the Fourier transform of the mode profile of any waveguide of the array, and is not scaled as in that reference.

Field at x + 3
Finally, the light illuminating the output plane of FPR 2 is partially coupled to the guided mode of the output waveguide, located at position d o .As commented previously, we are not interested in the field itself but rather in the coupling coefficient t(d i ,d o ,ν) which relates the complex amplitudes of the modes of the input and output waveguides (centered in d i and d o respectively).This can be calculated with the following overlap integral between the illuminating field and the output power normalized waveguide mode profile, i.e.
Introducing Eq. ( 22) into (23) we get where we have applied the same kind of approximation done in Eq. (19), i.e. we assumed that the radiation pattern of each of the individual waveguides of the array (B g (x 3 /α)) is approximately constant along the width of the mode profile of the output waveguide b 0 (x 3 -d 0 ).Now, using the properties of the convolution operation, and after some manipulations detailed in appendix B, we finally obtain where the new function q M (x) has been defined as the convolution of f M (x) with the mode field profile of the output waveguides b 0 (-x) , i.e.
( ) ( ) ( ) Equations ( 25) and ( 26), which to the author's knowledge have not been previously reported in the literature, directly give the transmission coefficient between any two arbitrarily located IOPs of the AWG.It can be observed that they share many similarities with Eqs. ( 22) and (21), respectively; the later gives the illumination in the plane x 3 -.The main differences between these two sets of equations are due to the different meaning of q M (x) and f M (x 3 ) which are related through convolution as seen in Eq. ( 26): q M (x) is a frequency response (field transfer) related function while f M (x 3 ) is an illumination one.To clarify the meaning of q M (x) we can refer to Fig. 2, which shows an optical imaging system comprised of two lenses with a focal length L f , and a screen with an opening of width Nd w .Two waveguides with power normalized modal field profiles b i (x) and b o (x), are placed at the input and output focal planes of the lenses, respectively, acting as IOPs.The IP is fixed at the centre position, while the OP is offset by an arbitrary value x from the axis of the system.At the design frequency ν o , this system can be seen as a continuous equivalent of the AWG, i.e., it shows the same imaging performance but removing all the sampling effects due to the discrete nature of the AW.Fourier analysis of this optical system predicts that the illumination at the output plane will be an image of the mode field profile of the input waveguide convolved with the sinc-type diffraction pattern of the screen opening.In fact, as depicted in the figure, the illumination in the output plane turns out to be equal to f M (x 3 ) as given by Eq. (21).Coupling to the OP can be calculated as the overlap integral of the output illumination f M (x 3 ) with the offset mode field profile of the output waveguide b o (x 3 -x), but this is just a sample of the convolution of f M (x 3 ) with b o (x 3 ) as done in Eq. ( 26).Thus it is easy to see that the function q M (x) defined in Eq. ( 26) is the transfer function between the centered IP and the x-offset OP of Fig. 2.
It can be noticed that, when the number of arrayed waveguides tends to infinity (N ∞) so that, as commented previously, the weighted Sinc function tends to a delta and the windowing effect disappears, Eq. ( 26) reduces to its left-most term ( ) ( ) and this factor measures the field mismatch existing between the input and offset output mode field profiles.Under such circumstances, if the input and output (power normalized) mode field profiles are equal (b o (x) = b i (x)), so that no input/ouput field mismatch occurs, then where it has been assumed supposed that the mode field profile had even symmetry, as is usually the case.

Gaussian approximation
In many cases it is useful to approximate the mode field profiles of all the waveguides in the AWG (input, output and arrayed waveguides) by power normalized Gaussian functions with mode field radius ω i , ω o , ω g , respectively.In this situation, as demonstrated in Appendix C, expressions (25) (26) can be further simplified.For the most frequent case, in which the input and output mode field profiles are equal, i.e. ω i = ω o , Eqs. ( 25) and ( 26) become Again, if no windowing effect exists (N ∞), the error function tends to unity and Eq.(30) reduces to

Discussion
Equations ( 25) and (26) (or equivalently Eq. ( 29), (30) under the Gaussian approximation) summarize in compact form the proposed Fourier optics model for the transmission coefficient between any pair of arbitrarily located IOPs of an AWG.Compared with previously existing models [2] [3], two are the major contributions of our approach: i) The proposed model gives a much more compact and interpretable expression for the end to end response of the device in which the different phenomena can be easily identified.Effectively previously existing approaches focused on obtaining closed expressions for the illumination at the output plane (for example, Eq. (31) in [3]), from which the transmission coefficient could only be obtained after further calculations (for example, through Eq. (37) in [3]) while, thanks to the definition of the new function q M (), we obtain a much compact expression for this transfer function in which all the characteristic phenomena of AWG behavior can be easily identified; ii) The proposed model removes the reciprocity inconsistencies of previously existing ones, which made them fail when applied to multiple input devices.Thus the results predicted for N o xN o wavelength routers or N o input 1 output multiplexers are necessarily more accurate than those obtained with previously existing approaches.We now discuss these two aspects separately.

3.1.Model interpretation
As stated previously, in the proposed model, summarized in Eq. ( 25) and ( 26), all the characteristic phenomena of AWG behavior can be easily identified: a) The frequency periodicity response and the spatial/frequency discrimination properties of the device appear as a consequence of the q M function in Eq. ( 25).Specifically, maximum signal transmission between any pair of IOPs will be achieved when the following condition is fulfilled where r can be any integer.From this expression the Frequency Free Spectral Range (FFSR) and the Spatial Free Spectral Range (SFSR) which give the space and frequency periodicity can easily be identified as SFSR=α/d w , and FFSR=ν o /m.This results totally agree with those previously reported [3].
b) The function q M (x) takes into account two different phenomena: i) the windowing effect due to the finite aperture of the AW, through the right-most term of Eq. ( 26) and, ii) the field mismatch between input and output waveguides, typically used for passband flattened designs [1], through Eq. ( 27).In the absence of the windowing effect and input/output mode field mismatch, it was shown that q M () is normalized to one, i.e., q M (0)=1.Figure 3 shows plots of the q M () function when input and output waveguides are identical slabs (λ 0 =1.55µm → ν 0 =193.1THz,refractive indexes n 1 =1.454 n 2 =1.451, width W=4µm, Gaussian modal field radius ω i =7.16µm) so that no field mismatch occurs.The focal length of the FPRs is L f =20.85mm (α=2.23•10 3 µm 2 ) and the separation between the arrayed waveguides is d w =26.61µm giving rise to SFSR=α/d w , =837.58µm =117ω i .The grating order is fixed to m=137 yielding FFSR=ν 0 /m=1.41THz.The spatial coordinate is normalized with respect to the modal field radius ω i .Figure 3(a) shows the comparison between the q M (x) function obtained using Eq. ( 26) with the real modes of the slab, and the Gaussian approximation Eq. (30).These results were obtained with N=200 (to replicate the situation in [3]) and show that the Gaussian approximation is only valid in the vicinity of x=0, but it fails to predict the level of the secondary lobes.These results closely resemble those of Fig. 10 in [3], taking into account that: i) there is an offset of approximately 3 dB because q M (x) does not account for the input and output coupling losses due to the discrete nature of the AW, ii) A frequency abscissas axis is used in [3] instead of a distance axis, however they are related through the frequency spatial dispersion parameter γ=FFSR/SFSR= 1.671 GHz/µm used in that reference.
Figure 3(b) shows the detrimental effect of windowing.Two curves are plotted: the first one corresponding to an ideal situation without windowing (q ∞ M ) and the second one corresponding to N=100.In both cases the real mode approximation has been used It can be seen that besides introducing high sidelobes, windowing also introduces losses, so this situation is usually avoided in practical designs.c) The minimum losses (maximum signal transmission) between any pair of IOPs is obtained under condition (32) and is given by where the sum has been approximated only by the dominant term corresponding to r=m.As can be seen in Fig. 3 this approximation holds true in most practical situations.The three right hand side terms in Eq. ( 33) account for the different loss mechanisms in the device: the last term ) 0 ( 1 M q accounts, as discussed previously, for windowing and the input/output field mismatch.In absence of these effects this term becomes one.The other two terms can then be identified as the input ) and output ) ( coupling losses.

Reciprocity and model accuracy
Besides the compactness and clearness of the proposed expressions to calculate the end to end transmission coefficient, the main differences with previous works [2] [3] involve the fact that the input and output section of the AWG has received the same treatment.This way the developed model clearly shows that the input losses, due to coupling from the input waveguide to the AW (through FPR 1 ), and the output losses, due to coupling from the AW to the output waveguide (through FPR 2 ), are, by reciprocity, identical in nature and are related to the existence of higher diffraction orders.This, of course, was a well known fact [1] [3] that was nonetheless unaccounted for in previously existing models.As a result, two benefits are obtained: i) The proposed model is reciprocal, i.e., if the role of input and output are interchanged and we calculate the transmission coefficient from the OP to the IP, the computed transmission coefficient does not vary.This follows immediately from the fact that interchanging the sub-indices 'i' and 'o' in Eq. (25) (26) does not change the computed transmission parameter; b) The developed model is more accurate than previous ones because the input coupling losses are now modeled with a higher order approximation.This approximation must be correct because is the one used to model the output coupling losses in previous approaches, and it is the responsible of the losses non-uniformity of the AWG which is a well documented fact.Thus, the accuracy improvement is a direct consequence of reciprocity, i.e. of using at the input the same modeling of the coupling mechanism that, in the past, has been successfully used at the output.
To clarify this question we notice that, the model in [3] can be recovered as follows.If we approximate B g (di/α)≈B g (0), which is exact for a centered IP, and under the Gaussian approximation (C.2) we get B g (di/α)≈(2πω g 2 ) 1/4 and putting this into Eq.( 22) we obtain: which is identical to that of Eq. (39) in [3], with the exception of a (probably typographical) error which makes the factor 1/d w missing in that work.With this approximation the transmission coefficient between any arbitrary IOPs can be calculated, following the same procedure as before, yielding: with q M () defined by Eq. ( 26).So, Eqs.(34) and (35) summarize, with our notation, the AWG model described in [3] and can be directly compared with our Eqs.( 22) and (25), respectively.It can be clearly seen that in expression (35) reciprocity does not hold because t(d i ,d o ,ν)≠ t(d o ,d i ,ν), and that this is due to the fact that in this model the coupling losses from any IP to the AW have been calculated as if the waveguide would have been centered at the origin and under the Gaussian approximation of the input waveguide mode field profile.It can also be observed that the application of the new model to an AWG demultiplexer (1 IP and N o OPs) yields nearly the same result as previously existing models, because in that case we can set d i =0 in Eq. ( 22) (25) and, under the Gaussian approximation (C.2), we obtain Eq. (34) (35).However, when applying the model presented herein to the study of a N o -input N o -output router or to a N o -input 1-output multiplexer, the results are quite different than those obtained in previous models.To highlight this fact, Fig. 4 shows a comparison of the output plane illumination f 3 (x 3 ,d i ,ν) of a 16x1 AWG multiplexer obtained with previous models (through Eq. ( 34)) and with the proposed model (through Eq. ( 22)).The waveguides, AW and free spatial region's data of the AWG are the same than those in Fig 3(a).which also coincide with those in [3].The red and blue curves correspond to Eq. ( 22) and (34), respectively, when the IP is situated in the center of FPR 1 (i.e., with d i =0) and at the design frequency ν 0 =193.1THz.The black and green dashed curves were also obtained by Eq. ( 22) and (34) respectively, and correspond to an IP situated at d i =239.67µm=33.48wi , and to a frequency shift of 0.4THz with respect to the design frequency ν 0 .These values correspond to the outermost port of the multiplexer with an input waveguide gap of 29.9µm and a channel spacing of 50GHz.
It can be observed that results predicted by the previous model (Eq.( 34)) for both input ports exactly coincide, i.e. no extra losses are accounted for the non-centered IP.However, it can be seen that the new model (Eq.( 22)) predicts for the outermost IP a 0.83 dB penalty with respect to the centered IP.This value exactly matches the non-uniform losses of the same device used as a demultiplexer (see ref. [3]) as should occur by reciprocity.Also, when comparing the results obtained by both models for the centered IP a 0.28 dB difference can be observed.As noted previously, it corresponds to the fact that the Gaussian approximation of the input field profile B g (0)≈(2πω g 2 ) 1/4 was done in Eq. (34), while is not necessary in the new proposed model Eq. ( 22).

Conclusions
An improved AWG Fourier optics model has been developed which provides general expressions to calculate the transmission coefficient between any pair of input/output waveguides located at arbitrary positions.As the approximations made at the input and output sections of the AWG are the same, the obtained equations are more compact and interpretable and, more importantly, they fulfill the reciprocity principle, which was lacking in previously existing models.Gaussian approximation of the input, output and arrayed waveguides mode field profiles, make it possible to further simplify the expressions, leading to very simple final equations which can be extremely useful for the first design process steps.The predicted results essentially coincide with those of the previous models when applied to 1 input N o output demultiplexers, but clearly give better results when applied to N o xN o wavelength routers or N o input 1 output multiplexers.

∑
Introducing Eq. (A.15) and (A.16) in (A.14) yields Now, introducing Eq. (A.17) into (A.13)we get Here, to get a closed expression for f 3 (x 3 ) we still need to calculate G 1 (f).This can be easily obtained by applying the Fourier transform to Eq. (A.10) yielding ( ) where F 1 (f) can also be easily calculated by Fourier transforming f 1 (x 1 ) as where the scaling ( ) and duality properties of Fourier transform have been applied.To make these results more interpretable we prefer working with the new function h(x 3 ) defined as introducing Eq. (A.20) into (A.19) and this into Eq.(A.23) we get where we made use of the scaling property of the convolution Note that, in Eq. (A.24), f 0 () is the input field in plane x 0 (see Eq. ( 3) of the text) so this equation can be written as ( ) With this definition Eq. (A.18) becomes We now start with eq. ( 24) from the main body of the paper which is reproduced here for convenience applying linearity of the integral we can rewrite it as The final integral in this equation can be expressed again as a convolution.Indeed, if we define

Appendix C
In many cases it is useful to approximate the mode field profiles of all the waveguides in the AWG (that is, the input, output and arrayed waveguides) by power normalized Gaussian functions, i.e. and finally, as the convolution of a Gaussian function with a Sinc function can be obtained using the complex error function, we obtain the following analytical expression for for q M (x)

Fig. 1 .
Fig. 1.Layout of the AWG showing one input and one output waveguides located at arbitrary positions

Fig. 2 .
Fig. 2. Optical imaging system showing the meaning of q M (x)

Fig. 3 .
Fig. 3. Coupling between input and output ports of the system of Fig.2 versus normalized position of the output waveguide.

Fig. 4 :
Fig. 4: Illumination at the FPR2 output plane for the centered and outermost input waveguide of a 16x1 multiplexer q M (x) can be finally calculated from Eq. (B.4) and the definition given for f M (x) in Eq. (21) of the text as the Fourier transform of a Gaussian function is a Gaussian function, we can easily calculate .(C.1) into Eq.(25) we obtain the convolution of two Gaussian functions is another Gaussian function, we get

6 )
In the most frequent situation, in which the input and output waveguides are equal, i.e. ω i = ω o , Eq. (C.6) reduces to