© The Author(s), published by EDP Sciences, 2024

1 Introduction

A multilayer dielectric grating (MLDG) consists of three parts: a periodically corrugated top layer (TG for top grating), a homogeneous connection layer (CL), and a high-reflectivity multilayer dielectric reflector (HR), as shown in Figure 1. In the grating literature, a polarization-independent MLDG (PIMLDG) refers to a MLDG whose diffraction efficiencies, at a given wavelength, in both transverse electric (TE) and transverse magnetic (TM) polarizations are close to 100%. PIMLDGs have found important applications in laser pulse compression [1] and spectral beam combining [2]; therefore, their design is of high practical importance [3, 4].

Figure 1 Schematic diagram of a Multilayer Dielectric Grating (MLDG). The MLDG is divided into three parts: a top grating (TG) with a symmetric trapezoidal profile, groove depth h, top width w, sidewall angle α, and grating period d; a connection layer (CL) with thickness tCL; a high-reflectance reflector (HR). The TG is etched into the top layer of thickness H.

It is generally believed that a high-efficiency MLDG operates on the principle of combining adequate diffraction by the TG and proper thin-film optical interference inside the CL; however, until recently, the detailed physical mechanism was not clearly understood. Although this lack of understanding has not prevented researchers from making good MLDG designs, because Maxwell-theory-based grating codes, which automatically take all physical effects into account, are used as design tools, a clear understanding definitely helps to achieve better designs or to speed up a design process. Recently, the internal mechanism of a MLDG to achieve high diffraction efficiency was expounded in detail [5]. Among the many results of [5], the role of the combined reflection phase is critical to MLDG design. By combined reflection phase, we mean the sum of the reflection phase at the top surface of the HR and round-trip path-length contribution for the propagating beam in the CL. The combined reflection phase, or reflection phase for short, inevitably has been included, in most cases unconsciously, in all previous MLDG designs. The rigorous grating codes numerically take care of everything.

The work in [5] reveals a hidden necessary condition for a MLDG to have polarization-independent high diffraction efficiency. The precise mathematical expression of the condition is given in [5] and Section 3; here suffice it to say that the condition concerns the difference between reflection phases in TM and TE polarizations and is referred to as the phase polarization difference condition (PPDC). To the best of our knowledge, all PIMLDG design works so far do not explicitly address the PPDC, except the work in [6] where the simplified modal method naturally leads to the condition. In a PIMLDG design work, if no means is provided to satisfy the PPDC, no matter how sophisticated the employed optimization tool is, the optimized result may not be the best.

In this work, we present a method to design PIMLDGs that takes advantages of the theoretical results of [5]. As a theoretical paper, reference [5] stresses on the ideal case of 100% diffraction efficiency, whereas in this fabrication tolerance incorporated design paper we consider the practical cases of η* ≤ η ≤ 1, where, and throughout the paper, η = min(η^TE, η^TM), η^TE and η^TM are the – 1st-order diffraction efficiencies in TE and TM polarizations, respectively, and η* is a threshold value. The semi-analytical MLDG theory of [5], especially its diffraction efficiency formula expressed as an elementary function of the reflection phase φ, is the foundation of this paper. To understand the theory of [5] and to implement the algorithm to be presented below, the reader would need to study [5]; however, to follow the present paper, digesting [5] is, hopefully, not absolutely necessary because we have tried to make the description of the diffraction efficiency formula and its dependence on φ self-contained in Section 3. The reflection phase φ as a thread enables us to design first the TG and HR of a PIMLDG separately and then, after the separate designs, combine them to function as a high efficiency PIMLDG. This is one of the computation time saving aspects of the present work. The semi-analytical theory also permits us to satisfy the PPDC explicitly. It turns out that to do so we have to add a few aperiodic thin-film layers atop a periodic stack, as to be explained in Sections 3 and 5.

2 Problem description

As in [5], we make the same five assumptions: (i) The MLDG is lossless and the reflection coefficients of the HR are ${\rho }^{\mathrm{TE}}=\exp \left(\mathrm{i}{\phi }_{\mathrm{HR}}^{\mathrm{TE}}\right)$ and ${\rho }^{\mathrm{TM}}=\exp \left(\mathrm{i}{\phi }_{\mathrm{HR}}^{\mathrm{TM}}\right)$, where ${\phi }_{\mathrm{HR}}^{\mathrm{TE}}$ and ${\phi }_{\mathrm{HR}}^{\mathrm{TM}}$ are reflection phases in TE and TM polarizations, respectively. (ii) The wavelength-to-period ratio λ/d of the TG renders in both cover and CL only two propagating diffraction orders, i.e., 2n_CL/3 < λ/d < 2, where n_CL is the refractive index of the CL and without loss of generality air is taken as the cover. (iii) The TG has the mirror symmetry as in Figure 1. (iv) The incident plane wave is in the −1st-order Littrow mounting. (v) The coupling between the TG and HR via evanescent orders of the TG inside the CL can be neglected. The last assumption in essence is an approximation, which has been shown numerically a very good one for practical design purposes. In addition, we assume that the TG is to be etched into the top layer.

The key parameters defining the grating problem are shown in Figure 1. Ion-beam etched grating profile shapes tend to be trapezoidal with a sidewall angle α. Instead of the upper base width w, it is often more convenient to use the derived, dimensionless quantity duty cycle (or filling factor) f = w/d to specify the trapezoid. In the vertical dimension, the TG-CL combination is over specified with three parameters h, H, and t_CL purposefully; we use two convenient ones as the situation selects. In the figure and the following, we assume the grating ridge and the CL are made of the same material, but this is not a restriction. The top layer thickness H is settled in the multilayer fabrication step and the groove depth h is controlled in the ion-beam etching step. The other parameters, not shown in Figure 1, include thicknesses t_H and t_L and refractive indices n_H and n_L of the high- and low-refractive index materials making up the HR.

Henceforth, the grating parameters that remain constant during design will continue to be referred to as parameters, while those that may vary will be referred to as variables. Thus, all refractive indices, the grating period, and the wavelength (hence incident angle due to the Littrow mounting assumption) are parameters. In this paper, we consider only rectangular and trapezoidal grating profiles, and for the latter the sidewall angle α is treated as a parameter. Among the variables, t_H and t_L obviously belong to HR variables; we lump H together with the rest into TG variables. We further classify the TG variables into two categories: tolerance variables and optimization variables.

In this work, diffraction efficiency plays a role different from that in many design programs: high η per se is not a design target and η* is only used to set up a restriction in the sense of optimization theory. Our design is aimed at finding a combination of TG variables that yields the largest tolerance range for groove depth h and duty cycle f. Consequently, h and f are tolerance variables, subject to restriction η ≥ η*, and only H is an optimization variable. Such a variable assignment is based on our grating fabrication experience and our understanding of the general thin-film coating process: precise control of h and f is more difficult than that of layer thicknesses. We also include ${\delta =\phi }_{\mathrm{HR}}^{\mathrm{TM}}-{\phi }_{\mathrm{HR}}^{\mathrm{TE}}$ as an optimization variable of the TG. This seems to contradict the fact that δ is a property of the HR; however, the inclusion is justified because the PPDC ties the TG and the HR together. This is a threading effect of the combined reflection phase.

The above defined design problem aims at maximizing the range of the tolerance variables that has been greatly reduced by the objective of polarization independence. This potentially necessitates a large number of time-consuming rigorous grating simulations. Therefore, in devising the optimization algorithm the number of grating simulations should be minimized. The semi-analytical theory of reference [5] has made it possible. How this is accomplished is explained in the next two sections.

3 Design method

3.1 Diffraction efficiency formula

In reference [5], under the assumptions stated above, a simple formula for the diffraction efficiency η^σ is derived,$${\eta }^{\sigma }={\sin }^2\left({\zeta }^{\sigma }/2\right),$$(1)where σ = TE, TM is a superscript, not a power, and$${\zeta }^{\sigma }=2\Delta {\theta }^{\sigma }+2\Delta {\chi }^{\sigma }$$(2)with$$\Delta {\theta }^{\sigma }=\arg \frac{ {\tau }_{+}^{\sigma } }{ {\tau }_{-}^{\sigma } }\hspace{1em}\mathrm{and}\hspace{1em}\Delta {\chi }^{\sigma }=\arg \frac{ [ 1-r_{-}^{\sigma }\exp (\mathrm{i} {\phi }^{\sigma }\left)\right] }{ [ 1-r_{+}^{\sigma }\exp (\mathrm{i} {\phi }^{\sigma }\left)\right] },$$(3) $${\phi }^{\sigma }={\phi }_{\mathrm{HR}}^{\sigma }+2\beta t_{\mathrm{CL}}.$$(4)

Equation (4) is the expression of the combined reflection phase that has been referred to several times in Sections 1 and 2. The symbols ${\tau }_{\pm }^{\sigma }$ and $r_{\pm }^{\sigma }$ are related to another set of symbols by ${\tau }_{\pm }^{\sigma }={\tau }^{\sigma }\pm {\tau }^{\text{'}\sigma }$ and $r_{\pm }^{\sigma }=r^{\sigma }\pm r^{\text{'}\sigma }$, where τ^σ and τ^'σ are the 0th- and −1st-order normalized, transmitted, diffraction amplitudes of the TG for the incident plane wave from the cover, and r^σ and r′^σ are the reflected counterparts for the incident plane wave from the CL. For more details, the reader is referred to [5]. In any case, {τ^σ, τ′^σ, r^σ, r′^σ} are the only quantities in this work that require relatively heavy numerical computation because no close-form formulas exist.

The appearance of ${\phi }_{\mathrm{HR}}^{\sigma }$ and 2βt_CL as summands in (4) justifies the use of word “combined” and shows that they are mutually complementary as far as contributing to the diffraction efficiency formula is concerned. A change of ${\phi }_{\mathrm{HR}}^{\sigma }$ can be compensated by an adjustment of t_CL. This is the key that allows the designs of TG and HR to be separated and computation time to be saved. This is also one of the reasons that we view the reflection phase as a thread in the overall MLDG design.

Because τ^σ etc. depend on a, where a stands for a combination of tolerance variables, the function ζ^σ in the above should have been written as ζ^σ(φ^σ; a). For simplicity, we drop and tacitly remember this a dependence, and write ζ^σ as ζ^σ(φ^σ). It follows from equations (2) and (3) that ζ^σ(φ^σ) has these properties: (a) It is a differentiable, 2π-periodic function of φ^σ; (b) it has one and only one pair of maximum and minimum per period; (c) its variation is less than 2π, i.e., max ζ^σ(φ^σ) − minζ^σ(φ^σ) < 2π. It follows further from equation (1) that η^σ being above the threshold η^* means in the plot of η^σ vs. φ^σ there is at least one interval in which the curve η^σ(φ^σ) enters the zone η^* ≤ η^σ ≤ 1 within one period of φ^σ. Figure 6 of [5] schematically illustrates all possible critical-point behaviors of η^σ(φ^σ), excluding the accidental concurrence of sinζ^σ = 0 and (ζ^σ)′ = 0, where the prime denotes derivative. To use the figure for the present purpose, one only needs to draw a horizontal line η^σ = η^* in each of the subfigures. The end points of the nonempty interval(s) can be elementarily determined. In Appendix, it is shown that there can be zero, one or two intervals of η^σ ≥ η^* per period of φ^σ, provided that the end points of the chosen period do not break an interval. We will denote by ψ^σ the set of intervals of η^σ ≥ η^* within a non-breaking period. These properties of ζ^σ(φ^σ) make it possible to predict, in the next subsection, if a TG in a MLDG is capable of delivering an efficiency η ≥ η^* without bringing in a HR.

3.2 TG design strategy

In view of the nature of the design problem, we decide to take an ergodic search approach with respect to the tolerance variable space. Take a rectangular grating as an example and consider the rectangular domain A = [h_min, h_max] × [f_min, f_max]. The choices of h_min, h_max, f_min, and f_max depend on the specific application case. In reference [5] it is shown that high diffraction efficiency of a MLDG more likely occurs when the transmitted −1st- and 0th-orders of the TG have comparable and greater than 25% efficiencies. This information and a preliminary numerical test can help to set the upper and lower limits. A domain A larger than necessary does no harm, except for wasting a little computation time.

In what follows, we use a = {h, f} to denote a pair of tolerance variables and ψ^σ(a) to denote the ψ^σ determined by a. Then, ψ^σ(a) being nonempty means we can find at least one φ^σ per period that makes η^σ(φ^σ; a) ≥ η^*. We collect together all elements a of A that have this property and denote the collection by B_σ. In terms of set-theoretic symbols,$${\mathbf{B}}_{\sigma }\equiv \left\{a\right|a\in \mathbf{A} \mathrm{and} \left(\exists {\phi }^{\sigma }\right)\left|{\eta }^{\sigma }\left({\phi }^{\sigma };a\right)\ge {\eta }^{\mathrm{*}} \right\}=\left\{a\right|a\in \mathbf{A} \mathrm{and} {\psi }^{\sigma }\left(a\right)\ne \mathrm{Ø}\},$$(5)where symbols ∈, |, and ∃ mean phrases “element of”, “such that”, and “there is”, respectively, and Ø denotes the empty set. B_σ is the admissible subset of A for η^σ ≥ η^*. The collection of all elements that are simultaneously in both B_TE and B_TM is then denoted by B and set-theoretically expressed as$$\mathbf{B}\equiv {\mathbf{B}}_{\mathrm{TE}}\cap {\mathbf{B}}_{\mathrm{TM}}=\left\{a\right|a\in \mathbf{A} \mathrm{and} \psi \left(a\right)\ne \mathrm{Ø}\},$$(6)where ∩ means set intersection and$$\psi \left(a\right)\equiv {\psi }^{\mathrm{TE}}\left(a\right)\times {\psi }^{\mathrm{TM}}\left(a\right)=\left\{\right({\phi }^{\mathrm{TE}},{\phi }^{\mathrm{TM}}\left)\right|{\phi }^{\mathrm{TE}}\in {\psi }^{\mathrm{TE}}\left(a\right) \mathrm{and} {\phi }^{\mathrm{TM}}\in {\psi }^{\mathrm{TM}}\left(a\right)\}.$$(7)

In equation (7), for two intervals [x₁, x₂] and [y₁, y₂] along the x and y axes, respectively, [x₁, x₂] × [y₁, y₂] means the rectangular domain in the Cartesian coordinate system Oxy. B is the admissible subset of A for η ≥ η^*. The relationships among A, B_TE, B_TM, and B are schematically shown in the left side of Figure 2a. In general, B may consist of many isolated regions along the h axis, because η^σ is pseudo-periodic in h [7, 8]. In the present analysis, we assume that [h_min, h_max] covers the first period that contains high η^σ or h_max is the upper limit of experimentally achievable etch depth.

Figure 2 Relationships among various sets. (a) The T mapping from point b ∈ B to set γ(b) in (φH, δ) space. (b) The R mapping from point ϕ ∈ D ∩ Φ to B(ϕ) in (h, f) space.

In Appendix, it is shown that for a b ∈ B, ψ^σ(b) can be written as$$\begin{array}{cc}{\psi }^{\sigma }\left(b\right)=I_1^{\sigma }\cup I_2^{\sigma },I_i^{\sigma }=[{\phi }_{2i-1}^{\sigma },{\phi }_{2i}^{\sigma }],& i=1, 2,\end{array}$$(8)where ∪ means set union (z ∈ X ∪ Y, if z ∈ X, or Y, or both). In equation (8), $I_1^{\sigma }$ is by definition nonempty. If $I_2^{\sigma }$ is empty, either $0\le {\phi }_1^{\sigma }<{\phi }_2^{\sigma }\le 2\pi$ or $0<{\phi }_1^{\sigma }<2\pi < {\phi }_2^{\sigma }<{\phi }_1^{\sigma }+2\mathrm{\pi }$; if $I_2^{\sigma }$ is nonempty, either $0\le {\phi }_1^{\sigma }<{\phi }_2^{\sigma }<{\phi }_3^{\sigma }<{\phi }_4^{\sigma }<2\pi$ or $0<{\phi }_1^{\sigma }<{\phi }_2^{\sigma }<{\phi }_3^{\sigma }<2\pi \le {\phi }_4^{\sigma }<{\phi }_1^{\sigma }+2\pi$. The existence of $I_i^{\sigma }$ means there exists a φ^σ such that ${\phi }_{2i-1}^{\sigma }\le {\phi }^{\sigma }\le {\phi }_{2i}^{\sigma }$, excluding the possibility that the two equal signs hold simultaneously (the same exclusion applies to the two inequalities below). Applying this inequality to both TE and TM polarizations and using (4) and the periodicities of ζ^σ and δ leads to$$\begin{array}{cc}\tilde{\phi }+{\phi }_{2i-1}^{\mathrm{TE}}+2\mathrm{m\pi }\le {\phi }_H& \le \tilde{\phi }+{\phi }_{2i}^{\mathrm{TE}}+2\mathrm{m\pi },\end{array}$$(9a) $$\tilde{\phi }+{\phi }_{2j-1}^{\mathrm{TM}}+2\mathrm{n\pi }\le {\phi }_H+\delta \le \tilde{\phi }+{\phi }_{2j}^{\mathrm{TM}}+2\mathrm{n\pi },$$(9b)where m ≥ 0 and n are integers,$$\begin{array}{cc}\tilde{\phi }\equiv 2\beta h-\overline{\phi },& {\phi }_H=2\beta H,\end{array}$$(10)and $\overline{\phi }$, taking the place of ${\phi }_{\mathrm{TE}}^{\mathrm{HR}}$, is a constant that will be determined when we design the HR. Without loss of generality, for the moment we set $\overline{\phi }$ = 0. The restriction m ≥ 0 is due to the physical requirement of H ≥ h. Note that $\tilde{\phi }$ and φ_H have no physical meanings; they are introduced only for notational convenience. The subscript of φ_H serves as a reminder of the linear dependence on H.

For a given b ∈ B, we denote by ${\gamma }_{m,n}^{i,j}\left(b\right)$, γ_m,n(b), and γ(b) the solution of equations (9a) and (9b) for a fixed possible combination of i, j, m, and n, the union of ${\gamma }_{m,n}^{i,j}\left(b\right)$ for all possible combinations of i and j, and the union of γ_m,n(b) for all possible combinations of m and n, respectively. Consider γ_0,0(b). There may be one, two, or four parallelograms in the (φ_H, δ) space. The four boundaries of each parallelogram are defined by the straight lines given by replacing the four ≤ signs in equations (9a) and (9b) by the = signs. The right part of Figure 2a shows the case when both ψ^TE(b) and ψ^TM(b) have two nonempty intervals. If γ_0,0(b) has two parallelograms, they must be mutually shifted parallel to the δ axis or to the line δ + φ_H = constant. The set γ(b) is constructed by replicating γ_0,0(b) horizontally in the positive φ_H direction and vertically in positive and negative δ directions, both in steps of 2π.

The relationship between b and γ(b) can be viewed as a point-to-region mapping T,$$\begin{array}{cc}b\underrightarrow{T}\gamma \left(b\right)& \mathrm{for} b\in \mathbf{B}\end{array},$$(11)which means each element b ∈ B is mapped to an infinite number of parallelogramic regions collectively labeled as γ(b) in the two-dimensional space of (φ_H, δ). Since B is a union of all its elements, we have$$\mathbf{B}\underrightarrow{T}\mathbf{\Phi },$$(12)where$$\mathbf{\Phi }=\underset{b \in \mathbf{B}}{\cup }\gamma \left(b\right).$$(13)

In other words, Φ is the image of B in space (φ_H, δ). Although each b ∈ B gives rise to a series of parallelograms, since B is a two-dimensional continuum, as b runs through B the boundaries of γ(b) are smeared so that Φ no longer resembles parallelograms.

From the construction of Φ we know that for every element ϕ ≡ {φ_H, δ} ∈ Φ there are some b ∈ B such that ϕ ∈ γ(b). In general, for each ϕ ∈ Φ, there are many γ(b)s of different bs. Thus, the set$$B\left(\varphi \right)=\left\{b\right|b\in \mathbf{B} \mathrm{and} \varphi \in \gamma \left(b\right)\}$$(14)is not empty. We may view the relationship between ϕ and B(ϕ) as a point-to-region mapping R,$$\begin{array}{cc}\varphi \underrightarrow{R}B\left(\varphi \right)& \mathrm{for} \varphi \in \mathbf{\Phi },\end{array}$$(15)as shown in Figure 2b. R is the inverse of T in some sense. One of the differences between T and R is that the former is multi-valued (one b is mapped into infinitely many γ_m,n(b)s) while the latter is periodic (periodically distributed ϕs are mapped into the same B(ϕ)).

If we use Σ[G] to denote the area of a set G, then our design task can be stated as to find the ϕ^* such that Σ[Ω(ϕ^*)] ≥ Σ[Ω(ϕ)] ∀ ϕ ∈ Φ, where ∀ means “for all” and Ω(ϕ) is a rectangle inscribed in B(ϕ) whose sides are parallel to the h and f axes. Once ϕ^* is found, we obtain the grating-preferred polarization phase difference δ^* and the top layer thickness H^* under the assumption of ${\phi }_{\mathrm{HR}}^{\mathrm{TE}}=0$.

On the one hand, before searching for ϕ^* we cannot exclude any b ∈ B. On the other hand, from its definition we know that B(ϕ) is 2π-periodic with respect to δ, but only semi 2π-periodic with respect to φ_H because B(2βH + 2π, δ) may be different from B(2βH, δ) if $H<h_{\max }^{\mathbf{B}}$, where $h_{\max }^{\mathbf{B}}$ is the largest h in B. Therefore, without loss of generality, we can choose $\mathbf{D}=\left[2\beta h_{\max }^{\mathbf{B}}, 2\beta h_{\max }^{\mathbf{B}}+2\pi \right)\times [0, 2\pi )$ as the domain of ϕ^*. Any ${\varphi }^{\text{'}}=({\phi }_H^{\text{'}}, \delta \text{'})\in \gamma \left(b\right)$, b ∈ B, obtained according to the above prescription can be folded back into D by using the following formulas:$${\phi }_H=2\beta h_{\max }^{\mathbf{B}}+\mathrm{mod}({\phi }_H^{\text{'}}-2\beta h_{\max }^{\mathbf{B}},2\pi ),$$(16a) $$\delta =\mathrm{mod}\left({\delta }^{\text{'}},2\pi \right),$$(16b)where 0 ≤ mod(u, v) < v is the remainder of u divided by v.

From the above description of the TG design strategy, it is evident that getting B(ϕ) for all ϕ ∈ D ∩ Φ is the most critical step. Note that in general D ∩ Φ is D minus some isolated regions and points, for example, in Figure 2b we have D ∩ Φ = D − E, where E is a region whose points are not covered by any γ(b). Unlike γ(b) that can be semi-analytically and accurately determined, B(ϕ) can only be obtained numerically. Surely, given a ϕ ∈ D ∩ Φ, B(ϕ) could be determined by scanning through all b ∈ B, but that would be a big waste of computation effort. There is a much more efficient way to determine B(ϕ): The traces left in the numerical processes of carrying out the T mappings are sufficient for constructing the inverse R mappings. To illustrate this idea, let us consider a fictitious case of B = {b₁, b₂, b₃}, i.e., B is a set of only three elements (in Sect. 4, both A and Φ will be sampled by equal-spacing grid points). Suppose T(b₁) = D ∩ γ(b₁) = {ϕ₁, ϕ₂, ϕ₃, ϕ₄}, T(b₂) = D ∩ γ(b₂) = {ϕ₃, ϕ₆}, and T(b₃) = D ∩ γ(b₃) = {ϕ₃, ϕ₄, ϕ₅, ϕ₆, ϕ₇, ϕ₈}, all modulo 2π. Then, provided that all mappings T(b_i), i = 1, 2, 3, not just the mapped results, are recorded, it is easy to find that R(ϕ₁) = R(ϕ₂) = {b₁}, R(ϕ₃) = {b₁, b₂, b₃}, R(ϕ₄) = {b₁, b₃}, R(ϕ₅) = R(ϕ₇) = R(ϕ₈) = {b₃}, and R(ϕ₆) = {b₂, b₃}. In fact, it is the numerical discretization that has made the construction of inverse mapping possible. The algorithmic implementation of this idea is presented in Section 4.

3.3 TG and HR unification

After δ^* and H^* have been found we move on, by using a thin-film design program, to design a HR with dual optimization targets: min(|ρ^TE|, |ρ^TM|) ≥ ρ^*, where ρ^* ≈ 1, and mod(δ, 2π) = δ^*. Note that the second target, which is the PPDC that we mentioned in Section 1, may be unreachable, if we confine ourselves to using only a periodic thin-film stack as the HR. It is well-known that under usual application conditions the δ of a quarter-wave periodic stack is π. We first learned from [9] and then verified by ourselves that with common optical thin-film coating materials the δ value of a non-quarter-wave periodic stack does not deviate from ±π by more than ~7% for 10 pairs of periodic layers and ~18% for 20 pairs of periodic layers, provided ρ^* = 0.99. However, a polarization-independent TG design may request a δ^* that is substantially different from ±π. To solve this problem, we add a few aperiodic layers atop the periodic stack using the same high-low refractive index pair. The periodic stack at the bottom assures the high reflectivity and the aperiodic layers at the top produce the needed δ. The minimum total number of needed aperiodic layers depends on the specific application. Our experience shows that 10 is usually enough. Section 5 gives more details.

Incidentally, the present context supports the statement that we made in the Introduction: “if no means is provided to satisfy the PPDC, …, the optimized result may not be the best”. Not knowing the importance of the PPDC, one might not have thought about using a sufficient number of top aperiodic layers. Thus, a multivariable full-scale optimization may produce a design that does not very well satisfy the PPDC; consequently, the fabrication tolerance may not be the best.

After the HR design is successfully completed, if we change $\overline{\phi }=0$ to ${\phi }_{\mathrm{TE}}^{\mathrm{HR}}$ and H^* to H, where ${\phi }_{\mathrm{HR}}^{\mathrm{TE}}$ has the designed value, then δ^*, B(ϕ^*), and Ω(ϕ^*) are all not changed, provided$${\phi }_{\mathrm{HR}}^{\mathrm{TE}}+2\beta H=2\beta H^{\mathrm{*}}+2\mathrm{k\pi }$$(17)for an integer k. From equation (17), we obtain the adjusted top layer thickness$$H=H^{\mathrm{*}}+(2k^{\mathrm{*}}\pi -{\phi }_{\mathrm{HR}}^{\mathrm{TE}})/\left(2\beta \right),$$(18)where k^* gives the minimum $H\ge h_{\max }^{\mathrm{\Omega }},$ $h_{\max }^{\mathrm{\Omega }}$ being the maximum h in Ω(ϕ^*). This completes the design of the PIMLDG.

4 Numerical algorithm

In Section 3, we have assumed that the grating shape is rectangular, so the physically sensible domain A for the two-dimensional variable (h, f) is rectangular. For other grating profile shapes, the shape of A is no longer rectangular. For example, the trapezoidal profile shown in Figure 1 with an α ≠ 0 has an A that is a pentagon (a rectangle with its upper-right corner cutoff by the straight line f + 2htanα/d = f_max). To make the following algorithm more general, we introduce an augmented domain A′, which is always rectangular and includes A as a subset. We then first sample [f_min, f_max] of A′ and [0, 2π] of D by J and N equally spaced points, respectively,$$\begin{array}{cc}{f}_j=f_{\min }+\left(j-1\right)\Delta f, \Delta f=\frac{f_{\max }-f_{\min }}{J-1},& 1\le j\le J,\end{array}$$(19) $$\begin{array}{cc}{\delta }_n=\left(n-1\right)\Delta \delta ,\Delta \delta =2\pi /N,& 1\le n\le N.\end{array}$$(20)

To make the samplings of h and φ_H commensurable, we choose h_min, h_max, and N such that for an integer I,$$\begin{array}{c}\begin{array}{cc}{h}_{\max }-h_{\min }=\left(I-1\right) \Delta h,& \Delta h=\Delta \delta /\left(2\beta \right).\end{array}\end{array}$$(21)

Then, we let$$\begin{array}{cc}{h}_i=h_{\min }+(i-1)\Delta h,& 1\le i\le I\end{array}$$(22) $$\begin{array}{cc}{\phi }_{H,m}=2\beta h_{\max }^{\mathbf{B}}+{\delta }_m,& 1\le m\le N.\end{array}$$(23)

Note that $h_{\max }^{\mathbf{B}}$ is determined and the array φ_H,m is set up only after the for-loop of the main program to be given below commences. Note that A′ is used to set up the sampling points; only points in A will be sampled.

Next, we linearize the double-index pair (i, j) as q such that$$\begin{array}{ccc}q=j+\left(i-1\right)J,& 1\le j\le J,& 1\le i\le I,\end{array}$$(24)and the double-index pair (m, n) as p such that$$\begin{array}{ccc}p=m+\left(n-1\right)N,& 1\le m\le N,& 1\le n\le N.\end{array}$$(25)

The inverses of (24) and (25) can be easily found. Let Q = IJ and P = N². We set up an integer array C of dimension P × (Q + 2), and initialize it so that C(p, 1) = p and C(p, q) = 0 for 1 ≤ p ≤ P and 2 ≤ q ≤ Q + 2 (in practice, the second dimension of C can be much less than Q). After these preparations, we step through all a_q ∈ A′ as follows.

Flag = − 1

Counter = 0

For q = Q: − 1:1

 If a_q ∈ A

  Find the amplitudes {τ^σ, τ′^σ, r^σ, r′^σ} for building function ζ^σ(φ^σ; a_q)

  Find ${\phi }_k^{\sigma }\left(a_q\right)$, σ = TE, TM, k = 1, 2, 3, 4

  If ψ^TE(a_q) = Ø or ψ^TM(a_q) = Ø

   continue

  Else

   Counter = Counter + 1

   If Flag = − 1

    let $h_{\max }^{\mathbf{B}}$ be the current h

    set up the array φ_H,m

    set Flag = + 1

   End

   determine the γ(a_q) in D;

   find all p such that (φ_H, δ)_p ∈ γ(a_q), and

   set t = C(p,2), C(p,2) = t + 1, C(p, t+3) = q

  End

 End

End

After the above loop terminates, Counter gives the total number of points in B, and C(p, l) ≠ 0, l ≥ 3 are the addresses of the sampling points of B(ϕ_p), from which the addresses of Ω(ϕ_p) can be found. Sort the rows of C in descending order of its the second column elements to obtain array C′. Only the leading K rows of C′, for a K ≪ P, will be needed. Find the k^* so that $\mathrm{\Sigma }\left[\mathrm{\Omega }\right({\varphi }_{C^{\text{'}}\left(k^{*},1\right)}\left)\right]\ge \mathrm{\Sigma }\left[\mathrm{\Omega }\right({\varphi }_{C^{\text{'}}\left(k,1\right)}\left)\right]$ for all 1 ≤ k ≤ K, where k and C′(k, 1) are the row number in C′ and the corresponding address p (the row number in C), respectively. Then, C′(k^*, 1) gives the address of ϕ^*, from which we find the H^* and δ^*.

5 Design examples

To validate the design strategy of Section 3 and illustrate the effectiveness of the algorithm of Section 4, we considered two design examples. The diffraction efficiency threshold was set to be η^* = 97%. The main computer program for carrying out the designs was written in MATLAB. The rigorous grating simulation code KAPPA [10] was used for calculating diffraction amplitudes {τ^σ, τ′^σ, r^σ, r′^σ} of the TG, and the commercial thin-film design software TFCalc [11] was used for designing the HR. Note that the harmonic time convention of TFCalc is exp(+iωt), differing from that of KAPPA and our design program. In this section all phase values output by TFCalc are already converted to the exp(−iωt) convention.

The grating structure of the first example (Grating 1) is depicted in Figure 1 with the sidewall angle α set to 0. The TG, CL, and the low-index layers of the HR are all SiO₂, having refractive index n_CL = n_L = 1.46. The high-index layers of the HR are all Ta₂O₅, having refractive index n_H = 2.14. The substrate has refractive index n_S = 1.51. The incident wavelength is 1065 nm and grating period is 900 nm. As a result, β = 2.6695 π/λ.

The ion-beam etcher of our lab limits the maximum etch depth of SiO₂ to be around 2000 nm for the chosen grating period. So, we provisionally set h_max ≈ 2000 nm. Based on our experience that it is difficult to achieve η ≥ η^* with an etch depth less than 500 nm in SiO₂, we set h_min = 500 nm, f_min = 0.15, and f_max = 0.8. The sampling steps were Δδ = π/45 and Δf = 0.01, resulting in N = 90, J = 66, and Δh ≡ Δδ/(2β) = 4.4328 nm. To make h_max − h_min divisible by Δh, we slightly adjusted h_max from 2000 nm to 2002.7 nm, therefore, A = [500, 2002.7] nm × [0.15, 0.8] and I = 340.

The matrix C was numerically calculated in MATLAB. For this example, we found $h_{\max }^{\mathbf{B}}=h_{\max }$. The algorithm of Seibold [12] was used to find the points of Ω(ϕ), then the value of Σ[Ω(ϕ)]. We set K = 500 to search for the global solution ϕ^* which appeared at 9th row of C′. This ϕ^* gave us $H^{*}=h_{\max }^{\mathbf{B}}+3$ Δh = 2016.0 nm, δ^* = 72 Δδ = 1.6 π, and Σ[Ω(ϕ^*)] = 570. The minimum and maximum groove depths in Ω(ϕ^*) were $h_{\min }^{\mathrm{\Omega }}=h_{\min }+112$ Δh = 996.48 nm and $h_{\max }^{\mathrm{\Omega }}=h_{\min }+149$ Δh = 1160.5 nm, respectively. The minimum and maximum duty cycles in Ω(ϕ*) were $f_{\min }^{\mathrm{\Omega }}=f_{\min }+39$ Δf = 0.54 and $f_{\max }^{\mathrm{\Omega }}=f_{\min }+53$ Δf = 0.68, respectively.

Next, the HR was designed aiming at two targets: min(|ρ^TE|, |ρ^TM|) ≥ ρ^* = 0.99 and mod(δ, 2π) = δ^*. We set ρ^* = 0.99 instead of 1 to put more weight on the δ target. For the periodic stack (HL)¹⁴H, we set the optical thicknesses of H and L to λ/4 to provide a broad reflection bandwidth centered at λ. This resulted in thicknesses t_H = 129.5 nm and t_L = 199.5 nm. We introduced 2m aperiodic layers atop the periodic stack, resulting in a layer formula S(HL)¹⁴HL₁H₁…L_mH_mC, where H_i and L_i, 1 ≤ i ≤ m, have subscript dependent thicknesses, and C and S represent CL and the substrate. We gradually increased m from 0 until we found that the design with m = 4 and the following layer thicknesses met both design targets: t(L₁) = 347.9 nm, t(H₁) = 121.2 nm, t(L₂) = 222.6 nm, t(H₂) = 187.8 nm, t(L₃) = 200.0 nm, t(H₃) = 135.3 nm, t(L₄) = 268.7 nm, and t(H₄) = 100.3 nm. This thin-film design gave a ${\phi }_{\mathrm{HR}}^{\mathrm{TE}}=-0.3827 \pi$. We also found that m = 5 could provide any required δ^*. It is worth mentioning that the design of HR is not unique. However, they exhibit the same effect in the design problems discussed in this paper. The design was completed by setting H to be 1294 nm according to equation (18).

The grating structure of the second example (Grating 2) has a symmetric trapezoidal profile, with the sidewall angle α in Figure 1 set to 15°. Both the TG and CL are of high-index material Ta₂O₅, and its refractive index and those of the other layers as well as the wavelength are the same as for Grating 1. The grating period is 714.3 nm (1400 lines/mm) and β = 4.0119 π/λ. Similar to Grating 1, according to the fabrication condition in our lab and our design experience with Ta₂O₅ as the etch layer, we provisionally set h_min = 200 nm, h_max ≈ 1000 nm, f_min = 0.15 and f_max = 0.9. The values of Δδ, Δf, and N were the same as before; therefore, Δh = 2.9496 nm and J = 76. To make h_max − h_min divisible by Δh, we adjusted h_max from 1000 nm to 1002.3 nm, therefore, A′ = [200, 1002.3] nm × [0.15, 0.90], and I = 273. Numerically, we found $h_{\max }^{\mathbf{B}}=h_{\min }+157$ Δh = 663.08 nm which in this case is substantially less than h_max, and the number of sampling points in A is 8299. Based on the points of Ω(ϕ) for the leading K = 500 rows of C′, the solution ϕ^*, which appeared in the second row of C′, gave $H^{*}=h_{\max }^{\mathbf{B}}+55$ Δh = 825.31 nm, δ^* = 56 Δδ = 1.24π, and Σ[Ω(ϕ^*)] = 1008. The minimum and maximum groove depths and duty cycles in Ω(ϕ^*) were $h_{\min }^{\mathrm{\Omega }}=h_{\min }$ + 65 Δh = 391.72 nm, $h_{\max }^{\mathrm{\Omega }}=h_{\min }+112$ Δh = 530.35 nm, $f_{\min }^{\mathrm{\Omega }}=f_{\min }+15$ Δf = 0.3, and $f_{\max }^{\mathrm{\Omega }}=f_{\min }+35$ Δf = 0.50, respectively.

To design HR, we introduced 2m aperiodic layers on top of the λ/4 periodic stack (HL)¹⁴ (t_H = 132.7 nm and t_L = 212.1 nm), resulting in a layer formula S(HL)¹⁴H₁…L_mC. The design with m = 2 met design targets, with thicknesses t(H₁) = 184.3 nm, t(L₁) = 297.3 nm, t(H₂) = 186.9 nm, and t(L₂) = 177.1 nm, and led to ${\phi }_{\mathrm{HR}}^{\mathrm{TE}}=0.5962 \pi$. The design was completed by setting H to be 746.2 nm according to equation (18).

In the above two examples, we have used an even number of aperiodic layers. The choice is somewhat arbitrary, depending on the refractive indices of the top layer of the periodic stack and the CL (when the two are the same or different, the number of added aperiodic layers must be odd or even, respectively). We have found that as long as the total number is large enough, even and odd numbers of aperiodic layers give roughly the same results.

The design results of Gratings 1 and 2 are shown in Figures 3a and 3b, respectively, where the color maps and contour lines represent η calculated by using the formulas in Section 3.1, and the two blue rectangular boxes represent the boundaries of the two Ω(ϕ^*)s. Note that the two Σ[Ω(ϕ^*)] values depend on sampling densities and grating parameters. They are meaningful only within the individual examples; a cross-grating comparison is not useful.

Figure 3 Design results: (a) for Grating 1 and (b) for Grating 2. The color map and contour lines depict η calculated by using the formulas in Section 3.1. The blue rectangular boxes delineate the boundaries of Ω(ϕ*), inscribed within the contour lines of η = 97% (η*). The upper-right white triangular area in (b) belongs to the domain where the lower base width of the trapezoid is greater than fmaxd = 0.9d.

We also computed the η vs. (h, f) maps like those in Figure 3 by using the rigorous code KAPPA. This pair of maps is not shown here because its differences from Figure 3 are almost undetectable. The maximum absolute values of the differences in the blue boxes are less than 0.0013 for Grating 1 and 0.0038 for Grating 2. Easily noticeable differences can only be seen near the lines h = H, where Assumption (v) of [5] is the least accurate.

6 Discussion

6.1 Numerical examples in Section 5

In designing the two example gratings in Section 5, we learned some interesting features of the design algorithm of Section 4. First, some observations about the integer matrices C and C′: 1) Out of P = N² = 8100 rows of C, there are 311 and 1014 rows for Gratings 1 and 2, respectively, whose the second column positions are 0. This suggests that the set Φ covers almost the entire infinite half space (φ_H, δ) with ${\phi }_H\ge 2\beta h_{\max }^{\mathbf{B}}$, leaving only a few small holes in a 2π × 2π unit cell, as schematically shown in Figure 2b. 2) For Gratings 1 and 2, Counter = 3315 and 2177, which means B is much smaller than A (Q = IJ = 22440 and 8299 sampling points in A, respectively). Therefore, choosing a not too large A can significantly save computation time. 3) The maximum value in the second column of C is 1640 and 1553 for Gratings 1 and 2, respectively; therefore, as stated in Section 4, in practice, the second dimension of C can be much less than Q. 4) Our CPU timing tests showed that calculating C and finding Ω(ϕ^*) from C′ took 1.03 s and 0.80 s for Grating 1, and 0.77 s and 0.51 s for Grating 2, respectively, only next to getting the diffraction amplitudes (314 s for Grating 1 and 3506 s for Grating 2) that takes an overwhelmingly large portion of the total computation time. We have seen that in both designs a relatively large K integer was used. If we were content with finding Σ[B(ϕ′)], the maximum-area B(ϕ), the work would be much easier; however, the shape of this B(ϕ′) would very likely be practically useless. For simplicity, we have chosen Ω(ϕ) as the upright inscribed rectangle in B(ϕ). Other useful shapes, such as arbitrarily oriented parallelogram and ellipse, are possible, but they do not take less computation time.

6.2 Computation-time saving

Because our design task of maximizing tolerance range potentially demands a large number of time-consuming rigorous grating efficiency calculations, we have made a serious effort to increase the computational efficiency of our design algorithm. In this effort the semi-analytical theory of [5] has been used to the fullest extent. Besides allowing the use of the reflection phase as a thread to achieve TG-HR separation and unification and direct handling of the PPDC, which we have elaborated, it has contributed to time saving in a number of other aspects. The reader may have noticed that throughout the design process, whenever possible, we avoid direct calculation of diffraction efficiencies, which requires full specification of the MLDG; instead, we work with intervals of η(φ) ≥ η^*. This is a big saving because determining the intervals does not require knowing the HR and each interval naturally contains a range of φ. Even when efficiency calculation is needed in testing if η(φ₀) > η^* (in Appendix) and in post-design verification (as in generating the η maps in Fig. 3), the calculation is done with the analytic formulas in Section 3.1, not a full-matrix code. Moreover, rigorous numerical calculation is needed only for getting the diffraction amplitudes {τ^σ, τ′^σ, r^σ, r′^σ} of the TG, and only once for each pair of (h, f) ∈ A. The rest calculations are all done with the analytic formulas. If a grating code is written cleverly (as in KAPPA), for a fixed TG, the two quadruplets {τ^σ, τ′^σ, r^σ, r′^σ}, σ = TE, TM, can be numerically obtained by launching only one incident plane wave in either the air side or the CL side, instead of separately launching four plane waves (two polarizations and two incident sides. For more details on this point, the reader is referred to reference [13]). Moreover, a code that explicitly takes advantage of Littrow mounting and the grating groove profile symmetry can further cut the computation time to a small fraction of the original cost [14]. And, the last but not the least, the algorithm of Section 4 is designed so that the inverse mapping R is constructed while the forward mapping T is built, and both mappings are saved in integer-address form.

6.3 Extension and limitation

In Sections 3–5, the number of tolerance variables of the TG is two, but the basic design strategy can be extended to cases of more complicated TG structures with more specification parameters, as long as the five assumptions stated in Section 2 are satisfied. For instance, referring to Figure 1, besides f and h, the sidewall angle α may be the third tolerance variable. The TG may consist of two or more etched-through planar layers. Only the dimension of set A is increased from 2 to 3 or more, hence more time is needed for the ergodic search through A.

As stated in Section 4E of [5], the semi-analytical model of MLDG of [5] is inapplicable to 2D-periodic structures or 1D-periodic structures in conical mounting. Consequently, the design method proposed in the present work is inapplicable to these types of structures.

The validity of the algorithm in Section 4 depends on all five assumptions. From a practical point of view, Assumptions (i)–(iii) can be well satisfied. Assumption (v) is a good approximation in most cases except when guided-wave excitation occurs; however, the exceptional case is precisely the intended applications of an MLDG should avoid and it can be evaded by either choosing an appropriate λ/d ratio or selecting a lower high-index material for the HR. Assumption (iv) is restrictive because PIMLDGs are often used with a deviation angle (at Littrow mounting it is 0) over a wavelength band. An extension of the semi-analytical theory of [5] to general off-Littrow mounting is currently underway.

7 Conclusion

We have presented a reflection-phase threaded approach to designing PIMLDGs aiming at maximizing the fabrication tolerance premised by η ≥ η^*. The analytic diffraction efficiency formula of [5] is used at various places of the presented optimization algorithm to save computation time. The threading effect of the reflection phases φ^TE and φ^TM allows the surface-relief grating at the top and the multilayer stack at the bottom first to be designed separately and efficiently, and then to be combined to perform as a PIMLDG. By introducing a few (4 and 8 for the two designed example gratings) aperiodic layers atop of a periodic stack, the PPDC is satisfied explicitly. Two numerical examples are provided to illustrate the design principle and how the optimization algorithm works.

Acknowledgments

The first author would like to thank Kairu Wei, the author of reference [9], for having insightful discussions with her about selected topic of the PhD thesis.

Funding

The work was funded by the National Key Research and Development Program of China (2021YFB2802100).

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

Necessary data are included in this article. Other data will be available from the corresponding authors upon reasonable request.

Author contribution statement

Shen carried out MATLAB programming, numerical calculations, data analysis, and wrote most parts of the manuscript. Li suggested the original idea and offered theoretical support to this work; he also wrote a part of the manuscript. Zeng provided initial guidance of Shen’s work and kept the work on the right track in its development. All authors have read and approved the final manuscript.

References

Appendix

Determination of ψ^σ(b)

For simplicity, we omit the superscript σ; the results to be derived apply equally to both polarizations. It is easy to show that the end points of the interval(s) of η ≥ η^* are solutions of$$\zeta \left(\phi \right)=\pi +\hspace{0.5em}\Delta {\zeta }^{\mathrm{*}},$$(A1a) $$\zeta \left(\phi \right)=\pi -\hspace{0.5em}\Delta {\zeta }^{\mathrm{*}},$$(A1b) $$\zeta \left(\phi \right)=- \pi +\hspace{0.5em}\Delta {\zeta }^{\mathrm{*}},$$(A1c) $$\zeta \left(\phi \right)=- \pi - \Delta {\zeta }^{\mathrm{*}},$$(A1d)where$$\Delta {\zeta }^{\mathrm{*}}=2\arccos \sqrt{{\eta }^{\mathrm{*}}} .$$(A2)

The properties of ζ(φ) stated in the text entails these conclusions: Among equations (A1a)–(A1d), no more than two can be satisfied by one ζ(φ). When two of them are satisfied, the equation pair must be one of the three: (a, b), (b, c), and (c, d). When only one equation is satisfied, it can be any one of the four. When none is satisfied, there are two possibilities: if η(0) > η^* the solution interval is [0, 2π), and otherwise it is empty.

Equation (2) can be rewritten as Δχ(φ) = ζ(φ)/2 − Δθ. Substituting any one of equations (A1a)–(A1d) into this new equation and taking the tangent of both sides, an equation of the form c₁ cosφ + c₂ sinφ + c₃ = 0 can be derived, where c₁, c₂, and c₃ are real constants. This equation has zero, one, and two solutions when $\mathrm{\Gamma }={\left(c_1^2+c_2^2-c_3^2\right)}^{1/2}<0, =0,\hspace{0.5em}\mathrm{and}>0$, respectively (for more details, see [5]). The case of Γ = 0 is theoretically possible but numerically never realized. In practice we exclude this exceptional case. Therefore, equations (A1a)–(A1d) may have two or four distinct solutions. We temporarily label them as ${\phi }_i^{\text{'}}$, and order and restrict them so that ${\phi }_i^{\text{'}}<{\phi }_{i+1}^{\text{'}}$ and $0\le {\phi }_i^{\text{'}}<2\pi$, while recalling that the final solutions are φ_i as in equation (8). Among the ${\phi }_i^{\text{'}}$s, the left end point(s) of the interval(s) can be distinguished from the right one(s) as follows. Let ${\phi }_0= \left({\phi }_1^{\text{'}}+{\phi }_2^{\text{'}}\right)/2$. When there are two solutions, if η(φ₀) > η*, then ${\phi }_i={\phi }_i^{\text{'}}$, i = 1, 2; if not, then ${\phi }_1={\phi }_2^{\text{'}}$ and ${\phi }_2={\phi }_1^{\text{'}}+2\pi$. When there are four solutions, if η(φ₀) > η*, then ${\phi }_i={\phi }_i^{\text{'}}$, 1 ≤ i ≤ 4; if not, then ${\phi }_i={\phi }_{i+1}^{\text{'}}$, 1 ≤ i ≤ 3, and ${\phi }_4={\phi }_1^{\mathrm{\prime }}+2\pi$.

All Figures

	Figure 1 Schematic diagram of a Multilayer Dielectric Grating (MLDG). The MLDG is divided into three parts: a top grating (TG) with a symmetric trapezoidal profile, groove depth h, top width w, sidewall angle α, and grating period d; a connection layer (CL) with thickness tCL; a high-reflectance reflector (HR). The TG is etched into the top layer of thickness H.
In the text

	Figure 2 Relationships among various sets. (a) The T mapping from point b ∈ B to set γ(b) in (φH, δ) space. (b) The R mapping from point ϕ ∈ D ∩ Φ to B(ϕ) in (h, f) space.
In the text

Figure 3 Design results: (a) for Grating 1 and (b) for Grating 2. The color map and contour lines depict η calculated by using the formulas in Section 3.1. The blue rectangular boxes delineate the boundaries of Ω(ϕ*), inscribed within the contour lines of η = 97% (η*). The upper-right white triangular area in (b) belongs to the domain where the lower base width of the trapezoid is greater than fmaxd = 0.9d.

In the text