1. Introduction
Most applications of holographic photopolymers require high diffraction efficiency but also low scatter. Examples include data storage, where scatter contributes to the bit error rate, and augmented reality where scatter contributes to haze that detracts from visibility. Unfortunately, it is commonly observed in material development that these two quantities are coupled. That is, high dynamic range materials capable of recording high diffraction efficiency also typically develop large scatter. It is thus a significant challenge to develop materials formulations and processing methods resulting in a high signal to noise ratio (SNR).
Previous theoretical studies of scatter development in holographic photopolymers have identified noise gratings seeded by initial scatter centers as the source of recorded haze [
1,
2]. The central concept of these studies is that recorded scatter must be considered as an ensemble of weak holograms, a small fraction of which are Bragg-selectively replayed by the reference and object beams. This amplifies a fraction of the scatter which then interferes with the writing beams to write further noise gratings. If grating development is fast in comparison to the recording rate, this process dynamically forms characteristic arcs and lines of scatter at distinct angles. In contrast, an approach where particles are thought of as seeds for self-written spatial filaments [
3] would be expected to scatter approximately isotopically. Reported experimental scatter patterns are generally consistent with the former theory, showing strong angular dependence in narrow arcs and lines consistent with Bragg-degenerate noise grating growth. Monte Carlo numerical studies are consistent with these results [
4].
The observed scatter patterns can be explained by the Ewald sphere analysis shown in
Figure 1. In the left two images, object and reference beams, shown as blue arrows, write a holographic grating with positive and negative side bands, shown as the two heavy black arrows. Scatter in a cone around the reference beam (top) or object beam (bottom), shown as a circular patch on the Ewald sphere, interferes with the object beam to write weak noise gratings. The positive and negative side bands of these noise gratings, shown as thin black lines, reside on two-dimensional spherical surfaces centered on the holographic grating K vector. The interference of the object beam with the reference (top) and object (bottom) scatter thus generates four two-dimensional noise grating manifolds, shown as shaded spherical patches. The right-hand plots show that a curved one-dimensional slice of two of these manifolds is Bragg matched when replayed by the reference beam. The slice, shown as heavy red lines, is found at the intersection of the Ewald sphere and the diffracted k vectors that occupy the curved spherical manifolds. Within these arcs, Bragg-matched diffraction of the reference beam off of the noise grating manifolds increases the magnitude of the noise. This larger scattered intensity writes new and stronger noise holograms and the process repeats, amplifying the scatter in a distinct pattern. The figure shows only half of the possible noise grating manifolds for clarity. Interreference of the reference beam with scatter follows the same process, mirrored around the
xz plane.
While this analysis and the previous studies explain the angular distribution of recorded scatter patterns, they do little to guide materials and process development that would maximize SNR. The amplitude of the noise gratings is expected to be strongly dependent on the statics of the scatter centers, the dynamics of the recording process, Bragg selectivity of the noise gratings, and recording parameters such as period and intensity. These are difficult to study experimentally since it is challenging to know or precisely control the details of the initial weak scattering. Additionally, experimentally exploring the parameter space of materials formulations and processing conditions is costly when it is not clear what variables should be probed and in what ranges they would be expected to impact the SNR.
The goal of this paper is to therefore to establish a numerical approach to study these questions. Because the observed Bragg-matched noise gratings have structure in both transverse dimensions and the amplification process depends on growth in both depth and time, such a study requires a simulation in three spatial dimensions and time. Conversely, numerical studies of predominantly local effects such as phase segregation in holographic polymers have typically been restricted to a slice through the hologram in the plane of the grating wave vector and time [
5].
For generality, we will simulate a specific model material broadly similar to widely available commercial formulations [
6]. The goal is to capture the minimal materials and processing physics in order to establish scaling relationships and elucidate possibly fruitful approaches to experimental reduction of haze. In the following sections, we illustrate typical haze dynamics and angular dependence in a model material. We then describe a multi-physics code suite created to simulate signal and noise growth in this material. The code is validated by comparing to a closed-form solution to two-beam grating growth via a harmonic expansion solution to the chemical reaction/diffusion equations. Finally, we explore noise dynamics with this tool to make predictions about noise growth versus experimentally accessible parameters.
2. Methods
To study the how coherent light illuminating a scattering, photosensitive medium drives the development of both signal and noise, we implement a multi-physics code suite. The technique models the holographic recording process by iterating four primary calculations as a function of exposure time. These are (1) the Fourier beam-propagation (FBPM) method used to calculate optical propagation and diffraction within the 3D index structure, (2) a discrete time-stepping solution of the first-order photochemical rate equations, (3) diffusion of mobile species via a fast-Fourier transform (FFT) solution of Fick’s law and (4) update of the material index via an a priori model for two-component holographic photopolymers. We describe each step in more detail below.
2.1. Sample Preparation Including Scatter Centers
The substrates and recording material are assumed to contain photo-insensitive scattering centers that initiate noise gratings. These may be scratches in the substrate, dust, or incompletely mixed components of the formulation, for example. These scatter centers can be described either by shape (e.g., microspheres) or via their angular scatter distribution (e.g., Rayleigh). In the latter case, the microscopic shape of the scatter center sampled on the computational grid is found by inverse Fourier transforming the specified angular scatter distribution. Identical copies of these shapes are randomly distributed throughout the 3D recording medium and/or substrates with a specified number density. Given the shape and number density, the cross-section of the scatter centers is adjusted to create a specified initial haze.
Figure 2a shows a representative cross section of the polymer refractive index illustrating the randomly distributed scattering particles before any exposure.
2.2. Holographic Exposure
For each time step of the exposure, two coherent Gaussian writing beams, symmetrically tilted relative to the surface normal, are propagated through the depth of the material using the non-paraxial FBPM method. The x and y spatial dimensions are at least Nyquist sampled relative to the optical wavelength in the material in order to capture scattering angles over the full forward hemisphere. The writing beam angle is chosen to yield a transmission hologram with period,
sufficiently larger than the grid size such that individual grating fringes are adequately sampled. The writing beams attenuate due to initiator absorption, coherently interfere, scatter off of both the static microparticles and dynamically recorded noise gratings, and diffract off of the developing hologram. The three-dimensional optical intensity at each time step is calculated from the magnitude squared of the electric field, as shown in
Figure 2b.
2.3. Photochemical Reactions
This intensity then drives a sequence of chemical reactions, calculated through time-stepping the first order chemical rate equations. The first reaction in the sequence is the photolysis of an initiator (assumed to be present in sufficient excess that it does not significantly deplete) by the optical intensity to generate primary radicals with rate constant
, which rapidly react to form a radical monomer
. These radicals diffusive with diffusivity
D and terminate via a unimolecular reaction with rate constant
. We chose a unimolecular termination reaction to match the intensity reciprocity observed for two-component holographic photopolymers [
8]. To our knowledge, the linear and reciprocal intensity response of these materials has not been adequately explained in the literature and is surprising given the otherwise universal sublinear intensity response of radically initiated photopolymers. With these assumptions, the radical concentration evolves as:
The termination rate constant
is assumed to encompass any first order process which prevents radicals from participating in the subsequent propagation reaction. In steady state, the solution to this equation under sinusoidal illumination,
, is found to be a sinusoid with finite visibility:
where
is the average intensity,
is the grating wave number and the radical visibility is given by:
where
is the diffusion rate of radicals across a period relative to radical termination rate. As shown in
Figure 3, the distribution of radicals is identical to the sinusoidal intensity (
) when this unitless ratio is much smaller than 1. However, when diffusion is sufficiently fast that radicals can diffuse into the dark region of the fringe before terminating
, the visibility of the radical distribution decreases inversely proportional to period squared. This behavior has previously been identified as the resolution limit in published materials [
9]. This explains the observation that significant improvement in dynamic range can be achieved by increasing the termination rate via radical traps on the matrix [
10]. Here we will assume the steady-state distribution of radicals with finite visibility given by Equation (3), governed by the unitless diffusion rate to termination rate ratio
.
These radicals then drive the propagation reaction with rate constant
, depleting the mobile monomer as:
The concentration of polymer
is found by the time integral of the final term in the Equation (4). This assumes that all polymer is immobilized, e.g., by hydrogen bonding, entanglement or an explicit reactive group on the matrix. Oligomeric species are not tracked, simplifying the description of the propagation reaction to a single rate equation. Note that polymer concentration is defined to be the concentration of monomer that has reacted such that, neglecting diffusion, the sum of monomer and polymer is conserved as
, where
is the formulated monomer concentration. Extensions of this approach to more fully describe the polymer development have been extensively studied [
11]. The goal here is to capture the essential chemical physics of the material with the minimal description that will reveal scaling relations governing the dynamics of recording noise. Details of the time-stepping scheme for Equations (1) and (4) are given in
Appendix A.
2.4. Monomer Diffusion
Mobile monomer must diffuse to implement the first term in Equation (4). The three-dimensional distribution of monomer is Fourier transformed via the FFT and multiplied by the transfer function of Fick’s second law in homogenous space, then inverse transformed. Assuming that diffusivity, D, is independent of space and time, that transfer function is where are the Cartesian spatial frequencies and the characteristic spatial frequency of diffusion in a single time step, , is given by . Note that enforces the conservation of mass.
2.5. Index Development
The index of the photopolymer is then calculated from these chemical distributions. The local volume fraction of polymer is first found from the fractional conversion and the known polymerization shrinkage. The volume fraction of the rubbery solid matrix is found assuming perfect mixing and incompressibility as one minus the monomer and polymer volume fractions. Finally, the index of the three-component mixture is calculated by the Lorentz–Lorenz equation:
where the refractive indices
of the monomer, polymer and matrix, respectively, are found using prism-coupler measurements of mixtures of the components. These are typically not the same as measurements of pure compounds. For example, vitrification typically limits complete shrinkage upon polymerization. Conversely, in the rubbery host, the writing monomer can typically fully convert with minimal shrinkage stress. This yields a greater shrinkage and refractive index than exhibited by neat polymer at full conversion. Taking the fractional volume shrinkage upon polymerization to be
, the associated volume fractions of each component are:
where
is the monomer volume fraction before exposure and
and
are the concentration of monomer and polymer normalized to the initial monomer concentration. A cross-section of the typical resulting index is shown in
Figure 2c.
An important aspect of this calculation is that index change and shrinkage upon polymerization result in finite index contrast immediately upon polymerization. Subsequent diffusion amplifies this response. Since the development of haze is hypothesized to be self-reinforced noise gratings that form during exposure, the time dynamics of the index formation process are critical.
2.6. Holographic and Haze Readout
Finally, the efficiency of the hologram is found by illuminating the material index. This index is taken to be the sum of the static (scattering) and dynamic (holographic photopolymer) index distributions. A Bragg matched reference beam incident on the surface z = 0 is propagated through the volume with FBPM. The diffracted intensity is integrated at the exit face in an angular range centered on the object beam. The readout reference beam at 405 or 660 nm is assumed to be sufficiently weak and/or outside the initiator absorption spectrum that it does not modify the index distribution. Next, a normally incident beam is used to measure the haze. Haze is calculated, per ASTM standard D1003, as the fraction of the power scattered by the sample out of a normally incident probe for all angles (measured from the normal) greater than
and less than
. Typical angular spectra of the intensity exiting the material in both cases are illustrated in
Figure 4.
2.7. Post Exposure Development
Finally, after the exposure is complete, monomer is allowed to diffuse, becoming spatially uniform. This step simulates dark development time and assumes no limits to monomer transport due to finite solubility in high conversion regions. The uniform monomer is then fully converted to polymer, simulating a perfectly incoherent and spatially uniform flood cure. Quantitative conversion is assumed, that is no monomer remains. A Bragg selectivity curve is then calculated by measuring the diffraction efficiency as a function of reference beam angle.
3. Validation
To ensure the appropriateness and proper implementation of this combination of methods, it is critical that the code suite be validated. There are limited theoretical methods for predicting scatter development in these materials and the problem is inherently statistical. Thus, we instead validate the technique using a fully theoretical model of two-beam Bragg holography. A comparison to theory is chosen rather than experiment to enable direct comparison with no fitting to find parameters that are not always easily accessible experimentally, e.g., diffusion and reaction rates.
Combining Equations (2) and (4) results in a single rate equation governing this simplified system. Restricting the analysis to one dimension in the plane of incidence gives:
where
is the average polymerization rate. Equation (7) can be normalized with the substitutions:
where
is the formulated monomer concentration. This transforms Equation (5) into a unitless form:
which admits a standard harmonic solution
, here expanded up to 4th order [
12]. The resulting set of coupled first-order ordinary differential equations is:
These equations describe how the sinusoidal radical distribution of finite visibility given by drives the evolution of monomer concentration as a sum of grating harmonics. Consider the ideal recording case where diffusion is much faster than propagation such that . The equations show that all harmonics of the monomer higher than zero would remain identically zero, consistent with fast diffusion enforcing a spatially uniform monomer distribution. The average monomer concentration, , would then decay exponentially like which is the ideal first order reaction with time constant where is the average rate of polymerization.
Equation (10) is a sufficiently simple set of coupled first-order differential equations that it can be solved in closed form with the help of a symbolic mathematics package. With the monomer distribution
known, the polymer distribution can then be found from the time integral of the last term in Equation (9):
where
is the polymer concentration normalized to
.
represents the local concentration of monomer units that have reacted to become polymer. The left-hand side of Equation (11) is harmonically expanded as
and terms for each harmonic
are set equal, generating integral equations for each polymer harmonic
. These can also be solved in closed form with a symbolic math package. The resulting solutions to Equation (10) for the time evolution of the monomer harmonics,
and the solutions to Equation (11) for the time evolution of the polymer harmonics,
, are a complete description of how the material responds to a sinusoidal intensity pattern. In the ideal recording case of perfect radical visibility,
, and very fast monomer diffusion,
, the all polymer harmonics greater than two are zero. The average and first harmonic of the polymer become
and, thus,
. This is the ideal recording case where the polymer records a perfectly copy of the intensity.
With the monomer and polymer distributions known in space and time, the volume fraction of the matrix,
, can be found from the assumption of incompressibility,
. This completely specifies the composition of the trinary mixture. The composition and refractive indices of the constituents are then substituted into Equation (5), yielding a complicated but completely closed-form expression for material index as a function of exposure time. This lengthy equation depends on measured material parameters such as refractive indices, diffusion and reaction rates through the terms
and
, hologram period and time. The first harmonic of this calculated index is then used to predict holographic efficiency as a function of reference beam angle via the beta value method [
13]. The result is a completely closed-form prediction of hologram efficiency as a function of material parameters, writing conditions, and playback conditions.
The reason for developing this fully theoretical expression for holographic efficiency is that it can be compared to the numerical code suite which ideally implements the same physics but in a completely numerical, discrete manner. Specifically, in the limit of no scattering, the code suite will reproduce the intensity pattern, chemical reactions, distribution of species and refractive index to predict Bragg diffraction efficiency via the numerical techniques described previously. A comparison of the numerical to the theoretical solutions thus serves to validate the implementation and numerical approximations of the code suite.
First, we use both the theory and numerical system to compare the evolution of the index grating itself. A 10 µm thick film of 20 wt% bisphenol A diacrylate in polyurethane matrix is exposed to two 405 nm Gaussian beams tilted to produce an unslanted 1 µm period transmission hologram. Radical diffusion is taken to be negligible
but monomer diffusivity is finite and sufficiently slow
to generate a significant second, but limited third and higher, harmonics. The normalized exposure time,
, is sufficiently long to consume over 90% of the monomer. The solid lines in
Figure 5c show the evolution of the fundamental harmonic recorded into the polymer,
, as well as the second harmonic,
, caused by the finite rate of monomer diffusion. The open circles show the fully numerical result, extracted by fitting a harmonic series to a slice through the 3D sampled grid. The excellent agreement between the closed form 1 + 1D and numerical 3 + 1D approaches validates that the numerical suite has faithfully reproduced the desired physics. Next, both the theoretical and numerically derived first harmonics of the polymer and monomer are combined through the Lorentz–Lorenz equation to predict the growth of the first harmonic of the refractive index, as shown in
Figure 5d. The solid line is a completely closed form prediction of grating strength and the numerical suite accurately reproduces this prediction.
To further validate that the code correctly implements volume diffraction of off the numerically calculated index structure, a Bragg detuning curve found from both methods is compared. In the theoretical case, the calculated first harmonic of the refractive index is used in the beta value method [
13] to predict the efficiency as a function of incident beam angle relative to the initial reference beam angle. In the numerical case, the three-dimensional space is illuminated with a tilted reference beam and the electric fields at the far side of the material are calculated via the FBPM method. The diffracted power is found by integrating the intensity in a range of spatial frequencies centered on the diffracted object beam. This is repeated for a series of reference angles. To compare the beta value method prediction of plane-wave diffraction off of an infinite grating to the numerical calculation of a Gaussian shaped hologram read with a Gaussian shaped reference beam, the theoretical efficiency is divided by
where
is the 1/
e radius of the electric field [
14].
Figure 6 shows that the two approaches agree extremely well.
4. Results
The validated code is now applied to investigate how noise and signal develop in a holographic photopolymer containing scatter centers. The holographic photopolymer uses bisphenol A diacrylate (BPADA) as a writing monomer in a polyurethane matrix initiated with 2,4,6-trimethylbenzoyl diphenylphosphine oxide (TPO) at 405 nm. At 633 nm, the refractive indices of the monomer, polymer and matrix are 1.545, 1.5789, and 1.476 respectively. Dispersion of the indices is not included in this study. The shrinkage of the monomer in matrix upon full polymerization is found via the slanted grating technique [
15] to be
, consistent with typical vales of 22–24 mL/mol for complete shrinkage of acrylate monomers. Initiator concentration is sufficiently low that the samples may be considered optically thin and, thus, no absorption is included. Rather than use measured diffusivity, reaction rates, initiator concentration and intensity to calculate a specific value for the diffusion to reaction rate ratios for monomer and radicals,
and
, we will take these to be free variables in order to explore their role in haze development. Note that these ratios are not fixed even for a specific material since they depend on hologram period and optical intensity, thus, they can be varied experimentally to exploit the trends shown next. Similarly, we take initial haze to be a free variable to understand its importance in haze development.
The writing beams are tilted, except where stated, to produce a grating period of one µm and have a 1/e2 intensity diameter of 200 µm. The readout beam used for nondestructive measurement of diffraction efficiency and haze is at 660 nm with a 40 µm 1/e2 intensity diameter, tilted at the appropriate Bragg matched angle or normally incident, respectively. Initial haze is generated by Rayleigh scatter centers with number density whose cross-section is adjusted to give the desired initial haze.
Typical growth dynamics are shown in
Figure 7. Even for this case with relatively high initial haze of nearly 0.1, haze is observed to grow much more slowly than diffraction efficiency. In this case, the hologram is overmodulated before significant haze growth occurs. This observation is consistent with the hypothesis that noise gratings are written by interference of writing beams and scattered light while the intended grating is written by the much stronger interference of reference and object beams. The first result of this numerical technique is therefore that significant haze growth arising from amplification of Bragg-matched noise gratings only occurs for holograms that are strongly saturated.
Haze growth via amplification of Bragg-matched noise gratings is further supported by recording in diffusion and reaction-limited materials as shown in
Figure 8. The distinct patterns visible in the spatial frequency spectrum of the scatter in the lower plots are consistent with Bragg degenerate noise gratings amplified during writing by fast diffusion. Conversely, when no diffusion occurs during writing, the haze growth is limited and shows no distinct patterns. This also serves to illustrate that haze growth in two-beam holography is distinctly different than self-written waveguides [
3] in which scatter centers seed filaments whose scattering would be approximately isotropic in the transverse plane.
The top and bottom plots of
Figure 8 reveal how recording dynamics impact grating and scatter development. In the top diffusion-limited regime of
, negligible haze has developed because the feedback process whereby noise gratings seed the growth of additional gratings has been suppressed. Simultaneously, the single strong hologram has been limited in its diffraction efficiency because monomer was unable to replenish that depleted in the bright regions of the fringe. Conversely in the bottom reaction-limited case of
, the experimentally observed scatter patterns (see e.g.,
Figure 1) are observed. Diffraction efficiency is in the over-modulated regime since nearly all the monomer has been consumed due to rapid transport into the bright fringes. A multiplexed set of weaker holograms with appropriate scheduling would yield the same total usable index contrast in both cases, but the upper case would exhibit greater SNR.
The significant dependence of total haze upon the diffusion to reaction rate ratio shown in
Figure 8 indicates
has a critical role in haze development.
Figure 9 plots the final haze versus
under the same conditions as
Figure 8. The plots show two distinct behaviors. When diffusion across a period is sufficiently slow that noise gratings cannot develop during the exposure, haze is limited to that which is recorded by the instantaneous index change. Conversely, when diffusion is sufficiently fast, noise gratings are self-reinforced during the exposure. The transition between these two regimes is seen to be approximately the boundary between diffusion-and reaction-limited systems at
. Critically, this indicates that the relevant spatial scale for scatter is approximately equal to the period,
, consistent with the hypothesis that haze is not spatially random but is dominated by growth of Bragg degenerate noise gratings whose scale is similar to the grating period. This is contrast to the alternative hypothesis that the spatial scale of haze is similar to the particle size. Since scatter centers are typically much smaller than the period, noise dynamics controlled by diffusion would be much faster than signal dynamics in this case. The important implication for recording high SNR holograms is that both signal and noise obey similar dynamics such that diffusion-limited grating development will largely suppress haze development. The material and process variables that influence this condition are related by the single parameter
.
Since haze is seeded by scattering, an obvious strategy for reduction of haze is to reduce the initial scatter, e.g., via resin filtering and clean-room processing. Given that recorded haze is an amplified growth process, it is possible that reductions in initial haze could have super-linear impact on final haze.
Figure 10 shows how recorded diffraction efficiency and final haze scale with initial haze, controlled here via the cross section of Rayleigh scatter centers of the same number density of
. These results show that final haze is a fixed multiple of the initial haze. The ratio depends on the material and recording conditions. This result indicates that clean sample preparation is important, but has only linear impact on the final haze value.
Finally, we consider the dependence of recorded haze on sample thickness. It is commonly observed experimentally that haze grows rapidly with sample thickness and this is confirmed in
Figure 11 where the initial haze has been subtracted to highlight the dynamics of haze growth caused by the medium. The left and right plots illustrate how uniformly distributed scatter (left) produces less haze than the same scatter concentrated on the substrate through which the light enters (right). Conversely, the distributed scatter (left) exhibits greater nonlinearity versus thickness, consistent with noise gratings being written by scattered intensity which grows linearly in depth (left) but is initially constant in depth (right). While we have fit this data to cubic polynomials to guide the eye, the thickness dependence of the haze does not appear to have any simple scaling relationship. This may be due to the complexity of the writing intensity pattern for haze which is generated by scatter from the static and dynamic haze interfering with the reference and object beams whose intensity varies in depth as they couple. Once significant monomer has been consumed (
), the haze and the hologram also compete for monomer, slowing haze growth. This highlights the fact that complete understanding of haze dynamics requires a theoretical treatment. The numerical results and approach presented here may provide input to such an approach.