Photoelastic techniques are used to make both qualitative and quantitative measurements of the forces within idealized granular materials. The method is based on placing a birefringent granular material between a pair of polarizing filters, so that each region of the material rotates the polarization of light according to the amount of local stress. In this review paper, we summarize the past work using the technique, describe the optics underlying the technique, and illustrate how it can be used to quantitatively determine the vector contact forces between particles in a 2D granular system. We provide a description of software resources available to perform this task, as well as key techniques and resources for building an experimental apparatus.
I. INTRODUCTION
Photoelasticity has long been used by engineers to quantify the internal stresses within solid bodies.^{1} The technique is based on the understanding, dating to Maxwell,^{2} that light can be polarized (the electric and magnetic fields have a welldefined orientation) and its speed depends on the medium’s index of refraction. Many crystalline materials are birefringent: their crystalline axes specify a fast and slow direction for the propagation of light.
Noncrystalline materials such as polymers and glasses can nonetheless be photoelastic, with the degree of birefringence at each point in the material depending on the local stress. When a photoelastic material is placed between two polarizers and subjected to stress, each region of the material rotates the polarization of light according to the amount of local stress and the stressoptic coefficient of the material. This creates a visual pattern of alternating bright and dark fringes within the material, and this pattern depends in a detailed way on the orientation of the polarizers, the shape of the material, and how it is stressed.
Before the creation of finiteelement computing methods, photoelastic analysis was a key technique for predicting regions of high compressive or tensile stress during the design of manufactured parts. As such, it has been an important part of engineering practice in both its qualitative and quantitative forms.^{3–6} It was therefore inevitable that photoelastic methods eventually be applied to the study of idealized conceptions of granular materials, as was first done by Wakabayashi using epoxy tracer particles in 1950s,^{7} followed by Dantu (using ordered/disordered arrays of glass cylinders^{8}) and later by Drescher and de Josselin de Jong^{9} (in polymer disks) and Liu et al.^{10} (in glass spheres). The networklike structures visible in Figs. 1 and 2 have come to be known as force chains.
These pioneering experiments played a key role in the physics community’s realization that heterogeneous force transmission identifiable at the particle scale is of great importance for understanding the mechanics of granular materials.^{10,11} A series of papers by Behringer and coworkers established the photoelastic method as the key technique for making quantitative advances.^{12–14} In the past two decades, many different experimental teams have used photoelastic force visualization to quantify myriad phenomena. For example, it has been possible to reveal the particlescale anisotropy of the contact force networks,^{15} identify interparticle contacts,^{16,17} measure fabric and stress tensors,^{18} examine particle shape dependence,^{19} test the validity of statistical ensembles,^{20} identify dilatancysoftening,^{21} measure force chain order parameters,^{22} evaluate the grainscale stresses caused by growing plant roots,^{23,24} observe the effects of fluid flow,^{25} follow slow dynamics under shear,^{13,26–28} and quantify fast dynamics.^{29–32} Because the photoelastic response is known analytically for circles/ellipses, several of these studies have measured the vector contact forces at every interparticle contact within the granular material.^{33–35}
This paper first provides a review of photoelasticity and how it can be used to create highquality images of stress within the granular material (Section II). These images can be used either semiquantitatively (Section III) or quantitatively (Sections IV and V), with the latter providing a method to determine the vector contact forces for each of many circular particles in the system. We provide a description of software resources available to perform this task. Finally, we describe key techniques for building an experimental apparatus (Section VI), including resources for the purchase of key components.
II. Photoelasticity
This section provides a brief overview of polarized light and photoelasticity, to aid the reader in building a qualitative understanding of the method. A standard optics textbook (such as Hecht^{36}) can provide quantitative details about polarized light, and the classic twovolume monograph by Frocht^{1} provides a comprehensive treatment of photoelasticity as a method. Several recent book chapters also provide an introduction to the method’s use in the context of granular materials.^{35,37}
A. Polarized light
Electromagnetic waves travel in a direction coperpendicular to the sinusoidal electric and magnetic fields of which they are comprised (see Fig. 2(a)). Linearly polarized light consists of waves that all have the same orientation of their electric fields. Two common ways to generate linearly polarized light are reflection from a metallic surface and transmission through a polarizing filter (also known by a brand name, Polaroid^{TM}). These filters consist of a thin sheet of iodineimpregnated polyvinyl alcohol, with polymer chains preferentially aligned along one axis. The mechanism through which they polarize the light is that the iodine provides conduction electrons which polarize along the polymer chains and thereby cancel the electric field of the incident light polarized in the parallel direction but not in the perpendicular direction. Typically efficiencies are to absorb >95% of light in the parallel direction and transmit <40% of light in the perpendicular direction.^{38}
Circularly polarized light consists of electric and magnetic fields that remain mutually perpendicular but rotate with constant magnitude with lefthanded (counterclockwise) or righthanded chirality (clockwise) with respect to the direction of propagation. To create circularly polarized light, light is first linearly polarized and then passes through a quarter wave plate that advances the phase of one polarization with respect to the other. A quarter wave plate is a birefringent material of chosen orientation (fast/slow axes) and thickness matched to the wavelength of interest. The orientation of the fast axis at $ \xb1 45 \xb0 $ (see Fig. 2) with respect to the linear polarization determines whether the output is leftcircularly polarized (LCP) or rightcircularly polarized (RCP). Typically, the two optical elements are purchased as a single, fused unit. Note that, unlike linear polarizers, a circular polarizer can be arbitrarily rotated by any angle around the transmission axis with no effect on the filtered light but cannot be flipped over.
B. Polariscopes
The basic principle of a darkfield polariscope is to place a photoelastic sample between a polarizer and an analyzer of the opposite polarization (e.g., vertical and horizontal linear polarizers or LCP and RCP circular polarizers), viewed in transmission. Sample configurations are shown in Fig. 4. Note that care must be taken during the installation of circular polarizers, to ensure that the polarizer and analyzer are both placed with their quarterwave plates facing the sample.
In each case, initially unpolarized light becomes polarized after passing through the polarizer. As the light passes through the sample, the amount of birefringence at each point depends on the local stress state (photoelasticity). This birefringence causes the light to rotate its polarization as it passes through the sample, and when it exits the sample the polarization has been rotated by a total of θ degrees. The light now passes through a final polarizer called the analyzer. The final intensity of light depends on the local stress in the birefringent sample as well as the orientation of the analyzer. In a stressed sample, this leads to an image of light and dark fringes (see Fig. 3).
While both linear and circular polariscopes are possible, there is a practical difficulty with linear polariscopes that limits their utility: In a linear polariscope, there are two possible reasons for fringes to appear: (1) the photoelastic response of the sample rotates the polarization to match/mismatch the analyzer; and (2) if the principal stress of the sample happens to be aligned with the polarizers, it will not be affected by the photoelastic response of the material. Therefore, for linear polariscopes, the fringe pattern depends on both the relative angles between the polarizer (fixed) and the principal stress (changes at different locations and times during experiment). For this reason, quantitative work with polariscopes should always be conducted with circularly polarized light, which provides azimuthal symmetry.
While the resulting image of light and dark fringes depends uniquely on the stress in the sample, there is an ambiguity in the phase. Various groups have proposed various methods to unwrap the phase map of the resulting image to extract the stress in the sample. However, many images of the sample are needed with different wavelengths of light^{39} and different analyzer orientations.^{40} These methods must resolve both isochromatics (associated with the magnitude of local stress) and the isoclinics (the angle of stress). In investigations of granular materials, we are solely interested in the contact forces of each sample rather than the local stress state itself. Therefore, it is sufficient for contactforce measurements to use an analyzer which has been arbitrarily rotated about its azimuthal axis (with respect to the polarizer); this results in an image of the isochromatics only. Throughout the remainder of the text, we assume images are acquired using the circular polarizers set up as shown in Fig. 4.
Generally speaking, the samples in Figs. 1 and 2 illustrate that particles with more force are brighter (darker) against a dark (bright) background because the light that passed through them had its polarization rotated. Often, this contrasting effect is sufficient to identify general features of stresstransmission within a granular sample. Methods for quantifying these features are described in Section III. Importantly, there is a quantitative, onetoone mapping between the stress field in the sample and the resulting pattern of fringes.^{1} Therefore, it is possible to solve the inverse problem and use polariscope images to determine the vector contact forces on each particle in the sample. Quantitative methods for performing this inversion are detailed in Sections IV and V.
III. SEMIQUANTITATIVE MEASUREMENTS
Originally, images of the type shown in Figs. 1 and 2 were used for semiquantitative investigations, without the knowledge of the vector contact forces between the particles. Even in modern experiments, there is much that can be learned from comparing the spatial patterns of the force chains. For example, images of the force chains can reveal history dependence via spatial correlations,^{41,42} the ensembleaveraged response to point forces,^{43} and the principal axis of the response.^{44} By subtracting images taken before/after an event, it is possible to examine spatial patterns of failure^{26,45} and the speed of sound.^{31}
To quantify how the number of dark/light fringes increases with the stress on each particle, Behringer and coworkers^{11,13,14,46} calculated the squared local gradient, averaged over all pixels in an image
They determined (see Fig. 5) that this quantity provides an empirical measure of the 2D stress on the system. More recently, this has come to be known as the gradientsquared (G^{2}) method and has been used on the systemscale, chainscale, particlescale, and contactscale to provide a semiquantitative way to track changes in granular forces.^{32,47–49}
The G^{2} method has several important advantages over the photoelastic inverse methods to be described in Section V. First, it works on lowresolution images (10 pixels/particle is enough) as long as there is sufficient intensity contrast. Second, image processing times are very fast: at the time of writing, it takes 0.05 s to calculate the pressure within a 1 megapixel image on a single processor. Finally, there are no parameters in Eq. (1). To translate G^{2} into pressure, only a single calibration experiment needs to be done under the same lighting conditions as the experiment, in order to properly account for the image intensity and contrast. For these reasons, the G^{2} method has remained popular, even as the speed and accuracy of vector contact force measurement codes have improved.
IV. PHOTOELASTIC THEORY
As illustrated in Fig. 6, a set of z contact forces on a disk determines the pattern of photoelastic fringes. Because the analytical solutions for the stress created by point contact forces on a disk are known,^{1} it is possible to determine the intensity field $ I ( x , y ) $ from the set of z vector forces. In analyzing images from experiments (to be shown in Section V), we will need to do the inverse: for a known $ I ( x , y ) $ recorded by a camera, determine the vector contact forces which created that image. In this section, we present the theory required to determine the mapping from $ f \u2192 \u2192 I $ and calibrate the measurements so that they can be used for the inverse problem in Sec. V.
A. Solution for diametric loading
The special case of diametric loading provides both an instructive example for understanding quantitative photoelasticity and the means of calibration necessary for any quantitative experiment. As illustrated in Fig. 7, the number of dark/light fringes increases with the diametric load of magnitude f (equal and opposite forces along the diameter of the disk). We can understand this behavior starting from the known solutions for the stress tensor at the center of the disk,^{1,50}
The intensity of the fringe pattern at any point is given by
where $ ( \sigma 1 \u2212 \sigma 2 ) $ is the difference in the principal stress at the point (x, y), h is the material thickness, λ is the wavelength of light, and C is the stressoptic coefficient (a λdependent material property).
The spatial variations in I, visible in Fig. 8, are present because $ \sigma 1 \u2212 \sigma 2 $ is spatially varying. We will consider the simpler case of monochromatic light, so that $ F \sigma = \lambda / C h $ (the stressoptic coefficient) is a constant. In the case of isochromatic fringes, Eq. (5) reduces to
Using Eq. (4), we have $ \sigma y y = \u2212 3 \sigma x x $, where the offdiagonal components $ \sigma x y $ and $ \sigma y x $ are zero and the shear stress is zero at the center of the disc. The difference in the principal stress is the difference in the eigenvalues, which then is given by $ \sigma 1 \u2212 \sigma 2 =  \sigma x x \u2212 \sigma y y  $. Then, at the center of the disc, the number of fringes N_{fr} is given by
For an unstressed disk, $ \sigma 1 = \sigma 2 = 0 $ everywhere in the disk, the difference in the phase is zero and thus the isochromatic image is dark. The first bright fringe occurs when the difference in the principal stress is $ \sigma 1 \u2212 \sigma 2 = F \sigma / 2 $. The bright and dark fringes on the disk arise from the periodic function $ ( sin 2 x ) $ in the phase retardation. As shown in Fig. 7, the total number of fringes increases as the sin function cycles through an increasing number of bright and dark fringes for increasing values of $ \sigma 1 \u2212 \sigma 2 $.
Thus, $ F \sigma $ can be empirically determined by diametrically compressing the disk by a known force and recording the fringe number as follows. For diametric loading, the fringe number is defined as the number of fringes counted between the edge of the disk and the center of the disk, as shown in Fig. 8(a). The fringe number at the boundary is zero. Moving toward the center, each bright/dark band increments the fringe number by 0.5. In order to measure C (the material parameter), a typical experiment proceeds by placing a disk of known radius and thickness into a load cell and illuminating it with monochromatic, polarized light. As the diametric load on the particle is slowly increased, a measurement is taken each time a new fringe appears (N_{fr} increases). A scatter plot of the diametric force and number of fringes allows for a fit to the left and right sides of Eq. (7). Since all other parameters are known, this provides a measurement of C for that particular material, as shown in Fig. 8(b). Most of the commonly used materials (see Section VI A) have $ C \u2248 3000 $ Br.
B. Solution for arbitrary contacts
Ultimately, we will determine all contact forces acting on a particle based on the pattern of isochromatic fringes recorded in images of photoelastic disks. This goal is illustrated for a single particle in Fig. 6. To achieve this will require a formula for the solution I(x,y), as a function of an arbitrary set of z vector contact forces, acting at arbitrary contact locations around the perimeter of a disk. In this section, we derive the solution to the stress field $ \sigma ( x , y ) $ inside an elastic disk in mechanical equilibrium with z point forces on the boundary, along with a short review of the necessary theory of elasticity. The solution presented here largely follows the work by Refs. 1, 50, and 51.
For any linear elastic solid with boundary forces, there are three basic assumptions: (a) the condition of equilibrium (force and torque balance); (b) a constitutive relation (here, linear stressstrain); (c) SaintVenant’s condition of compatibility (no gaps or overlaps). In two dimensions and in the absence of body forces, the condition of equilibrium for a static solid is given by
Equivalently, $ \u2207 \u22c5 \sigma = 0 $. The second assumption is a linear relation between the strain and stress tensors,
The third assumption is SaintVenant’s condition: the solid deforms continuously leaving no gaps and creating no overlaps.^{52} For plane strain in the absence of body forces, the compatibility of infinitesimal strains is given by
The solution to Eq. (10) in two dimensions is solved using Airy stress functions, φ, where the components of the stress tensor are derivatives of φ,
Then the compatibility equation of the stress function is biharmonic
In general, the solution to Eq. (12) is guessed by using symmetries or adjusting a known solution to the boundary conditions.
Consider a stress on a disk, as shown in Fig. 9, with radius R and thickness h. We define the origin at point O_{1} where the direction of pressure is O_{1}O_{2}, with angles $ \theta 1 $ and r_{1} as shown in Fig. 9. Points on the circumference of the circle have a radial stress $ \sigma r r = \u2212 2 f \pi h cos \theta $ in the direction r_{1}. The angle between $ \sigma r r \u2032 $ and the tangent is $ \theta 2 $. The components of the stress tensor are
Substituting $ \sigma r r $ and using $ \theta = \theta 1 $ and $ r 1 = 2 R sin \theta 2 $ everywhere on the circumference, then
These two tensors are the superposition of three stresses at point A on the circumference as follows:

normal stress: $ \u2212 f \pi R h sin ( \theta 1 + \theta 2 ) $.

tangential stress: $ \u2212 f \pi R h cos ( \theta 1 + \theta 2 ) $.

uniform stress in the direction of r_{1} with normal and tangential components: $ \u2212 f \pi R h sin ( \theta 2 \u2212 \theta 1 ) $ and $ \u2212 f \pi R h cos ( \theta 2 \u2212 \theta 1 ) $, respectively.
If the disk is force and torque balanced, torque balance implies $ \u2211 i z f i cos ( \theta i , 1 + \theta i , 2 ) = 0 $, and the sum of the uniform tension must also be zero. Therefore, the remaining term is the sum of the normal tensions at the boundary $ \u2211 i z f i \pi R h sin ( \theta i , 1 + \theta i , 2 ) $ along the circumference, which is simplified using the inscribed angle theorem $ \theta 1 + \theta 2 = \pi 2 \xb1 \alpha $ and the identity $ sin ( \pi 2 \xb1 \alpha ) = cos ( \xb1 \alpha ) = cos \alpha $. To ensure that the boundary of the disk is stress free, a uniform tension of magnitude $ f \pi R h cos \alpha $ is added to the radial stress distribution. The solution for z concentrated forces on a disk where $ \theta i $ is the angle between fi and r_{i} and α as shown in Fig. 9 is
All other components of σ are zero.
The solution is in zpolar coordinates. To write the solution in a single coordinate system, we use a Cartesian coordinate system with the origin at the center of the disk. Let β be the angle between the positive Yaxis and the line segment $ C O 1 \xaf $, as shown in Fig. 9. We compute $ \sigma r r $ as in Eq. (15) and transforming coordinate systems with a rotation by $ \theta i = \u2212 ( \beta i + \alpha i ) $ as
where $ T ( \theta i ) $ is the rotation matrix.
The eigenvalues of σ are
where the subscripts 1 and 2 correspond to the + and − signs, respectively.
We have presented the solution to zforces on a disk and given a calibration of $ F \sigma $ for the photoelastic disk. Section V details the numerical methods involved in taking this solution and finding the contact forces in an array of disks.
V. AUTOMATED INVERSION
The ultimate goal of the forceinversion process is to start from images such as Fig. 10, find all of the particles, and determine the vector contact forces between them using the theoretical framework described in Section IV. While the equations for the image intensity due to a known set of particle positions and contact forces are known analytically, the inverse of the problem is not. Several research groups^{15,20,28,35} have created computational solutions to this inverse problem, including in C++^{53} and MATLAB,^{54} a technique pioneered more than a decade ago^{33} by Majmudar.
Here, we describe our most recent implementation, PEGS, inspired by the earlier work of Majmudar^{33} and Puckett.^{34} PEGS is written for the MATLAB platform and is available for download and development on the GitHub open source software platform.^{54} It makes use of Matlab’s builtin parallelization toolbox, making it suitable for running on highperformance computers to reduce computational time.
In general, it is necessary to have two images for each experimental configuration: one taken in unpolarized light for detecting particle locations, and one taken in polarized light for measuring the forces on those particles.^{16} Several possible ways to do this are discussed in Section VI D.
Throughout this section, we illustrate the process using the sample set of images shown in Fig. 10. The top image is presented as taken in RGB color by the camera. This comes from a setup illuminated by two monochromatic light sources: red unpolarized light and green circularly polarized light (see Fig. 16). This configuration allows the camera to automatically collect the two necessary images within a single color frame.
The PEGS inversion process comprises three programs. The preprocessor program detects particles, identifies neighbors, and validates which of those neighbors are contacts. This program operates on both the polarized and unpolarized images. The positions of the particles and contacts are then used by the fringeinverter program, which performs the imageinversion process by which the photoelastic fringe patterns from a circular polariscope are used to determine all vector contact forces. This process operates on polarizedlight images only. Finally, the postprocessor program contains a suite of tools for errorcorrection and visualization. This provides a means to check the output for quality and then use this feedback to adjust the program’s parameters before submitting batch jobs to analyze a full dataset.
A. Preprocessor
The first step is to obtain the locations and radii of all particles in the packing, from an image such as Fig. 10(b). From experience, we find that this needs to be done with an accuracy of $ \u2272 0.05 d $ in order to successfully invert the fringe pattern; the quality of the fringeinversion fits will improve significantly with an increased centroid accuracy. We have found that the MATLAB imfindcircles() function provided in the Image Processing Toolbox is able to achieve the necessary level of precision; other Hough transforms provide similar performance. Several algorithms are available within that function, and for the results shown below we used the Hough transform (see Ref. 55 for a review), applied separately for each particle size. The results shown in Fig. 11(a) use red/blue circles to mark the locations of all particles with large/small radii, respectively. Many other particlefinding routines, including convolution methods, could alternatively be used. For a review, see the work of Shattuck.^{35}
In order for the fringeinversion process to work efficiently and effectively, it is important to have a highquality list of forcebearing contacts. The initial list of candidatecontacts is created by first identifying all neighbors by a threshold distance between the centers of a pair of particles of known radii, located at $ x \u2192 i $ and $ x \u2192 j $. A naïve definition of a contact would be to call two particle contacts if their separation was within some tolerance d_{tol} of the sum of their respective radii,
However, even with subpixel resolution, this method is inadequate for determining whether the contacts are forcebearing. For each disk, the number of neighbors n satisfying Eq. (18) is greater than or equal to the actual number of contacts z, and d_{tol} must be set generously to overcome uncertainty in particle positions.
Therefore, we perform a secondary screening step to determine which of the neighbors should be included in the subset of z contacts. The criteria for inclusion are that the photoelastic response in an area adjacent to the point of contact must be above a threshold value.^{16,33} In the analysis that follows, we use a threshold based on the $ \u27e8 G 2 \u27e9 $ response (see Eq. (1)), rather than the intensity. If both particles show sufficient response near the contact point, then the contact is declared to be loadbearing.
A sample comparison is shown in Fig. 11, where all of the small white and green circles indicate possible contacts, but only the green circles meet the threshold criteria for being a contact. This significantly reduces the possibility, shown in Fig. 11(b), that a fringe located near the outer edge of one particle cause a false positive. This thresholding method is sensitive to having sufficient resolution and contrast. As a practical measure, it is necessary to adjust the contact threshold and the size of the evaluated area around the contact point for each specific dataset. In the PEGS software,^{54} we provide a verbose option in the program allowing the user to get a visualization of the detected contact network, as shown in Fig. 11(c).
In addition, uncertainty in which contacts should be included in z can be overcome after the fact: the forceinversion code can later determine that $ f i \u2248 0 $ (a null contact). This simply increases the parameter space over which the inversion process must search, from 2z − 3 to 2n − 3. Therefore, it is better to set d_{tol} and the G^{2}threshold to be generous in including likely contacts on the list. While this will result in an increased computational load, the force solver is free to set the force of any contact to zero, but it is not allowed to create new contacts.
B. Fringeinverter
Using the list of particle positions, radii, and contacts, we can iterate over all particles in the system to determine which set of $ f \u2192 = ( f , \alpha ) $ forces from the general solution (Section IV B, Refs. 1, 50, and 51) would produce the observed pattern of isochromatic fringes. In the text below, we use the same naming of the location of the contacts β, forces f, and impact angles α as used in Fig. 9.
The inversion process proceeds one disk at a time, for which we already know the locations of the z contacts from Section V A. Importantly, the location of all contacts must be along the line between the two particle centers and is specified by the polar angle β. This is satisfied by default in the algorithm described in Section V A.
We generate an initial guess for each force magnitude f by first applying the gradient squared method from Eq. (1) in Section III to the entire disk. We then distribute this total force among the z contacts in proportion to the value of G^{2} within a small area around each contact point. To make an initial guess for the angle α of each force, we select a value for which the given contact locations β and guessed values of f would collectively be in force balance (see Section V C). This is included as an optional step, since force balance will not be applicable for systems for which there is significant dynamics.
Using the set of z values of $ ( f , \alpha , \beta ) $, we solve for the stress field $ \sigma ( x , y ) $ on a grid chosen to match the resolution of the original image. Optionally, it is possible to downscale the grid for faster computing time at a reduced resolution. Together with Eq. (6), this creates a pseudoimage $ I ps ( f , \alpha , \beta ) $ of the isochromatic fringes. We compare the guess to the original image I_{obs} by the residual function
where x and y are the spatial coordinates of the images. This residual function is suitable for optimizing via any standard algorithm, with the values of $ f 1 \u2026 z , \alpha 1 \u2026 z $ are varied in order to find the minimum value for e, as was first done by Refs. 16 and 33. We have had good success using the LevenbergMarquardt^{56,57} implementation contained within the MATLAB lsqnonlin() function. Note, however, that this will only find a local minimum in the $ ( f , \alpha ) $ landscape and is therefore sensitive to the quality of the initial guesses. Because the stress field calculation is computationally expensive, we calculate the σ solutions for each gridpoint in parallel. This option, available natively in MATLAB, enables the code to run efficiently on multicore/threaded computers or high performance computing clusters.
A sample outcome for the image in Fig. 10 is shown in Fig. 12. The left column contains the original image, and the right column contains the pseudoimage resulting from the optimization procedure. Visually, the agreement is excellent; Section V D will show how to use Newton’s third law to provide a quantitative evaluation. Note that the number of fringes rapidly increases in the vicinity of strong contacts, as shown in the singledisk pseudoimages in Fig. 10. Because this information is not available for comparison in the observed images due to the optical quality of the particles, we find that the accuracy and efficiency of the results are improved by fitting only the innermost 0.95r of each disk rather than using all pixels up to the outermost edge.
C. Force and torque balance
For both the initial guesses (mentioned above), but also for all subsequent optimization steps, it is possible to improve the quality of the guesses by making use of force and torque balances. Because these additional constraints reduce the number of unknowns in $ ( f , \alpha ) $ from 2z to 2z − 3, we have found that this optional step greatly improves both the accuracy of the fit and the time to convergence.
For z contacts, the force and torque balance constraints are given by
For z = 2, the above equations reduce to
By examining Fig. 9, α can be determined geometrically by the internal angles in the triangle O_{1}O_{2}C. The magnitude of the contact force f is the only free parameter.
The initial test of the algorithm for minimizing $ e ( f , \alpha ) $ is used on a set of diametrically compressed disks, where I_{obs} is the image of the disk and f_{applied} is known. In Fig. 8, the result of the best fit $ I ( f , \alpha ) $ is compared with I_{obs}. For z = 2, initial guesses quite far from the true $ ( f , \alpha ) $ will still converge; however higher fringe numbers require better initial guesses. The error in the magnitude of contact force was typically $ \u2272 5 % $ and did not depend on the magnitude of the force (see Fig. 13).
The solution for $ z > 2 $ requires a bit more work. The force and torque constraints in Eq. (20) reduce the number of free parameters by 3. To solve for these constrained parameters, one could do a coordinate transformation from $ ( f , \alpha ) $ to the Cartesian (f_{x}, f_{y}); however this is not necessary. The solution requires a coordinate transform of the contact locations β by rotating each contact so that $ \beta z = 0 $ (i.e., taking $ \beta i = \beta i \u2212 \beta z $). After, $ ( f , \alpha ) $ are found, we simply rotate the coordinates back. Solving Eq. (20) for f_{z−1}, f_{z}, and $ \alpha z $, we have
This will ensure a set of $ f 1 \u2026 z $ and $ \alpha 1 \u2026 z $ that is both force and torque balanced.
D. Postprocessing
The forceinversion code provides the full set of $ f \u2192 = ( f , \alpha ) $ values for each contact, and the residual score e (Eq. (19)) for each disk. Using Newton’s third law, it is possible to quantitatively evaluate the reliability of the results shown in Fig. 12. To illustrate this on the small number of example particles, Fig. 13 compares the force fit on one side of a contact to the (presumably identical) force fit on the other side. Ideally, all contacts would have $ f A \u2192 B = f B \u2192 A $; therefore deviations from this line indicate errors in the inversion process. We observe that these are of a consistent relative magnitude throughout the range of forces observed.
In addition, since each contact is subject to two independent optimizations, it is possible to flag contacts with unreliable results. Individual values, or entire disks for which the residual score e is too high, can be automatically or manually replaced by the value from the oppositecontact. This feature also allows the data to be transformed so that all contacts (rather than all particles) experience forcebalance, by simply assigning each contact a force which is the vector average of the two independent measurements.
Visualizations such as Fig. 12 also provide feedback about improving the resolution and lighting. For example, the forceinversion technique is limited by resolution and contrast in I_{obs} and also governed by the number (and quality) of initial guesses provided to the algorithm. In the low force limit, I_{obs} is dark in the center of the disk and contacts appear as small dim spots near the perimeter. Greater image resolution and contrast can improve the accuracy of the contact forces found during optimization. Similarly, the resolution of the image also limits how well large forces are calculated, where fringes can become pixelated. Large forces (and large z) are more computationally expensive, and more/better initial guesses are required to find the correct local minimum for $ e ( f , \alpha ) $.
VI. BUILDING AN EXPERIMENT
A. Making particles
Several commercially available materials are available with a uniform photoelastic response. These generally come in various stiffness and are sold in one of the two main formats: sheets, or castable liquids. Methods suitable for both formats are discussed below. Because of the expense associated with all types of particles, it is advisable to make only a dozen particles in the first batch and check them between polarizers to determine whether the quality is good enough. For instance, frozenin stresses often take the form of a bright ring around the outer edge of the particle when viewed in a polariscope. It is sometimes possible to eliminate/reduce these stresses by slow, gentle heating over a few days or weeks.
In designing the particles, several geometrical considerations are important. The granular materials will be most amendable to experimentation when all particles (disks/cylinders) stay in a single plane. This can be achieved by making particles which are wider than they are tall so that they do not tip over when sheared. When calculating the amount of material needed, it is helpful to cut at least 25% more particles than seem necessary; particles inevitably become lost or damaged over time, and it is difficult to later make identical particles. These are then best replaced by particles made in the same batch as the others to avoid the necessity of particlespecific calibrations.
In humid climates, photoelastic particles will experience some degree of adhesion to themselves and their substrate due to liquid bridging. One method of reducing this effect is to dust the experiment with a very small amount of baking powder, a technique pioneered by the Behringer Lab.
1. Sheet cutting
The following are common sources of sheet stock. Vishay PhotoStress Plus^{58} is a proprietary material with excellent uniformity. Nearly clear urethane sheets are available from various industrial sources^{59} and are a cheaper alternative. Both of these options come in sheets of various stiffness and thickness and have optical clarity and uniformity suitable for quantitative measurements. Sheets are available with a small amount of dye, which can improve particletracking. Although clear cast acrylic (Plexiglas^{TM}) is widely available, it is generally too stiff (low $ F \sigma $) for most granular experiments. In addition, it is weakly dichroic (it exhibits polarizationdependent light absorption) and is therefore not suitable for quantitative photoelastic analysis.
For any sheet material, several cutting methods are possible depending on the technical capabilities of the machine shop. For circular particles, a skilled machinist can create a custom tool (see Fig. 14) which acts as a spinning disk cutter to create circular particles of a specified diameter and with little waste. Alternatively, a smallsized computercontrolled mill can trace the outside of individual particles and make repeated custom (noncircular) shapes. However, the sidemill tool will waste the material in proportion to its diameter. In either case, the machinist will need to lubricate the cutter with dish soap rather than the usual machine oil, which degrades the polymer. With good lubrication, it is possible to create particles with smooth edges all the way down, and no frozenin stresses. It is also possible to use a water jet cutter to make custom shapes with little waste if care is taken to minimize the presence of a notch or flap at the starting/ending location.
Two simple cutting options, which initially seem promising, are less practical than they might appear. The use of a punch is illadvised since these leave a lip at the bottom of each particle. Unfortunately, both PhotoStress and urethane break down due to the heat of a laser cutter and release toxic gases, so these do not work. For a lowcost option, it is instead advisable to simply use a band saw to cut strips and then particles as was done in early experiments.^{43}
2. Casting
Both Vishay PhotoStress Plus^{58} and urethane^{60,61} are also available as a twopart castable liquid. Both formulations are available in a variety of stiffness and can be clear or lightly dyed. The companies provide protocols for making and releasing molds. Note that in order to achieve bubblefree particles, it is necessary to cure the material under vacuum. Castable liquids have allowed for the creation of particles of unusual shapes and also with magnetic inclusions.^{62}
B. Polariscope configurations
Two main polariscope configurations are possible: the light can pass through the granular material once on its way to the camera (transmission) or twice (reflection). The most common configuration is the transmission polariscope, shown schematically in Fig. 15. If the experiment is horizontal, the upper plate could be omitted but is often included to prevent outofplane buckling of the granular layer. The pair of polarizers may be either of opposite chirality (darkfield, more common) or the same chirality (brightfield).
In principle, it does not matter whether the polarizers are inside or outside of the supporting layers. In practice, however, if the supporting layers are themselves made of a birefringent material (e.g., acrylic), then the polarizers must be placed directly adjacent to the particles to avoid measuring photoelastic effects present in the supporting layer. The camera and lightbox will be placed on either side of this sandwich. Such a configuration has the additional difficulty, however, that the polarizers will become scratched over time.
There are several advantages in arranging the system so that the cameraside polarizer is attached directly to the camera. First, this avoids the expensive purchase of a second, large polarizing filter. Second, easytouse polarizing filters are available for many commercial camera lenses (see Section VI C for details).
In some cases, a transmission polariscope is undesirable when the supporting layer must be opaque. This is the case, for instance, for a porous membrane (frit) providing airlevitation as in Ref. 20. To accommodate such situations, it is possible to build a reflection polariscope which requires optical access from only one side of the granular material. This is shown schematically in Fig. 16.
To create reflecting particles for use in a reflection polariscope, a simple method is to coat one side of the particles with silver spray paint (e.g., RustOleum Mirror Effect, used by the Behringer Lab). These particles will reflect the incident light and simultaneously flip its polarization. Therefore, a darkfield polariscope results from configuring both the illumination source and the camera with the same chirality polarizer.
C. Polarizers
Two formats of polarizers are typically available, large sheets and camerabased filters. Care must be taken to identify which polarization is specified and to select left/right circularly polarized products as desired. For large sheets, current sources include API Optics, http://www.apioptics.com/,^{38} Edmund Optics, http://www.edmundoptics.com/,^{63} and Alight Polarizers, http://www.polarization.com/.^{64} Photography shops sell camerabased filters that are sized and mounted for use on standard camera lenses.
Care must be taken when installing the filters in the polariscope, so that their orientation is correct. A simple check for a darkfield (lightfield) polariscope is to place them facetoface and rotate one polarizer with respect to the other. When the quarterwave plates are both facing inwards, rotating either polarizer will always maintain a dark (light) view. If either polarizer is flipped inside out with respect to the other, then the view will go alternately light and dark as one polarizer is rotated with respect to the other. Note that circular polarizers built for standard photographic purposes will have their quarterwave plate facing the camera sensor; this is the wrong orientation for a polariscope analyzer (see Fig. 2). To remedy this, it is usually possible to unscrew the mounting frame and flip the filter over within its holder so that the quarterwave plate faces outwards from the camera. Recall, also, that for a reflecting polariscope, having a polarizer and analyzer with matched chirality (LCP vs. RCP) gives a darkfield image. A final difficulty is that the specific chirality of photographic filters may or may not be provided at the time of purchase.
D. Lighting
Since the fringe pattern for birefringent particles depends on the wavelength of light through both the stressoptic coefficient and rotation of the light phase (see Eq. (5)), it is best to record monochromatic light for quantitative results. There are two main ways to do this, either by using a monochromatic light source or by filtering polychromatic light at the source or at the camera itself. A very simple method, suitable for nonquantitative studies, is to simply visualize only one of the three RGB (redgreenblue) channels of a fullcolor image, a technique pioneered by Ref. 16, by applying multiple filters in succession, and now possible to achieve using singleimages. Because these channels have a bandwidth of ∼100 nm, quantitative studies using this type of color separation are wise to still use monochromatic light. This could also be achieved via a monochromatic filter which screws directly onto the photographic lenses, available at any camera supply store.
Different lighting solutions are appropriate depending on whether the system is a transmission or reflection polariscope. In either case, it is better to achieve uniform lighting before collecting data rather than requiring a postprocessing step to equalize the brightness levels across all images. In all cases, it is preferable to use lowtemperature lights such as lightemitting diodes (LEDs) or fluorescent bulbs rather than halogens or incandescents which will heat the experiment and change the material properties of the particles and/or their supporting layers.
For a transmission polariscope, several commercial solutions are available. Common types of lightboxes include artists’ copyboards, sign and menu boards, and medical xray lightboxes. LEDbased boards are typically quite thin and have no AC flicker when running on rechargeable batteries. A doctor’s lightbox, while commonly available used at low cost, uses fluorescent lightbulbs which have significant 60 Hz flicker unless they have been specially designed to eliminate this or use an electronic ballast that upconverts the frequency to 20 kHz. It is also straightforward to build a custom lightbox out of lowcost strips of LEDs which can be cut to length to suit any shape experiment. For a reflection polariscope, many other options are possible, including colored LED spotlights (either decorative or theatrical) which provide monochromatic light.
To perform particletracking and photoelastic measurements at the same time, it will be necessary to have at least two different sources of light. Typically, a polarized, monochromatic light source provides the photoelastic image via a polarizing filter placed on the camera. A nonpolarized light source, which will still be predominantly transmitted through a cameramounted polarizing filter, provides other measurements such as the particle positions. There are numerous methods for achieving this goal: multiple cameras with different filters and postprocessing to align the images, electrically or mechanically triggering different lights in sequence; placing/removing one of the polarizing filters; or using RGB colorseparation of simultaneous sources within a single camera. This last option is illustrated in Fig. 10, for which the two images came from the red, green LED configuration shown schematically in Fig. 16.
This configuration leaves the third color channel (blue) available for supplemental use. For example, it is possible to draw features with an ultraviolet marker that will glow blue when exposed to black light. This technique has been used to monitor the rotation of particles,^{18,28} to identify and track a subpopulation of particles,^{20} and to outline the edges of very clear particles to improve their detection.
ACKNOWLEDGMENTS
K.E.D. is indebted to Bob Behringer, in whose lab she was first introduced to the use of photoelastic materials and who has championed the technique for many years. The photoelastic experiments by the three authors have been supported by the National Science Foundation (Nos. DMR0644743 and DMR1206808) and the James S. McDonnell Foundation. This review article began as a resource sheet and talk for the NSFfunded (No. EAR1537902) workshop “Analog Modeling of Tectonic Processes” held in Amherst, MA in May 2015 and was further developed for the Spring School “Imaging Particles” held in Erlangen, Germany in April 2016, funded by the German Science Foundation (DFG) through the Cluster of Excellence “Engineering of Advanced Materials” No. EXC315. The authors thank Arne te Nijenhuis for his work on Fig. 2, Amalia Thomas and Zhu Tang for their work on Fig. 8(b), and Amalia Thomas for her work on Fig. 14(b), and to Amalia Thomas and Nathalie Vriend for describing their techniques for casting urethane particles.