Optical response of jammed rectangular nanostructures

1 Introduction

Random packing persists to be an alluring topic, pertinent to fundamental questions in physics, chemistry, and biology [1], [2], [3]. Within the field of optics and photonics, in particular, understanding light–matter interactions in random packed media is crucial and urged by the growing usage of optical sensors and imaging systems in probing complex living cells, liquids, and granular media. In addition, the thriving genre of disordered photonics domesticates randomness toward various applications in light trapping [4], radiative cooling [5], and random lasing [6].

In this work, we investigate the optical response of packed rectangular nanostructures on a surface, as they are commonly employed as dipole scatterers in optical devices [7], [8], [9] for various applications including light harvesting [10] and biosensing [11]. However, a detailed electromagnetic simulation of such an ensemble is a computationally expensive task to perform and one rather seeks the effective medium description as an approximation. Many homogenization theories with varying degrees of applicability and complexity have been developed toward this aim [12], [13]. Bruggeman’s theory models aggregate structure with constituents that are treated on an equal footing and therefore cannot be applied in this case [14]. Maxwell–Garnett approximation, on the other hand, models inclusions dispersed in a continuous host medium. Its analytic simplicity arises from the consideration of only the one-point probability function (density) where convergence is assured under the dilute and long wavelength limit. However, at large packing density, the positional and orientational correlations between the dipoles are not negligible anymore and can drastically alter the effective dielectric constant of the ensemble. The strong contrast expansion method presented in the study by Rechtsmanand and Torquato [13] is rather a generic and exact approach that includes the contribution of high-order point probability functions and thus captures the correction due to the emerged correlations.

In this article, we investigate the influence of the packing dynamic on the optical response of jammed rectangular nanostructures. A comprehensive workflow chart can be found in the supplementary material (section S1). In section 2, we define the packing protocol used in the study and compute the maximum achieved packing densities at various aspect ratios. We proceed in section 3 with point process statistical analysis to unravel the short-range ordering and spontaneous alignment between the packed rectangles. In section 4, we model the ensemble as a two-phase isotropic medium and estimate the effective dielectric constant (εeff) using the strong contrast expansion method. The last section concludes and summarizes the work.

2 Random close packing

We consider the packing of hard rectangles of length l and width w on a square substrate, where the interaction potential ϕ(x) is infinite inside the rectangle region and zero otherwise. Therefore, a stable state is a state with no overlapping elements. The maximum packing density and the state characteristics are protocol dependent. In this work, we implement a collective rearrangement packing protocol which models an abrupt surge or sedimentation of the rectangle concentration near an adsorbing edge which is common in many chemical synthesis techniques [15], [16]. The algorithm starts by placing N identical rectangles of aspect ratio α=l/w at random positions and orientations on a flat substrate of area A_s with an initial packing fraction ϕi=Nlw/As, as shown in Figure 1A. In an iterative procedure, each overlapping rectangle is individually displaced from its initial position with a random radial distance and orientation that are uniformly distributed on the range of [0, 2l] and [0,π], respectively. The procedure persists until the rectangle avoids overlapping or n attempts has been reached after which the process is terminated, and the rectangle is removed from the packing process. To calculate the maximum packing density (ϕmax), we configured the initial density to ϕi=0.6, number of attempts to n = 100,000, and particle number to N = 10,000. Periodic boundaries were used to avoid finite size effects, and the results were averaged over multiple realizations to ensure convergence. The maximum packing density achieved for different aspect ratios is shown in Figure 1B. We note that ϕmax is aspect ratio dependent with a unique cusp around α=1.5. The decrease of ϕmax to the right of the cusp can be explained through the lens of Onsager’s exclusion area principle [17]. The principle states that each rectangle has an excluding area A(l,w,θ) within which other similar rectangles with a relative angel θ are forbidden to occupy if overlapping is to be avoided, as illustrated by the inset in Figure 1B. The ratio of the rectangle area lw to what it excludes on average 〈A〉 across all probable angels [0,π] has been shown to give an accurate estimation of the percolation threshold [18] and is given by the following equation:

where c1=2/π and c2=(2+8/π) in the case of penetrable rectangles. However, our protocol induces a nontrivial positional and orientational correlation that makes the excluded area principle difficult to be derived for ϕmax estimation. In fact, this is the parking problem in 2D which is still an open question to be answered. Nonetheless, we found that the behavior can still be well captured by Eq. (1), using a fitted value of the c₁ and c₂ coefficients, as shown in Figure 1B. However, the model fails to fit the cusp and ϕmax at low aspect ratios. Similar behavior has been observed for packing hard ellipsoids, explained through the isostatic conjecture, and verified through the famous M&M’s experiment [19], [20]. The conjecture states that the mean contact number between packed elements is on average twice the number of degrees of freedom in a jammed configuration. Consequently, deviating from squares to rectangles introduces an additional degree of freedom (orientation) to the packing process, which increases the average contact numbers and consequently ϕmax. Strictly speaking, our packing protocol does not lead to a jammed state rather than a saturated state with an average contact number of zero. Nonetheless, the state can still be treated as a jammed one since the element cannot move once they settled on a nonoverlapping position [20]. Therefore, there are two competing effects occurring in the close packing protocol. On the one hand, adding an extra degree of freedom to squares gives extra flexibility to pack at higher densities (isostatic conjecture), yet on the other hand, increasing the degree of anisotropy prevents short-range ordering and limits the maximum packing density instead (exclusion area principle). These behaviors have a direct influence on the effective dielectric constant of the random jammed media. For example, engineering a random medium with high dielectric constant requires, in general, a high density of packed rectangles. Unintuitively, this can be optimally approached by choosing rectangles of 1.5 aspect ratio, as can be deduced from Figure 1B.

Figure 1:

Random close packing of hard rectangles:

(A) illustration of the collective rearrangement packing protocol with ϕ(x) as the interaction potential and I(i)(ω) is the indicator function of phase i in realization ω and (B) is the computed maximum packing density (ϕmax) as a function of the aspect ratios (α). The blue dots represent the simulated data, and the red curve is a fitting line based on the excluded area principle given by Eq. (1) with c₁ = 0.047 and c₂ = 1.725. The inset illustrates the excluded area of a rectangle of length l and width w with another identical rectangle with a relative tilt θ.

3 Point process statistical analysis

The collective rearrangement packing protocol produces a statistically homogeneous medium that we assume ergodic (any single realization of the ensemble is representative of the ensemble in the infinite area-limit). We start our investigation by performing a stochastic point process analysis. Each rectangle is represented by its two midpoint coordinates (x, y) and the angle (θ) that the longer axis makes with respect to a fixed global reference. We calculate two important statistical descriptors. The first is the radial pair correlation function, which is defined as follows:

where ρ is the number density (number of points per unit area in the infinite area limit) and ρ2(r) is the number of points within a distance of r and r + dr from a reference rectangle. Thus, a deviation of g₂(r) from the unity provides a measure of positional correlation or anticorrelation between the rectangles. The second statistical descriptor is the orientational correlation function gθ(r) defined in the study by Ma and Torquato [21] as follows:

which is an average measure of the degree of alignment between two rectangles within a distance of r and r + dr. Thus, gθ>0 suggests statistically parallel rectangles, gθ<0 suggests statistically perpendicular rectangles, and gθ=0 suggests the lack of any preferential orientation. The two statistical descriptors were applied to the states of maximum packed density that were presented in Figure 1B. The width is normalized (w = 1) without loss of generality, and the descriptors are plotted as shown in Figure 2A and B. We observe the emergence of three distinct features on both surfaces that can be fitted by three lines (p₁, p₂, p₃) that intersect at (α=1,r=1). Transitioning from squares to rectangles causes an interesting trifurcation of the first appearing peak at α=1, as clarified by the top view insets of Figure 2A and B. The trifurcation is a clear sign of the isostatic conjecture. The conjecture states that in a jammed configuration, the mean contact number between packed elements is on average twice the number of degrees of freedom (DOFs) [19], [20]. For simplicity, we can approximate squares as circles and rectangles as ellipsoids by ignoring the sharp edges that have a negligible influence on the packing dynamic. Building on this approximation, there are two DOFs for squares (x, y) and three DOFs for rectangles (x,y,θ). According to the isostatic conjecture, we can conclude that the average contact number increases from 4 to 6 when deforming squares to rectangles. Consequently, the ensemble contains two extra possible configurations, that is, in addition to the original configuration, and forms the trifurcation observed in Figure 2A and B. The low g₂ correlation in the triangular region between p₁ and p₂ is a direct consequence of the excluding area principle. In other words, as the aspect ratio increases, it becomes statistically difficult to pack rectangles within close proximity. Figure 2C shows a cross-sectional plot of both g₂ and gθ for two different aspect ratios. For the square case (α=1), we note a sharp increase in g₂ at r>1 indicating a high probability of occupancy. On contrary to the case of circles packing where the increase is abrupt, rectangle packing has a finite slope that is attributed to their radial asymmetry. In addition, the high probability of occupancy of the first neighbor square induces a negative correlation on the next adjacent regions, a repeated process that explains the oscillatory behavior of g₂(r) that is damped with distance. We also note that gθ is approximately zero, suggesting the lack of any preferential alignment between the packed squares. On the other hand, for α=3, gθ is positively correlated in the region between p₂ and p₃, whereas g_r is negatively correlated. The statistical interpretation indicates that it is highly constrained to place two rectangles in proximity, yet if it is deemed necessary, they must be well aligned. However, such constraint is lifted at p₂ and further relaxed at p₃. The three constraints can be traced geometrically, as illustrated in Figure 2D. The p₂ line equals to the minimum distance when two perpendicular rectangles are not overlapping, that is, when rp2=(α+1)/2. The p₃ line equals to the minimum distance of two stacked rectangles on their longer axes, that is, when rp3=α. We also note that at each relaxation point (p₁, p₂, p₃), a transition in the sign of gθ occurs.

Figure 2:

Point process analysis of maximally random packed rectangles:

(A) radial pair correlation functions g2(r,α) given by Eq. (2) and (B) orientational correlation function gθ(r,α) given by Eq. (3). Three different aspect ratios are highlighted on both surfaces with different colors; the blue line corresponds to α=1, green to α=1.3, and red to α=3. The insets on both surfaces correspond to a top view perspective with three important features highlighted by p₁, p₂, p₃ lines. (C) A cross sectional plot of g₂ and gθ surfaces at α=1 and α=3 and (D) geometrical illustration of the origin of p₁, p₂, p₃ lines, respectively. For each configuration, the dashed circle with radius r_pi represents the shortest distance between rectangles’ center for overlap to be avoided.

We conclude from this analysis the lack of long-range translational or nematic order in the ensemble. The effective permittivity in the 2D plane is thus macroscopically isotropic and polarization independent at all aspect ratios. In addition, high aspect ratios have a destructive behavior on short-range positional order, and therefore, their scattering features will be weaker. Furthermore, the average rectangle orientation after lifting the p₂ constrain is approximately θ≈π/4. In the simple dipole picture, this suggests a spectral redshift for the resonance mode supported along the longer axis (h), as illustrated in the middle configuration of Figure 2D.

4 Strong contrast expansion of the effective dielectric constant

The statistical properties of phase i in two-phase heterogeneous media can be specified by an infinite set of n-point probability functions Sni [22]. In a homogeneous and isotropic medium, the first term reduces simply to the density of phase i(S1(i)=ϕi), whereas the second term S2(i)(r) is interpreted as the probability of finding both endpoints of a line segment of length r in phase i. In the following discussion, we will drop the superscript (i) and implicitly refer to the rectangle phase. The scattering behavior of the ensemble is captured by the spectral density function χv(k) which is the Fourier transform of the autocovariance function χv(r)=S2(r)−S12. Figure 3A shows the spectral density for three different aspect ratios (α=1,1.3,3) at their maximum packing density limit. Structures with low aspect ratios exhibit a clear attenuation in their scattering behavior at small k wave vectors. This is attributed to the suppression of long-range density fluctuations due to the positional ordering of the packed elements. The attenuation becomes weaker as α increases due to the destructive role of the addition of an extra DOF to the packing process, as discussed previously. The scattering behavior is reminiscent of hyperuniformity where χv(k)→0 as k→0 [23]. Hyperuniform structures have been shown to exhibit unique optical properties including the formation of isotropic photonic bandgap that can be used for light guiding and confinement [24], [25]. In addition, the suppressed scattering can lead to transparency in a dielectric medium [26] and enhanced optical absorption in a lossy medium [27]. Although, strictly speaking, our packing dynamic does not lead to hyperuniform structure (χv(k)≠0 as k→0), shared properties are expected. Given that the scattering is weak below the percolation threshold due to the low number of scatterers and similarly weak at maximum packing densities due to the suppression of long-range fluctuations, we expect that there should be an intermediate regime where the scattering events are maximum. We evaluated the zeroth wavevector spectrum density function χv(k=0) for different aspect ratios (α) and located the densities at which the function is maximum, as shown in Figure 3B. Indeed, such an intermediate regime exists between the percolation and maximum packing density limits. In addition, its aspect ratio dependence can also be fitted by the exclusion area principle using Eq. (1).

Figure 3:

Two-phase heterogeneous analysis:

(A) spectral density images for α=1,1.3,3 at ϕmax followed by a radial average plot and (B) density plot as a function of the aspect ratio, showing three curves fitted using the excluded are principle given by Eq. (1): (i) percolation threshold curve (yellow) fitted with c1=2/π and c2=(2+8/π) where ensembles below exhibit weak correlation; (ii) maximum scattering density curve (blue) fitted with c1=0.3780 and c2=1.5255 represents the ensembles at which the zeroth wavevector spectral density function χv(k=0) exhibits a maximum; and (iii) the maximum packing density curve (red) fitted with c1=0.047 and c2=1.725.

From the calculated autocovariance function, we can proceed in calculating the effective-dielectric constant by the strong contrast expansion method. The expressions presented in the study by Rechtsman and Torquato [13] were formulated for 3D random structures. We rederive the method for two-phase medium in 2D and truncate the expansions up to the second order to include the 2-point probability function S2(i)(r). In addition, we focus on the effective dielectric constant experienced by a transverse magnetic plane wave propagating parallel to the 2D medium plane (kz=0,Ez=0). The effective permittivity can be obtained from the following expression (see section S2 for full discussion):

where ϕp is the area density of the packed rectangles, βpq=(εp−εq)/(εp+εq) is the rectangles polarizability (phase p) with respect to the environment (phase q), and βeq=(εeff−εq)/(εeff+εq) is the effective polarizability with respect to the environment. A2p is the coefficient of the second-order correction for a 2D isotropic medium

where H0(1) is the Hankel function of the first kind which acts as propagator function for χV(r). It is important to emphasize that Eq. (5) is approximately valid when λq≫rmax, where λq is the wavelength in the background medium and rmax is the radius at which the positional correlation is negligible (χV(r>rmax)≈0). In other words, all positionally correlated dipoles are assumed to be excited in phase by an external plane wave of wavelength λq. However, the condition is not stringent, and an extension of the applicable wavelength range has been recently shown possible [28]. It can be noted that when A2pβpq2 approaches zero, the expression reduces simply to the well-known Maxwell–Garnett approximation. This is approximately valid for large wavelength (kq→0) and short-range autocovariance function χV(r). In this limit, pure dielectric constituents (εp,εq) result in a pure effective permittivity (Im[εeff]=0). The conditions are approximately met for subwavelength rectangles below the percolation threshold, as illustrated in Figure 3B. This is because, above the percolation threshold, the formation of an infinite sized cluster of overlapping rectangles requires many adjustments of the rectangles’ positions and orientations to reach a nonoverlapping state, resulting in a long-range autocovariance function χV(r). The real and imaginary parts of εeff were evaluated as a function of wavelength for three different aspect ratios (α=1,1.3,3) at their maximum packing densities (ϕmax), as shown in Figure 4. In the studied configurations, χV(r) becomes in the order of 10⁻⁶ for r > 5. Therefore, εeff calculated by Eq. (4) is approximately accurate for λq≫5. It can also be noted that the complex permittivity approaches the one calculated by the Maxwell–Garnett approximation as λq→∞. Furthermore, large aspect ratios damp the resonance and spread the scattering for larger wavelengths. This is consistent with power spectrum density shown in Figure 3A since there exists a proportionality between Im[εeff] and χV(k) as k→0 [13].

Figure 4:

The complex effective dielectric constant (εeff) for jammed rectangular nanostructures calculated by the strong contrast expansion method and plotted as a function of wavelength (λq). The 2D rectangles have width w = 1, height h=α, permittivity εp=9, and are maximally packed (ϕmax) in a background medium with permittivity ϵq=1. (A) The real part of εeff showing an explicit redshift as α increases, and (B) the imaginary part of εeff showing larger scattering as α increases for λq≫rmax. In the long wavelength limit, both the real and imaginary parts of εeff approach the Maxwell–Garnett approximation.

5 Conclusion

Random packed media are a ubiquitous and natural outcome of various chemical synthesis techniques. In the subwavelength limit, the complex inhomogeneous medium can be described by an effective homogeneous one with great accuracy. In this work, we statistically analyze jammed rectangular dipoles under the random close packing protocol for various densities and aspect ratios. The arising microscopic correlations were traced and shown to have direct and indirect consequences on the effective dielectric constant of the medium. Statistical tools and concepts such as Onsager’s excluded area principle, the positional correlation function, and the orientational correlation function, are of great utilities in describing the state of the ensemble and deduce some of the optical characteristics such as polarization dependence and spectral shifts. To study the influence of structural correlations on the macroscopic optical response, we accommodate the strong contrast expansion method to two-dimensional structure and use it to estimate the effective dielectric constant for the generated ensembles. This allows us to capture various effects beyond what Maxwell–Garnett approximation can, such as scattering enhancement and suppression as well as correlation-induced spectral shift. This work paves a systematic path toward engineering random medium with tailored optical properties.