Full-field Ultrahigh-speed Quantification of Dynamic Shear Ruptures Using Digital Image Correlation

Producing dynamic ruptures in the laboratory allows us to study fundamental characteristics of interface dynamics. Our laboratory earthquake experimental setup has been successfully used to reproduce a number of dynamic rupture phenomena, including supershear transition, bimaterial effect, and pulse-like rupture propagation. However, previous diagnostics, based on photoelasticity and laser velocimeters, were not able to quantify the full-field behavior of dynamic ruptures and, as a conse-quence, many key rupture features remained obscure. Here we report on our dynamic full-field measurements of displacement, velocities, strains and strain rates associated with the spontaneous propagation of shear ruptures in the laboratory earthquake setup. These measurements are obtained by combining ultrahigh-speed photography with the digital image correlation (DIC) method, enhanced to capture displacement discontinuities. Images of dynamic shear ruptures are taken at 1-2 million frames/s over several sizes of the field of view and analyzed with DIC to produce a sequence of evolving full-field maps. The imaging area size is selected to either capture the rupture features in the far field or to focus on near-field structures, at an enhanced spatial resolution. Simultaneous velocimeter measurements on selected experiments verify the accuracy of the DIC measurements. Owing to the increased ability of our measurements to resolve the characteristic field structures of shear ruptures, we have recently been able to observe rupture dynamics at an unprecedented level of detail, including the formation of pressure and shear shock fronts in viscoelastic materials and the evolution of dynamic friction.


Introduction
Earthquakes occur as dynamic frictional ruptures along preexisting interfaces (or faults) in the Earth's crust. In our laboratory, earthquakes are mimicked by dynamic rupture propagating along the inclined frictional interface separating two pieces of analog material, Homalite-100, pre-stressed in compression and shear [1,2] (Fig. 1). Studies of dynamic ruptures using earlier versions of this setup have addressed a number of important issues in earthquake dynamics, such as confirming the possibility of supershear transition ( Fig. 2(a)), that is rupture propagation at speeds faster than the shear wave, demonstrating the change of rupture mode from crack-like to pulse-like with decreasing fault pre-stress, investigating the importance of the bimaterial effect for the rupture speed and directivity, and studying rupture interaction with the free surface [1][2][3][4][5][6][7][8]. However, these previous studies were unable to provide full-field mapping of displacements, velocities, strains and strain rates during rupture propagation. In those studies, the diagnostics used to capture the rupture behavior were based either on temporally accurate but spatially sparse laser velocimeter measurements or on temporally sparse full-field photoelastic images, which recorded the maximum shear stress field. Our first step towards spatially continuous mapping was to quantify the static fullfield displacements and strains of arrested ruptures [9].
Here we present dynamic measurements capable of resolving the spatial and temporal features of shear ruptures in our laboratory earthquake setup, using ultrahigh-speed photography combined with digital image correlation, enhanced to capture displacement discontinuities on the interface (Fig. 1). This metrological advance allows us to make the leap from qualitative observations of key rupture features, such as that of shear Mach cones in supershear ruptures with photoelasticity ( Fig. 2(a)), to fully quantitative field measurements ( Fig.  2(b)). It also allows us to observe and quantify new features that were undetectable by previously used diagnostics, such as the evolution of dynamic friction coefficient [10] and pressure shock fronts [11]. An example of selected snapshots of displacement, velocity and shear stress full-field maps, obtained Fig. 1 Schematics of the laboratory earthquake experiment. Dynamic shear ruptures spontaneously evolve along the frictional interface, inclined at an angle α, of two Homalite plates under a static prestress load P. Ruptures are initiated by the small burst of a NiCr wire placed across the interface and connected to a capacitor bank. In previous versions of this setup, coherent laser light was transmitted through the stress-birefringent specimen and collected by a gated-intensified high-speed camera to produce a sequence of photoelastic images of the propagating rupture. In the current version, white light produced by a flash light source is reflected by the specimen's surface and captured by a low-noise high-speed camera, typically at 1-2 million frames/s. The portion of the specimen to be imaged, the field of view, is coated by a flat white paint and decorated by a characteristic speckle pattern. The two fields of view discussed in the text are shown in the schematics. Next, the textured images are processed by digital image correlation algorithms to produce a temporal sequence of full-field displacement maps. The displacement fields are then post-processed to produce velocity, strain, strain rate, and stress maps. The displacement sequence is modified from [10]  with the ultrahigh-speed digital image correlation for one of the ruptures studied here, is shown in Fig. 3. Digital image correlation is an optical technique that produces full-field measurements of displacements by analyzing digital images containing a characteristic grey-level content [12,13]. This technique can be used to measure 2D or 3D displacement fields. Here we use its 2D version. In 2D-DIC, surfaces are assumed to be planar and perpendicular to the camera axis and displacements are assumed to be in-plane, to minimize perspective distortions [12,13]. Digital image correlation algorithms determine the displacement field providing the best match between a reference and a deformed image. Correlation algorithms can have either local or global approaches [14,15]. In the local approaches, to regularize the non-uniqueness of the pixel-by-pixel correlation problem, the image matching is performed considering small windows, or Bsubsets^, separated by a distance, referred to as Bstep^, which can be less than half a subset size, i.e. subsets can overlap [12,16,17]. In the global approaches, the pattern matching is performed globally, typically using the finite element method [18,19], which allows enforcing compatibility of the displacement field. In recent years, full-field measurements enabled by DIC have been performed across a wide spectrum of applications (e.g. [12,20]). Most common digital image correlation approaches assume continuous displacement fields and cannot capture discontinuous fields associated with cracks or ruptures. Several techniques have been developed to deal with discontinuous displacement fields [21][22][23][24][25][26][27], however most approaches require applying constraints involving a theoretical interpretation of the experiment, limiting the range of applicability to cases where the theoretical assumptions are valid. In our work, we use the commercial DIC software Vic-2D (Correlated Solutions, Inc.) enhanced to resolve displacements along an interface featuring a displacement jump (see BDigital image correlation of dynamic shear ruptures^).
A significant metrological challenge in quantifying the fullfield behavior of shear ruptures is due to the time-and lengthscales involved in the rupture process. Spontaneously propagating shear ruptures with the speeds in excess of 2 km/s require a minimum temporal sampling of the order of 1-2 million frames/s to capture their temporal evolution, as well as an adequate spatial sampling to resolve their features. In the last decade or so, high-speed digital image correlation applications have increased [28][29][30][31][32][33][34][35][36][37], but their development has been limited by the high-speed camera technologies available.
A key requirement for the digital images to be analyzed with DIC is their low noise level [38][39][40][41][42][43]. This requirement, together with that of adequate spatial and temporal resolution, makes high-speed camera selection not a trivial task. A high-speed camera survey of only a few years ago [39] reveals that Experimentally determined snapshots of (a) interface-parallel displacement, (b) interfaceparallel velocity, and (c) shear stress caused by a supershear rupture (P = 23 MPa, α = 29°). Displacement time-histories are produced by processing the sequence of ultrahigh-speed images with digital image correlation and by employing the symmetryadjustment procedure described in the text. Velocity fields are computed as the displacements time derivatives. The strain fields (not shown here) are obtained by spatial derivatives of the displacement fields. The stress components are computed using linear elastic plane-stress constitutive equations with the dynamic Young's modulus of Homalite to account for its strain-rate dependent behavior cameras with a frame rate above 1 million frames/s, referred to as ultrahigh-speed cameras, can achieve a high spatial resolution (in the megapixel range) but typically have a low record length and rather high electronic noise levels, while cameras that can attain a large number of recorded frames can do so at a low sampling rate. Recent advances in high-speed camera technologies have enabled a higher number of recorded frames in the ultrahigh-speed range [44]. While the spatial resolution of these cameras is still limited, their low noise level makes them good candidates for DIC [32,35]. Ultrahigh-speed photography combined with digital image correlation tailored to capture dynamic ruptures has allowed us to discover new phenomena. We have recently been able to visualize two sets of shock fronts, pressure and shear, emanating at the tip of very fast, spontaneously propagating, supershear ruptures; this previously unexplored phenomenon was found to be related to the viscoelasticity of the material [11]. While shear shock fronts were expected at the tips of supershear ruptures and have indeed been observed before experimentally using photoelasticity (e.g. [1,7,45,46] and Fig. 2(a);), pressure shock fronts were not thought to be possible. Resolving full-field quantities that are sensitive to the dilatation/compression, such as velocity and volumetric strain rate fields, has enabled us to image pressure shock fronts and attribute their formation to the high strain-rate sensitivity of the polymers inducing a spatially heterogeneous material stiffening [11]. Another recent development in the study of dynamic ruptures, enabled by our ultrahigh-speed DIC method, is the ability to track the time evolution of dynamic friction at any point along an interface traversed by a propagating rupture and thus to monitor its dependence on relevant quantities such as slip, slip rate, and normal stress [10]. The measured transient friction response can be used to validate and inform friction laws that are a key input in numerical simulations of earthquake ruptures.
The outline of the present paper is as follows. We start by describing the experimental setup, including the laboratory earthquake experiment, the ultrahigh-speed diagnostics, the digital image correlation approach, and the post-processing procedure to turn the displacement fields into velocity, strain, and strain rate fields. We then discuss characteristics of the full-field structure at the vicinity of a dynamically growing shear rupture. In particular, the question we aim to answer is: what is the best selection of experimental parameters (such as field of view, speckle pattern etc.) and analysis parameters (such as the subset size) to accurately resolve the spatial and temporal features of dynamic ruptures? We anticipate that the selection depends on whether we prioritize the imaging of rupture features in the far field or achieving a high level of accuracy in a small area very near the rupture tip. Finally, we discuss the error analysis that allows us to select the optimal speckle pattern that minimizes the measurement uncertainties, among the infinitely many alternative speckle patterns.

The Laboratory Earthquake Experiment
The laboratory earthquake setup reproduces the main features of a fault in the Earth's crust loaded in compression and shear by the frictionally held interface of two Homalite quadrilateral plates [2] ( Fig. 1). Specimens are obtained by cutting square plates of Homalite-100, with the dimensions 200 mm × 200 mm × 10 mm, using computer-numerical-control (CNC) milling, and producing an interface of inclination angle α (Fig. 1). To obtain consistent surface roughness and repeatability of the dynamic frictional rupture experiments, the following protocol is observed: (i) the mating surfaces of the specimen are polished to a near-optical grade finish; this procedure erases any manufacturing marks due to the CNC cutting; (ii) the surfaces are roughened by employing a micro-bead blasting procedure with abrasive glass media having diameters in the range of 104-211 μm [10,46]; (iii) new test specimens are used in every test. The specimen assembly is mounted in a rig, positioned in a hydraulic servocontrolled loading machine, and compressed with a uniaxial load P (Fig. 1). It is possible to regulate the level of shear τ 0 = P sin α cos α and normal σ 0 = Pcos 2 α pre-stress on the fault by controlling the applied loading P and by setting the inclination angle α. The non-dimensional pre-stress is given by τ 0 /σ 0 = tan α. Dynamic ruptures are nucleated by a local pressure release provided by a rapid expansion of a NiCr wire filament due to an electrical discharge.
One important aspect of our laboratory earthquake setup is that it allows us to produce well-developed dynamic ruptures in samples of tens of centimeters, instead of several meters as would be required for natural rocks [47]. This is enabled by the use of Homalite as analogue material, which has a significantly (~20 times) lower shear modulus compared to rocks. The lower shear modulus significantly decreases all relevant critical length scales, such as the rupture nucleation size needed for dynamic rupture growth [48].

Surface Patterning for Image Matching
In order to provide a characteristic pattern for image matching to the transparent Homalite sample, the specimen's surface is first coated with a uniform layer of Krylon® flat white paint for plastics and then a random black-speckle pattern is applied to the area to be imaged with a dot-on-dot technique (Fig. 4). The speckle patterns shown in Fig. 4 minimize the measurement error, based on a comparative analysis of several alternative patterns (see section BCharacterization of the Measurements Error^and Fig. 18). The target speckle size is 3-6 pixels [12], since this size compromises between resolving sharp displacement gradients, which is achieved with a small speckle size and avoiding aliasing due to under-sampling of the pattern, which occurs when too small speckles are used. In order to be consistently in the range of 3-6 pixels, we adjust the size of the speckles according to the dimension of the field of view. The average feature size is 3-4 and 5-6 pixels for the images of Fig. 4(a) and (b), respectively.

High-Speed Diagnostics
The high-speed setup features an ultra high-speed camera system (Shimadzu HPV-X), capable of recording up to 10 million frames per second, a high-voltage capacitor (Cordin 640) to discharge the NiCr wire and nucleate the rupture, and a highspeed white light source system with two light heads (Cordin 605) (Fig. 1). In selected experiments, we use the welldeveloped technique of laser velocimetry [46,49] to verify the accuracy of the full-field DIC measurements by comparing the velocity time-histories at corresponding locations. The Cordin flash, used in its brightest setting, ramps up to the maximum light intensity within 100 μs and maintains a constant level for about 1 ms. To ensure stable light intensity, the flash is triggered 200 μs ahead of the discharge of the NiCr wire. The HPV-X high-speed camera is triggered in synchronous with the trigger for the wire discharge, or with a delay, depending on the position of the imaging area with respect to the rupture nucleation site. In our experiments, the camera records a sequence of 128 images of the patterns distorted by the propagating rupture with a resolution of 400 × 250 pixels 2 , at temporal sampling ranging from 1 to 2 million frames/s, depending on the experiment, and exposure time of 200 ns.
The HPV-X camera is equipped with prime telephoto lenses, with focal lengths of 105 mm (Nikon Micro-Nikkor 105 mm f/2.8, for larger fields of view) or 200 mm (Nikon AF Micro-Nikkor 200 mm f/4D IF-ED, for smaller fields of view), depending on the field of view (Table 1). A longdistance microscope (model Infinity K2 DistaMax™) is used for the smallest field of view. The camera is mounted on a precision rotation stage (Arca-Swiss C1 Cube), capable of three-axis rotation, which allows us to position the camera axis perpendicular to the specimen's surface and tilt the camera at an angle α to align the pixels of the camera sensor with the specimen's interface. The geared head is in turn mounted on a system of precision translation stages, which allows us to align the camera with the location of the field of view on the specimen. The assembly is fixed on a rail system, which runs in the direction perpendicular to the specimen and enables adjustment of the camera distance from the specimen, depending on lens used and the field of view.

Selection of the High-Speed Camera for Dynamic Digital Image Correlation
The high-speed camera system is the keystone of a dynamic digital image correlation setup. The dynamic events produced in our laboratory setup last 60-120 μs, so to capture the evolution of our dynamic ruptures, the high-speed camera needs to have a frame rate on the order of 1-2 million frames/s (fps), with ideally a large number of recorded frames (e.g. 100 Table 1 Fields of view (FOV) analyzed in this study. The Shimadzu HPV-X camera is equipped with prime lenses of 105 mm or 200 mm focal distance, or a long distance microscope lens (Infinity K2 DistaMax™). The camera distance from the specimen is adjusted depending on the size of the area to be imaged. In the text we discuss results for two FOVs: no.2 (denoted as Blarge^) and no. 5 (denoted as Bsmall^)

FOV (mm 2 )
Pixel size (μm) Lens frames), a small exposure time to avoid blurred images, and an adequate spatial resolution to appropriately sample the field of view of interest. Most importantly, capturing our dynamic events require the ability to at least resolve displacements on the order of a micron, particle velocities of 0.1 m/s, and strains of 100 με. An important requirement to achieve this goal is that the digital images acquired by the high-speed camera need to have a low level of electronic noise.
To find a high-speed camera capable of capturing digital images with the above specifications, we have tested state-ofthe-art camera technology, including gated-intensified, rotating mirror, and on-chip memories cameras, with either chargecoupled device (CCD) or complementary metal-oxide semiconductor (CMOS) sensors. Gated-intensified cameras are based on micro-channel plate (MCP) intensifiers that amplify the light signal of the incoming photons. This allows gatedintensified cameras to take images at very short exposure times (on the order of nanoseconds). A series of images is obtained by splitting the incoming light beam by means of a beam splitter directing the light signal to individual CCD sensors, typically 4 to 8 sensors, depending on the model. As a consequence, gated-intensified cameras offer the highest frame rates (up to 100 billion fps) and the lowest exposure time (0.2 ns). The spatial resolution available with gatedintensified cameras is also among the largest, up to 2000 × 2000 pixels 2 for cameras with a frame rate of 200 million frames/s). However, due to the intensification of the light beams before reaching the CCD sensors, the electronic noise level is comparatively large. Another disadvantage of such systems is that these cameras can record a limited number of frames, which results in a sparse temporal resolution of the entire event. Gated-intensified cameras have been used in the past for digital image correlation applications [38,40,50] with image pre-processing helping to reduce noise and optical distortions [38]. However, the gated-intensified systems we have tested do not allow us to meet the stringent displacement, velocity and strain requirements outlined above, and their use for DIC applications seem to be more appropriate for the highest frame rates were no other technology is available. An example of the displacement fields associated with the dynamic ruptures produced in our setup, and obtained through processing the raw images with Vic-2D, are shown in Fig. 5(a) and (b). These images were acquired with a gated-intensified camera at 1 million frames/s, 50 ns exposure time and MCP gain setting minimized to reduce noise. Though these images provide some, limited, insight into the rupture process they are hard to use in a quantitative sense and they are not amenable to temporal or spatial differentiation to obtain particle velocity or strain fields.
Rotating mirror cameras are high-speed systems based on a rotating mirror, driven by either an electric motor or a gas turbine, splitting the incoming light beam to multiple CCD sensors arranged around a stator. These cameras are also capable of high frame rates and large spatial resolution, for example they can reach 4 million fps with 2000 × 2000 pixels 2 resolution and they keep resolution constant at the highest frame rate. These cameras can acquire a large number of images (on the order of 100) and can attain much lower noise levels compared to gated-intensified cameras. Rotating mirror cameras have been successfully used for dynamic digital image correlation applications [28][29][30]40] and other full-field techniques [51]. However, rotating mirror cameras can pose technical difficulties to operate at the highest frame rates, as these frame rates are achieved with the highest rotating speeds, which cannot be maintained for prolonged period of times. For instance, a rotating speed of 12,500 rev/s may be required to attain 1 million fps; maintaining this rotating speed for a prolonged period of time, such as during calibration or while waiting to trigger the experiment, can significantly reduce the fatigue life of the mechanical components. Further, rotating mirror systems, like gated-intensified cameras, achieve highspeed image acquisition by having multiple CCD sensors (one for each frame). Digital image correlation cannot be performed sequentially on a series of images including frames from different CCDs, as nominally identical images taken by different CCDs are slightly different and correlating them would produce important systematic errors. Typically the dynamic series of images is correlated with a matching reference series taken before dynamic deformation by the same CCDs. In a rotating mirror camera, selected CCDs can be purposely skipped to achieve a desired frame rate for a given mirror rotating speed, so it is not always straightforward to match two image sequences, as two subsequent images in time may not be taken by sequential CCDs.
High-speed camera systems capable of image acquisition and recording on an individual sensor have been developed in the last few years. One such camera is the Shimadzu HPV-X, which has a CMOS sensor featuring multiple onchip memories, and can record up to 128 images at full resolution (400 × 250 pixels 2 ), or 256 at half resolution [52,53]. The HPV-X has been successfully employed in combination with digital image correlation in a number of recent studies [10,11,32,35]. In this study, we employ the HPV-X at full resolution. This technology allows us to correlate each image with the chosen reference image within the same dynamic set of images. Each pixel in the CMOS sensor is also provided with noise reduction circuits. As a result, dynamic images have a low noise level, which enables us to process them digitally with pattern matching algorithms and to produce a temporal evolution of displacements, velocities, strains and stresses. One potential downside of this system is the low resolution (400 × 250 pixels 2 ). However, the actual spatial resolution of a displacement field obtained with DIC depends on a combination of the physical dimension of the subset size (see section BDigital image correlation of dynamic shear ruptures^) and the camera resolution. Displacement fields obtained analyzing with DIC low-resolution digital images, which contain a low noise level, can result in finer spatial resolution compared to displacement fields obtained from high-resolution but noisy digital images. Another potential downside of this system is the large exposure time (200 ns is the shortest), compared to gated-intensified systems (on the order of nanoseconds). Yet, this exposure time is adequate to capture sharp images of our dynamic ruptures without significant blur. An example of displacements fields obtained processing images acquired with the HPV-X is shown in Figs

Digital Image Correlation of Dynamic Shear Ruptures
The digital image correlation (DIC) method is used to analyze the sequence of images acquired with the HPV-X high-speed camera and to produce evolving displacement maps. We employ the correlation software Vic-2D (Correlation Solutions Inc.) enhanced with the BFill boundary^algorithm to capture displacement discontinuities at the interface. The pattern matching analysis is performed between each dynamic image in the sequence and an image of the specimen, belonging to the same sequence, taken before rupture arrival. The components of displacement are computed with respect to the selected reference configuration. Two important parameters in the correlation analysis are the subset size and step size. For any given subset, the 2D-DIC algorithm provides the two in-plane displacement components at the subset center. Using smaller subset sizes produces displacement maps with finer spatial resolution. In turn, smaller subsets result in larger errors as they contain less gray-level information. The choice of the subset size is also influenced by the signal-to-noise ratio (SNR); smaller subset sizes can be used for tests with larger SNR. The step size is the distance between the centers of two nearest subsets. Smaller step sizes increase the density of displacements maps, but also increase the computation time. In this study, we present correlation results obtained using different subset sizes to show their effect on spatial resolution and measurement accuracy of dynamic shear ruptures (see section BEffect of subset size in capturing dynamic rupture features^). An example of displacement fields obtained with subsets of sizes 9 × 9 and 21 × 21 pixels 2 , is given in Fig. 5(c)-(d) and (e)-(f), respectively. As noted earlier, most correlation algorithms are based on the assumption of continuous displacement fields [12] and are not well suited to resolve displacement discontinuities. Among the various strategies proposed to treat crack problems, finite element based formulations featuring nodes with enriched functions along the discontinuity have been proposed [22,23,26,27,54]. These methods have been successful in capturing crack behavior, especially for cases where the crack path is not known Ba priori^. When there is a predefined interface that confines the displacement discontinuity, it is possible to employ two independent meshes on opposite sides of the interface; the mesh size can be substantially smaller than a corresponding subset size, if the fields are regularized with mechanical constraints (e.g. [26]), and this allows to achieve fine spatial resolution. Generally speaking, mechanical regularization schemes produce accurate results for quasistatic problems, for which there are obvious mechanical constraints such as static equilibrium, but they do not have a straightforward extension to highly transient dynamic phenomena. In contrast, the subset-based approach that we employ here does not introduce any mechanical interpretation of the experiment, using only gray-level matching as a criterion, and achieves regularization with the subset size, as illustrated above. However, performing a correlation using a domain containing displacement jumps, such as on the interface of our sample, results in averaging displacements on opposite sides of the interface, as subsets are placed across the interface, and this prevents us from capturing the discontinuities across the interface ( Fig. 5(g) and (h)). In our previous study using DIC for static laboratory earthquake measurements, we used a low-frame-rate CCD camera with high resolution (2048 × 2048 pixels 2 ) and low noise levels [9]. That imaging setup allowed us to use small subset sizes (31 × 31 pixels 2 ) compared to the image size and therefore it was possible to resolve the discontinuous displacement field. In the present case, because of the lower resolution of the high-speed camera (400 × 250 pixels 2 ), typical subset sizes (e.g. 41 × 41 pixels 2 ) are much larger compared to the image size and averaging of the displacement on opposite sides of the interface compromises the proper physical interpretation of the displacement fields (Figs. 5(g) and (h)). Note that simply dividing the field of view in two sub-domains, separated at the specimen interface, would still result in missing the displacement jump as the standard Vic-2D algorithm would only be able to produce the displacement map up to half a subset away from the interface.
In the correlation analyses reported here, we use the BFill boundary^algorithm, developed by Correlated Solutions Inc. with our contribution. The BFill boundary^algorithm enables us to perform correlations over two independent domains, separated at the specimen interface, computing displacement fields up to half a subset away from the interface, and using affine transformation functions to extrapolate the displacements from the center of the subset up to the boundary set by the specimen interface. An example of displacement fields obtained with the BFill boundary^algorithm is provided in Figs. 5(c) to (f) (see also section BThe full-field displacement and velocity structure of dynamic ruptures^). One limitation of this approach is that the displacements on the interface are extrapolated and are not the result of an actual correlation. This approach also produces discontinuities at the correlation-to-extrapolation boundary (half a subset away from the interface, on each side) in the spatial derivatives ∂u i /∂x 2 , with i = 1, 2. The shear stress is one such quantity revealing these artifacts ( Fig. 3(c)). The band around the interface, especially visible at the bottom panel of Fig. 3(c) (t = 48 μs), has a thickness of one subset size, which in this particular case is 41 pixels or 1.9 mm. Future developments may include subsets for which the correlation is performed at edge pixels, rather than at the central pixel. Another limitation is that traction continuity is not ensured at the interface; we are currently developing a technique that would ensure normal and shear stress continuity on the interface (Tal Y., Rubino V., Rosakis A.J., Lapusta N., manuscript in preparation). Despite its limitations, the current method does characterize the fullfield behavior of dynamic ruptures without introducing theoretical interpretations (Fig. 3).
The displacement fields of the dynamic shear ruptures analyzed in this study should have an anti-symmetric and symmetric pattern for the components parallel and perpendicular to the interface, respectively, due to the symmetries of the experimental setup. The fields produced by Vic-2D largely exhibit these properties but with small deviations, likely because the two correlation regions (above and below the interface) are analyzed independently and owing to measurements errors. In order to enforce symmetry/anti-symmetry of the analyzed fields and further reduce the effect of the measurement error on the visualization of the full-field features, here we consider a Bsymmetry-adjustment^procedure where we redefine the displacement field in the direction parallel to the interface u 1 (x 1 , x 2 ) to be anti-symmetric as follows: , for x 2 < 0; and the displacement field in the direction perpendicular to the interface u 2 (x 1 , x 2 ) to be symmetric as: , across the entire field.

Post-Processing of the Displacement Fields
In order to compute velocity and strain fields, we filter highfrequency noise from the displacement fields. Typical local Gaussian filters result in locally averaging displacement information and smoothing out sharp or discontinuous variations. Here, we employ the Non-Local-Means (NL-means) filter [9,[55][56][57] to smooth displacement fields. In contrast with local filters, which smooth each pixel with neighboring pixels regardless of their content, the NL-means filter accounts for the Bcontext^surrounding each pixel. This is achieved by considering windows (neighborhoods) around each pixel and comparing them to neighboring windows. The windows are then averaged with Gaussian weights, where larger weights are assigned to windows that express a higher degree of similarity. This procedure enables efficient image de-noising, preserving sharp features and large gradients. The NL-means filter operates with the following input parameters: the size of the neighborhood N, the search area dimension Ω, which defines the span over which the search of similar neighborhoods is computed, and the noise parameter h, related to the noise level of the signal. In all cases analyzed in this study, we use: N = 3 × 3 pixels, Ω = 21 × 21 pixels, and h = 0.5. We found that a second iteration of the NL-means filter with the same parameters helped to further smooth the displacement fields without loss of information. We have also investigated the effect of filtering parameters both on the displacement and strain fields. The above reported parameters achieve displacement smoothing yet maintaining intact the original signal pattern. An example of filtered displacement sequence is shown in Fig. 3(a). The particle velocity maps are produced by time differentiation of the displacement fields, using the central difference scheme. An example of velocity fields sequence is provided in Fig. 3(b).
Strains are computed from the filtered displacement fields using the finite difference approximation [10], as detailed in Appendix 1. Strain rates are obtained through time differentiation of the strain fields, using the central difference scheme: Stress fields are obtained from strain fields using the standard plane-stress linear elastic constitutive equations, and as detailed in Appendix 1. The symmetry-adjusted velocity, strain, strain rate, and stress fields are computed, using the symmetry-adjusted displacement fields. An example of snapshots of the symmetry-adjusted shear stress fields is given in Fig. 3(c).

Capturing the Full-Field Behavior of Dynamic Ruptures
In this section, we demonstrate how the selection of the field of view and subset size can be adjusted to the level of desired resolution and accuracy and tuned to capture dynamic rupture features (Figs. 6, 7, 8, 9, 10, 11, 12 and 13). We have performed tests with fields of view ranging from 145 × 91 mm 2 to 5.6 × 3.5 mm 2 ( Table 1). All fields of view in Table 1 are imaged employing Nikon prime telephoto lenses (see section BHigh-speed diagnostics^), except for the smallest field of view (FOV no. 7), which is imaged using a long distance microscope (model Infinity K2 DistaMax™). In this study, we discuss the field of view 131 × 82 mm 2 (FOV no. 2 in Table 1), which we refer to as Blarge FOV^in the following, and the field of view 19 × 12 mm 2 (FOV no. 5 in Table 1), which we refer to as Bsmall FOV^. These two fields of view have been chosen as each represents an example of a FOV that either is capable of capturing the rupture features in the far field or allows achieving a higher level of accuracy in a small area in the proximity of the rupture tip. Note that the full-field maps are cropped and their size is slightly smaller than the reported FOV size. Digital images are acquired at 1 and 2 millions frames/s, for the large and small FOV, respectively.

The Full-Field Displacement and Velocity Structure of Dynamic Ruptures
Here we analyze the full-field behavior of two experiments conducted using the large and small FOV and under similar experimental conditions: the far-field applied loads are P = 25.2 MPa and P = 23 MPa for the experiment with the large and small FOV, respectively, while the interface inclination angle is α = 29°in both specimens. These experimental configurations have been chosen as they are known from previous experiments to produce dynamic ruptures with strong signal-to-noise ratio; and therefore such ruptures can be used to evaluate the ability of digital image correlation to capture rupture features. Since the ruptures are triggered using the same nucleation procedure and develop under similar far-field loading conditions, they display similar behavior. Both experimental conditions result in ruptures propagating at supershear speeds, as expected by previous tests performed in analogous conditions and as evidenced by tracking the rupture tip position along the interface. Snapshots of the displacement and velocity fields are given in Figs. 6-10. The two frames shown are taken at 66 and 46 μs after rupture nucleation for the large and small FOVs, respectively. Interface-parallel displacement ( m) ( m) Large field of view Small field of view Interface-parallel displacement for similar ruptures captured with the two fields of view shown in Fig. 1 and Table 1. (a-b) Full-field maps obtained by using a subset size of 41 × 41 pixels 2 and by filtering the resulting fields as discussed in the text. (c-d) Displacement vs. position along the interface, tracked one pixel below the interface, for three different subset sizes (the legend marks the width of the subset in pixels) and unfiltered fields, together with a curve obtained from filtering the field of subset size 41 × 41 pixels 2 . (e-f) Displacement vs. position along path perpendicular to the interface, for the case of 41 × 41 pixels 2 subset size and three processing approaches: no BFill boundary^algorithm, with BFill boundary^but no filtering, and with BFill boundary^and filtering. Note that left the panels have different axes ranges from the right panels. The filtered and unfiltered curves obtained with the largest subset size are on top of each other, indicating that the unfiltered fields are already well smoothed by the subset size. However, small differences not visible in these fields become more important for strain and strain rate fields The displacement fields of Figs. 6 and 7 have been obtained performing correlation analyses using a subset size of 41 × 41 pixels 2 for both image sets. Due to the different pixel sizes, this subset size results in the different physical area sizes: 13.4 × 13.4 mm 2 and 1.9 × 1.9 mm 2 for the large and small FOV, respectively ( Fig. 4 and Table 2). The full-field images of the interface-parallel displacement show that the displacement discontinuities are well captured in both cases ( Fig. 6(a) and (b)). The interface-parallel displacement maps show the bottom half moving in the positive x 1 direction (rightward) and the top half moving in the negative x 1 direction (leftward), consistent with a mode II shear rupture. The interface-parallel displacement, tracked one pixel below the interface, is plotted vs. position in Fig. 6(c) and (d). These curves show the rapid increase of interface-parallel displacement behind the rupture tip, up tõ 180 and 48 μm in the selected snapshots, for the large and small FOV, respectively. To show the difference between the enhanced Vic-2D algorithm and the standard formulation in resolving displacement discontinuities across the interface, we plot the interface-parallel displacement vs. position on a coordinate perpendicular to the interface (Figs. 6(e) and (f) for the larger and smaller fields of view, respectively). The plots indicate that the BFill boundary^algorithm is able to capture the expected sharp displacement discontinuities at the interface, with two rectangular domains separated at the interface, i.e. no subsets across the interface. This is an important development since, without this feature, it would not be possible to study the rupture behavior along the interface. In contrast, the standard correlation with no BFill boundary^algorithm has one domain encompassing the whole field of view and has subsets across the interface that contain information about both halves of the field of view moving in opposite directions. As a result this algorithm averages the displacement on opposite sides of the interface, underestimating the net displacement jump. Note that the interface-parallel displacement fields without and with the symmetry adjustment are quite similar and hence the fields with the symmetry adjustment are not shown.
The interface-normal displacement maps are also generally consistent with shear rupture deformation indicating that the two sides of the interface are moving together in the interfacenormal direction as dynamic sliding is occurring ( Fig. 7(a) and (b)). The interface-normal displacement obtained with the large field of view indicates that the interface moves as much as 16 μm downward behind the rupture tip ( Fig. 7(e)), with the magnitude of downward displacement increasing symmetrically away from the interface, up to about −25 μm (Figs. 7(a) or (c)). The interface-normal displacement map produced with the small field of view highlights another important feature: an upward motion of the interface in front of the rupture tip, ahead of the downward motion. The positive displacement is localized on a small region around the interface with the peak of about 5 μm on the interface (Fig. 7(f)). Note that such downward and upward movements of the interface are expected for shear ruptures due to compressional and dilatational quadrants created by the shear rupture around its tip (e.g., [2]); our full-field imaging procedure allows us to capture and quantify this movement.
Close inspection of the field reveals small interfacenormal displacement discontinuities along the interface (Fig. 7(a) and (b)) in some places, which is not expected for ruptures in our experimental setup captured before the reflections from the specimen boundaries come in (e.g., [5]). This minor apparent discontinuity is likely due to the fact that the two domains are correlated and extrapolated towards the interface independently in the BFill boundary^algorithm, without applying continuity constraints along the interface. (We are currently developing algorithms that would ensure continuity of normal and shear stress along the interface; Tal et al., manuscript in preparation.) To correct this artifact, we employ the Bsymmetry-adjustment^procedure described in the section BDigital image correlation of dynamic shear ruptures^. The interface-normal displacement obtained by using the Bsymmetry-adjustment^procedure is shown in Fig. 7(c) and (d) for the large and small fields of view, respectively. A comparison of the interface-normal displacement along the interface for the symmetryadjusted and standard fields is shown in Fig. 7(g) and (h). The symmetry-adjustment procedure also results in further smoothing the fields. Note that the symmetryadjusted fields result in qualitatively, and often quantitatively, similar conclusions about the along-interface properties compared to the original fields; for example, we checked that the conclusions about dynamic friction in [10] are the same for both fields. Note also that there are cases where normal displacement discontinuities are physically possible and indeed expected in rupturing systems, in which case it would not be appropriate to apply the symmetry adjustment. Indeed, for faults meeting the free surface at an angle, normal stress reduction and even opening is expected (depending on the fault angle and pre-stress level) as the rupture approaches the free-boundary due to the asymmetric rupture interaction with the specimen's surface [8]. Also, in bimaterial systems, the different material behavior on the two sides of the interface may also result in normal stress reduction or even in complete interface opening depending on the levels of both bimaterial contrast and the far field applied pre-stresses [3].
By now, it is clear that, using the large field of view, we cannot resolve all features of the displacement fields obtained using the small field of view. To further illustrate this conclusion, we compare the displacement maps obtained with the small field of view to a zoomed-in version of the large field of view. The close up is taken near the rupture tip, in the rectangular white box marked in Figs. 6(a) and 7(a) so that the size of the zoomed-in field is the same as that of the small field (Fig. 8). The interface-parallel displacement maps of the two FOVs, obtained by using both the standard and symmetryadjustment procedure, are in good agreement ( Fig. 8(a) and (b); Fig. 8(e) and (f)). However, the upward motion of the interface-normal displacement captured by the small FOV is missed by the large FOV (Fig. 8(c) and (d); Fig. 8(g) and (h)), which only shows the downward motion of the interface. Subsequent snapshots of the zoomed-in field reveal hints of positive displacement in the interface-normal displacement, though not as well resolved as in the small FOV. Note that the top panels of Figs. 6 and 7 are plotted using different displacement ranges, while all panels of Fig. 8 use the same displacement range for comparison. In addition, the magnified version of the interface-normal displacement for the large field of view, obtained with the standard correlation approach reveals more marked deviations from symmetry (Fig. 8(c)) compared to the small field of view ( Fig. 8(d)), making it likely that they are the result of correlation and extrapolation errors which are larger for the larger FOV. These deviations from symmetry are corrected by employing the Bsymmetry-adjustment^procedure ( Fig. 8(g) and (h)).
The interface-parallel and interface-normal velocities are shown in Figs. 9 and 10. The interface-parallel velocity is characterized by the formation of two distinct shock fronts, pressure and shear [11]. While shear shock fronts are expected for supershear ruptures and have been observed experimentally before using photoelasticity, e.g. [1,7,45,46], pressure shock fronts were not thought to be attainable until recently [11], as spontaneously propagating ruptures are believed not to be able to travel faster than the pressure wave speed, at least not according to classic theories of rupture developed for isotropic linear-elastic solids [60][61][62][63]. We have recently demonstrated that such features are indeed shock fronts by showing that they satisfy well-known kinematic relationships and that the apparent violation of classic theories is due to local material stiffening associated with the strain-rate-dependent material behavior of polymers [11]. The strain-rate dependence of Homalite-100 [11,58] results in increased wave speeds near the high-strain-rate region of the rupture tip compared to the wave speeds of the material in the far field, which is deforming at lower strain rate. The pressure shock front is formed as the pressure waves emitted at the rupture tip of the spontaneously propagating rupture coalesce in the lowstrain-rate region away from the rupture tip. Because of the formation of the pressure shock front, these ruptures can be regarded as supersonic with respect to the far-field wave speeds of the solid. Tracking the rupture tip along the interface (Fig. 9(e)  Several experimental studies of shear ruptures in polymers [2,7,45,46,59,[64][65][66], including Homalite and PMMA, have reported rupture speeds similar to those reported here and in [11], but their analysis only accounted for uniformly increased elastic constants due to strain rate effects which is not expected to result in the formation of pressure shock fronts. Previous experiments conducted under similar experimental conditions [7,46] as the ruptures presented here may have produced supersonic ruptures, but it was not possible to recognize them as such since photoelasticity, used in previous versions of the setup, is sensitive to the maximum shear stress and cannot capture pressure shock fronts ( Fig. 2(a)). In fact, even when using DIC to plot the maximum shear stress, pressure shock waves are not visible (Fig. 2(b)). Imaging and quantifying the heterogeneous strain rate field (see section BHigh-Speed Measurements of Strains and Strain Rates^and Fig. 13) inducing spatially varying effective elastic properties is a crucial step in explaining the formation of the pressure shock front and could not be achieved with previous diagnostics.
The interface-parallel velocity fields of Fig. 9(a) and (b) also show remarkable similarity with the molecular dynamics simulations of shear cracks [67][68][69][70]. In these numerical simulations, supersonic ruptures are obtained by modeling the material as stiffening with strain. Such material modeling also leads to higher elastic properties and wave speeds near the rupture tip region compared to the far-field elastic properties, and results in ruptures becoming supersonic. However, in our experimental ruptures, the stiffening mechanism is viscoelasticity rather than large-strain hyperelasticity, as the polymers tested are highly strain-rate sensitive but our experiments do not have large strains near the rupture tip [11].
The interface-normal velocity maps captured by the two FOVs show the interface-normal velocity being negative in two symmetrical wedge regions, peaking at about~3 m/s ( Fig. 10(a) and (b) for the standard procedure; Fig. 10(c) for Fig. 7 Interface-normal displacement for similar ruptures captured with the fields of view shown in Fig. 1. (a-b) Full-field displacement maps obtained by using a subset size of 41 × 41 pixels 2 and by filtering the resulting fields. (c-d) Symmetry-adjusted fields obtained using the procedure described in the text. (e-f) Displacement vs. position along the interface for three different subset sizes (the legend marks the width of the subset in pixels) and unfiltered fields, together with a curve obtained from filtering the field of subset size 41 × 41 pixels 2 . (g-h) Comparison between the standard and symmetry-adjusted solution for two different subset sizes. The interface moves up and down, as captured by the small field of view, due to the compressional and dilatational off-interface lobes created by the rupture. The large field of view analyzed with the largest subset size (41 × 41 pixels 2 ) only captures the downward motion, because of the comparatively large size of the subset compared to the localized region of upward motion, as emphasized in Fig. 8. The field produced with the smallest subset size (21 × 21 pixels 2 ) does capture the upward interface motion but its higher noise level affects the computations of strains and strain rates. Note that left the panels have different axes ranges from the right panels the symmetry-adjusted fields). These wedge regions of negative interface-parallel velocity form between the shear and pressure shock fronts. The interface-normal velocity field imaged with the small FOV additionally captures a region of positive velocity near the rupture tip (peaking at 1.5 m/s on the interface), followed by a smaller swing (~0.9 m/s) in the µ µ µ µ negative direction (Fig. 10(e)). This behavior is highlighted by plotting the interface-normal velocity vs. position along the interface (Fig. 10(f)). The large FOV lacks the resolution to capture this finer structure of the velocity field ( Fig. 10(a) and

Zoom-in of large field of view Small field of view (a) (b) (c) (d)
Interface-parallel displacement Interface-parallel displacement

(g) (h)
Interface-normal displacement Interface-normal displacement  Table 1) show an increasingly finer velocity resolution that enables capturing the positive velocity swing. This highlights the need of performing experiments with imaging windows tailored to capture different structures of the rupture fields.

High-Speed Measurements of Strains and Strain Rates
The full-field strain maps are given in Figs. 11 and 12. Since the DIC analysis is performed using as reference configuration the image of the loaded specimen, taken immediately before rupture arrival, the strain maps characterize the strain change with respect to the state of deformation prior to rupture propagation and they do not provide the absolute strain level (see section BCapturing Dynamic Ruptures with the High-speed Digital Image Correlation Method^). The strain change maps have the advantage of isolating the effect of dynamic rupture propagation from the deformations due to static preloading.
The interface-parallel strain field indicates that the portion behind the rupture tip of the upper (lower) half of the imaging window is under a state of tension (compression), as shown in Fig. 11(a) and (b) for the large and small FOV, respectively, in line with the left-lateral propagation of the rupture. The strains obtained from symmetry-adjusted displacement fields are given  Figure 9(a) is modified from [11]. Note that the small deviations from anti-symmetry in the fields of Fig. 9(a) and (b) are corrected by the Bsymmetry-adjustment^procedure ( Fig. 9(b) and (d)) in Fig. 11(c) and (d). It can be observed that there is a qualitative similarity between interface-parallel strain and interface parallelvelocity ( Fig. 9). They both have an anti-symmetric field but with inverted polarity. This is consistent with the relation u : 1 ¼ −V r ε 11 , valid under the assumption of rupture propagation at steady state, where u : 1 ¼ ∂u 1 =∂t is the interface-parallel velocity, ε 11 = ∂u 1 /∂x 1 the interface-parallel strain, and V r is the rupture speed. The interface-parallel strain map also displays the pressure and shear shock fronts ( Fig. 11(a) and (b)). Tracking the interface-parallel strain along the interface shows a peak in strain immediately behind the rupture tip and then a drop to a lower strain level (Fig. 11(e) and (f)), while tracking it perpendicularly to the interface shows the discontinuity across the interface line ( Fig. 11(g) and (h)). The strain behavior captured by the two fields of view is qualitatively the same, but the higher accuracy of the small FOV allows us to better interpret strain variations in the near field. In particular, the peak strain levels measured by the large FOV are lower compared to those measured by the small FOV, against the expectation of them being the same or larger, due to the slightly higher level of static preload. This is probably the effect of averaging displacements over the comparatively larger physical dimension of the subset size.
The full-field map of the interface-normal strain shown in Fig. 12(a) and (b) (Fig. 12(c) and (d)) indicates that the region behind the rupture tip and immediately above (below) the interface is in compression (tension). The large FOV also and velocity (Fig. 10(a)). The negative (downward) interfacenormal displacement loads in tension the strain wedge above the interface, and in compression the strain wedge below the interface. At the same time, the interface-normal displacement being negative in the off-interface region and near zero around the interface and behind the rupture tip leads to a compressive (tensile) region above (below) the interface. The small FOV displays only some traces of the antisymmetric wedge features ( Fig. 12(b) and (d)), since they fully develop outside the imaged area, so the small FOV is not adequate for capturing the complete structure of the interface-normal strain. The full-field map of the interface-parallel strain rate displays an increase (decrease) in strain rate at the rupture tip above (below) the interface, followed by a decrease (increase). This pattern is consistent with the region above (below) the interface undergoing a state of tension (compression) as the rupture propagates. Behind the peak in the interface-parallel strain, levels of strains are lower (Fig. 11) and this leads to the strain rate changing sign. Towards the back of the rupture strains are nearly constant which result in nearly zero strain rates. This structure can also be inspected by tracking the interface-parallel strain rate along (Figs. 13(e) and 14(f)) and perpendicularly ( Fig. 13(g) and (h)) to the interface. The interface-parallel strain imaged with the large FOV displays the two sets of pressure and shear shock fronts ( Fig. 13(a)) but the pressure shock front is more prominent since ε 11 is sensitive to compression/dilation. The faint trace of the shear shock front is due to some coupling in the deformations. While the large FOV allows us to visualize the formation of these shock fronts, the small field of view ( Fig. 13(b)) allows us to better resolve the level of strain rate. In fact the peak of ε 11 along the interface obtained with the small FOV is~4.5 • 10 3 s −1 (Fig. 13(f)), while by the large FOV is only~0.7•10 3 s −1 (Fig. 13(e)), due to the spatial smoothing deriving from a larger physical size of the subset.

High-Speed Measurements of Stress and Friction
The total shear stress is computed by adding the (nonuniform) shear stress change inferred from DIC to the resolved (uniform) shear stress level computed from the applied far-field load P and inclination angle α, as detailed in BPost-processing of displacement fields^and Appendix 1. Selected snapshots of the total shear stress are shown in Fig. 3(c) for the case of small FOV, and the maximum shear stress is shown in Fig. 2(b), for the same case. The stress fields shown in these figures are computed from the symmetry-adjusted displacements. The shear stress map indicates an increase in stress level ahead of the rupture tip up to a peak level and then a drop to a dynamic, residual level. The ratio of shear to normal stress gives the friction coefficient. Using this technique, we can compute the friction evolution for spontaneously propagating ruptures, at any point along the rupturing interface, and to study its dependence on any other variables, such as slip (i.e. relative displacement across the interface) or its rate, slip velocity [10]. Slip and slip velocity functions are obtained by tracking, along the interface, the interfaceparallel displacement (Fig. 3(a)) and velocity (Fig. 3(b)), respectively. For example, we have recently been able to observe the initial (transient) frictional strengthening with rapidly increasing slip rate, called the Bdirect effect,ŵ hich results as the rupture tip traverses the point of observation in the interface and creates a sudden increase in interfacial sliding from almost zero sliding speed to more than ten meters per second. This transient increase of frictional resistance associated with a positive jump in sliding speed (slip velocity) is predicted by rate-and-state friction formulations available in the literature [10,71,72]. As the rupture tip completely traverses the observation point, steady-state sliding conditions are eventually established in its wake. During that process, rapid weakening of friction with slip velocity is observed, down to a residual level, due to the shear stress drop and constant normal stress level at the interface. This behavior has been observed for a range of experimental conditions and has revealed enhanced dynamic weakening with slip velocity at steady state, consistent with the activation of flash heating [10,[73][74][75][76]. The dependence of friction on slip velocity, combined with other effects, has been already shown by several laboratory experiments [64,72,75,[77][78][79]. However, typical friction experiments impose slipvelocity histories and measure the resulting averaged friction resistance, assuming uniform sliding at the interface under consideration. In contrast, our dynamic measurements, enabled by ultrahigh-speed DIC, allow tracking of local friction evolution for individual spontaneously propagating ruptures and do not depend on the assumption of uniform sliding at the interface [10].

Effect of Subset Size in Capturing Dynamic Rupture Features
In order to explore the role of subset size in resolving the fullfield maps and reducing noise, we perform correlation analyses employing a range of subset sizes ( Table 2). The interfaceparallel displacement is tracked along the lower side of the interface for the larger and smaller fields of view in Fig. 6(c) and (d), respectively. Likewise, Fig. 7(c) and (d) trace the interfacenormal displacement profile along the interface. The displacement vs. position curves shown are obtained from correlation analyses performed with three subset sizes: 21 × 21 pixels 2 , 31 × 31 pixels 2 , and 41 × 41 pixels 2 . The displacement profiles obtained from correlations performed with smaller subset sizes (21 × 21 pixels 2 and 31 × 31 pixels 2 ) are able to capture the rupture behavior, and they follow closely the displacement profiles obtained with the larger subset size (41 × 41 pixels 2 ). However, as expected the correlations with smaller subset sizes display significantly noisier curves compared to the larger subset size. Containing the noise level is especially important when differentiating the displacement fields in time and space to obtain velocities, strains and strain rates (see Figs. [8][9][10][11][12][13]. Hence, we chose the larger subset size (41 × 41 pixels 2 ). To smoothen the displacement fields and reduce noise levels further, we filter displacement components with the NL-means filter (see section BDeveloping the High-speed Digital Image Correlation Method in the Laboratory Earthquake Setup^). The effect of filtering is noticeable in the displacement fields (Figs. 6 and 7), but it becomes more apparent in the velocity (Figs. 9 and 10) and particularly in the strain (Figs. 11 and 12) and strain rate fields (Fig. 13). An important question to ask when selecting the subset size is whether a given subset size is capable of resolving all relevant physical features of interest, while providing at the same time fields smooth enough to be amenable for velocity and strain computation. The interface-parallel velocity vs. position along the interface plot ( Fig. 9(c)) shows the velocity rapidly increasing behind the rupture tip, reaching a peak, and then decaying to near-constant level for all correlations shown in the plot, with the correlation performed using a subset size of 21 × 21 pixels 2 showing larger oscillations. One peculiar feature of the interfaceparallel velocity curve corresponding to a subset size of 21 × 21 pixels 2 is a smaller peak ahead of the main peak, not captured by the larger subset sizes, which may at first appear to be the result of noise in the correlation analysis. Previous observations of our laboratory raptures indicate that the Bdouble peak^is actually a physical feature of the rupture due to the reflection of the shear shock front on the rear surface of the specimen [10,46]. To confirm that this is a physical feature and not an artifact due to noise, we compare the interface-parallel velocities obtained by DIC measurements to simultaneous measurements at the same locations performed with laser velocimeters. In a recent study, we presented a comparison of DIC with velocimeters measurements for an intermediate imaging area size (50 × 31 mm 2 , FOV no. 4 of Table 1), which indicated excellent agreement between the two measurements [10]. However, it is not clear whether the large physical dimension (Table 2) of the subset size (13.4 × 13.4 mm 2 ) employed in the large FOV (131 × 82 mm 2 , FOV no. 2 of Table 1) is adequate to resolve important rupture features. To address this point, we show a comparison between velocimeter measurements and velocity time histories (for a point at the center of the FOV), obtained form DIC analysis of the large FOV and employing two subset sizes, at the opposite ends of the size range analyzed to emphasize the differences. The velocimeter measurements are performed at two locations, immediately above and below the interface, respectively. Here we report the results for the velocimeter measurement just below the interface, as the one above is symmetric and leads to the same conclusions. The velocimeter trace clearly shows the presence of two peaks in the particle velocity time history (Fig. 14(a)). The interfaceparallel velocity time history obtained from DIC analysis using the subset 9 × 9 pixels 2 well captures the double peak, but it exhibits a large noise level with oscillations of amplitude on the same order of magnitude as the two peaks ( Fig. 14(a)) making it difficult to differentiate noise oscillations from real signal, in the absence of an independent velocimeter measurement. The smallest subset size that enables us to clearly tell apart the two peaks from noise oscillations is 15 × 15 pixels 2 (not shown). Yet, the noise level is still important and the strain and strain rate fields obtained with this subset size are severely affected, making their interpretation difficult. Increasing the subset size reduces noise level but compromises spatial resolution. The velocity time history obtained using a subset of 41 × 41 pixels 2 has a significantly lower noise level, but it fails to capture the double peak ( Fig. 14(a)).
To illustrate this point further, we compare the particle velocity measured by the velocimeter at the two peaks (marked as A and B in Fig. 14(a)) and at the trough in between (marked as C in Fig. 14(a)) with corresponding DIC measurements as a function of the subset size (Fig. 14(c)). The plot shows that the double peak is resolved for subset sizes smaller than 21 × 21 pixels 2 (6.9 × 6.9 mm 2 ). The three velocity peaks measured by DIC are also reported for the case of small FOV in Fig. 14(d). It is not possible to perform the comparison with velocimeters in this case as the retro-reflective tape needed as target for the laser velocimeters would be too large with respect to the field of view and would compromise the correlation. Nonetheless, we observe that the velocity time histories obtained using subsets of different sizes for the small FOV all lead to recover the double peak in particle velocity (Fig. 14(b) and (d)). The correlation performed with the smallest subset size (9 × 9 pixels 2 ) displays some discrepancy compared to the larger subset sizes, due to higher noise level occurring for this small subset. We can use the velocity spatial profile on the interface of this well-resolved FOV to determine the spacing between peaks and get an insight on the maximum subset size to be used to resolve these features. The peak-to-peak distance is 14 mm and the peak-to-trough distance is~7 mm ( Fig. 14(b)). We observe that subset sizes whose physical dimension is larger than the smallest relevant length scale (~7 mm) fail to resolve the double peak in the velocity field (Fig. 14(c)). Only subsets smaller than 21 × 21 pixels 2 (6.9 × 6.9 mm 2 ) in the large FOV can resolve the double peak structure. This is consistent with the fact that it is not possible to resolve sinusoidal displacements with a period smaller than the subset size [80]. On the other hand, all subset size considered for the small FOV can resolve it, since they are all smaller than the spacing between the characteristic length scale ( Table 2).
In summary, using a subset size smaller than the size of the physical feature to capture ensures appropriate resolution of the feature, but selecting larger subset sizes, at the cost of filtering out small-length-scale features, allows us to obtain smoother fields and to highlight large-lengthscale features that otherwise would be obscured by noise. For example, the plots of the interface-parallel strain vs. position (Fig. 11(e)-(h)) reveal that the lower subset size choice of 21 × 21 pixels 2 is not able to capture the strain behavior as the physical interpretation of the curves is compromised by noise, while the subset size of 41 × 41 pixels 2 effectively filters out noise for both cases of field of view shown. The same conclusion is true for all other components of strain and strain rates (Figs. 11-13). This

Characterization of the Measurement Error
The measurement accuracy is affected by a number of factors including the electronic noise of the camera, lighting variations during the observation time window, suboptimal speckle pattern, subset size employed for the correlation analysis, and the correlation algorithm [81-84, 80, 85]. There are several ways to quantify the measurement error and to single out separate error sources. For example, an approach to evaluate correlation algorithms consists in analyzing images deformed experimentally or numerically by a known (uniform or otherwise) displacement field and quantifying the difference between the known fields and the correlation results [86]. Using this approach it is possible to determine the spatial resolution that can be achieved by a given correlation algorithm for a given subset size by finding the minimum length scale that can be resolved of a known displacement field of variable spatial frequency [86]. In this section, we assess the measurement error by correlating nominally identical images and using the resulting displacement fields as a quantification of the uncertainties. While this error quantification cannot be used to assess the goodness of the correlation algorithm to capture the highly dynamic transients associated to the propagation of our dynamic ruptures it can be used in a comparative sense, to assess the uncertainties due to the speckle pattern employed, high-speed camera noise, environmental factors (such as lighting variations through the image sequence), and analysis parameters (such as subset and step size, and filtering parameters).

Selection of Speckle Pattern
We produced six different types of speckle arrangement to assess the pattern that minimizes the measurement error. All patterns are obtained by coating the specimen surface with a uniform layer of Krylon® flat white paint for plastics and then applying a random black-speckle pattern. The patterns differ by the average speckle size, size distribution and density of the speckle features. The black speckles are deposited by either a dot-on-dot application, where the size of the feature is set by the pen point size, or by spraying a fine mist of black paint. We consider four patterns whose speckles are obtained with the dot-on-dot technique and two patterns obtained with the spray paint technique. For the dot-on-dot patterns, the speckles can be either small (sampled by a patch of approximately 3 × 3 pixels 2 ) or large (sampled by a patch of 6 × 6 pixels 2 ); and the feature distribution can be either dense or sparse, based on whether the distance between features is the same or larger than the average feature diameter. The four combinations of patterns are denoted in Fig. 15 by their feature size and density as: large dense ( Fig. 15(a)), small dense ( Fig. 15(b)), large sparse (Fig. 15(c)), small sparse ( Fig. 15(d)). The two patterns obtained with the spray paint technique are produced by either a coarse or fine mist of paint. The two patterns are denoted as: spray coarse (Fig. 15(e)) and spray fine (Fig. 15(f)). Since it is difficult to fine-tune the speckles size when using the spray paint technique, even with an airbrush, the feature size distribution is much broader with this speckling technique than in the case of the dot-on-dot technique.
In comparing different types of speckle patterns we employ the Shimadzu HPV-2 high-speed camera, which is a similar system to the Shimadzu HPV-X used in all other tests of this study. Since both high-speed camera sensors are based on multiple on-chip memories, we expect conclusions found by comparing patterns with the HPV-2 to hold with the HPV-X though the actual error levels would be different. For this reason, the error measurements presented in this section are normalized with respect to the reference values of the pattern that produces the lowest error levels.
Two nominally identical images of each type of speckle pattern are taken with the HPV-2 camera. For each pair of images we perform a correlation analysis with a subset size of 21 × 21 pixels 2 . We use the displacement fields obtained from these correlations as an estimate of the measurement error produced by each pattern. We conduct a statistical analysis over the displacement fields to determine the bias and standard deviation for each pattern. The speckle pattern that provides the smallest values of bias and standard deviation over those analyzed is the BSmall dense^pattern ( Fig. 15(b)), produced with the dot-on-dot technique. Note that the standard deviation levels are typically an order of magnitude larger than the bias (e.g. see Fig. 16 for the speckle pattern of Fig. 4(a)), and therefore are more important to characterize the error. The actual values of bias and standard deviation are not given here, for each case of speckle pattern of Fig. 15. Instead, to facilitate the comparison between patterns, Table 3 presents the error estimates normalized with respect to the BSmall dense^pattern. The patterns BLarge dense^and BSpray fine^also provide relatively low error levels. One shortcoming of the BLarge dense^pattern is actually not linked to the speckle size but to the clustering of several particles with consequent loss of gray scale gradient over a length scale comparable to that of the subset size. This results in a standard deviation 20% larger compared to that of the BSmall dense^pattern. Declustering of features improves the performance of this pattern and makes its error level comparable to the small dense pattern. In fact, the error levels of the images in Fig.  2(a) and (b) are very similar. The BSpray fine^pattern has the advantage of being relatively easy and fast to produce but producing a pattern by spray paint has limited capability of controlling particle size and density, resulting in a higher error variance (30-40% larger compared to the BSmall dense^pattern). Since the Shimadzu cameras have multiple on-chip memories, a smaller part of each pixel is photosensitive. The photosensitive portion of the sensor is referred to as fill factor. The CMOS sensor of the HPV-X camera has 128 on-chip memories/pixel and a fill factor of 37% [53]. By comparison, conventional CCD sensors have rather higher fill factors, typically larger than 90% [41]. As a consequence, features that may be detected by sensors with a higher fill factor may not be captured by this sensor. This is why to avoid spatial aliasing it is important to make sure that speckles are well sampled, by at least 3-4 pixels across. We prefer the dot-on-dot technique as it gives the smallest error and allows us to control the speckle size within a very narrow range as opposed to the spray paint, which can easily result in features smaller than 3-4 pixel and produce spatial aliasing. The remaining three patterns perform poorly because they do not contain dense enough information within the length scale of the subset size chosen. 21

Quantification of the Measurement Accuracy for Displacement, Strain and Stress Fields
Since the BSmall dense^pattern provides the smallest error, we produce speckle patterns on our specimens based on the same characteristics as this pattern ( Fig. 4(a)) or a de-clustered version of the large dense arrangement ( Fig.  4(b)). Both patterns of Fig. 4 have well separated features (no clustering) and a speckle size sampled by patches 4 -8 pixels in diameter. To characterize the error associated with the measurements presented in this paper, we take a sequence of nominally identical images of both the speckle patterns of Fig. 4 with the Shimadzu HPV-X camera, the same camera used for the dynamic sequence, and we use the same optics and lighting settings as for a dynamic test. The sequence comprises 128 images taken at 1 million frame/s and with 200 ns of exposure time. Images are correlated taking the first frame in the set as reference and using the subset sizes of Table 2.
The correlation analysis produces a sequence of interfaceparallel and interface-normal displacement fields whose deviation from zero displacement provides an error estimate. Performing the statistical analysis provides a time series of standard deviations and means (or biases) for each displacement component (Fig. 16). To obtain one combined value of bias and standard deviation for the whole sequence and for each component of displacement, we compute the bias of the sequence as a mean of the means, and the combined standard deviation as the squared root of the combined (or overall) variance (Appendix 2).
To show the magnification effect on the measurement error, we compare the overall probability distribution functions obtained for the two fields of view studied above. When the displacement uncertainty is measured in pixels the error characteristics of the two fields of view are similar, indicating that the speckle pattern used, as well as other experimental conditions, are similar for the two cases. The probability distribution functions of the interface-parallel displacement (measured in pixels)  characterizing the full sequence of images obtained for the large and small FOV are shown in Fig. 17(a), for the unfiltered correlation analyses performed with a subset size of 41 × 41 pixels 2 .
The probability distribution function of the interface-parallel component of displacement (not shown in the plot) is similar. While measuring the error in pixels is a useful non-dimensional assessment of the effect played by the several factors influencing the measurement, it is important to quantify the error in the relevant physical length scale to assess the accuracy of our dynamic ruptures measurements. The pixel error measurements are converted to micrometers multiplying them by the pixel size, which is 327.9 μm/pixel for the large and 46.5 μm/pixel for the small fields of view, respectively ( Table 1). As a result, the displacement error for the large FOV is approximately 7 times larger compared to the case of small FOV (Fig. 17(b)), when measured in microns. The subset size plays a key role in controlling the measurement error. To quantify its effect, as well as the role of magnification, we plot the overall standard deviation vs. subset size for both displacement components and fields of view, with the measurement error expressed in pixels and microns in Fig. 17(c) and (d), respectively. The plots show the standard deviation rapidly decreasing with increasing subset size. The standard deviation expressed in pixels is larger for the case of small FOV, as its average speckle size is larger and subset sizes contain less gray level variations. However, this difference only becomes important for subset sizes smaller than 21 × 21 pixels 2 . For example, using a subset size of 9 × 9 pixels 2 , the interfacenormal standard deviation of the small FOV is 22% larger than the large FOV, while for a subset size of 41 × 41 pixels 2 , it is only 4% larger. When the measurement error is expressed in microns, the standard deviation is much larger for the case of large FOV, due to the magnification effect. For instance, the interface-normal standard deviation of large FOV is larger by a factor of 5.7, using a subset size of 9 × 9 pixels 2 , and it is 6.7 times larger for a subset size of 41 × 41 pixels 2 . Another way of expressing the effect of the magnification factor is to track the interface-normal displacement error along the interface for both large and small FOV and using correlations with three different subset sizes for comparison (Fig. 18). As expected, the displacement error measured in pixel is similar for the two FOV ( Fig. 18(a) and (b)), while it is significantly different when it is expressed in microns ( Fig. 18(c) and (d)).
In order to quantify the strain measurement uncertainties, we use the displacements fields obtained from nominally identical images to compute strain fields (see section BPost-   processing of the displacement fields^) and we take the deviations from zero of these strain fields as an estimate of the strain error. We analyze the sequence of strain errors statistically analogously to the treatment of displacement errors. As expected, increasing the subset size reduces the measurement error, and filtering the displacement fields before computing strains results in even smaller errors (Fig. 19). For example, the pooled standard deviation of the interface-parallel strain is 152.0 microstrains for a subset of 31 × 31 pixels 2 , and drops to 95.2 microstrains for a subset of 41 × 41 pixels 2 . Filtering the displacement fields before computing strains lowers the standard deviation to 18.6 microstrains.

Conclusions
Quantitative full-field visualization of highly transient phenomena such as dynamic shear ruptures has required us to tackle a threefold challenge. First, the stringent specifications on temporal and spatial resolution of displacements, velocities, strains and strain rates make the high-speed camera choice not a trivial task. Only recently the advances in highspeed camera technology have made it possible to record a large number of images at a high enough frame rate with a low enough noise level, making them suitable for ultrahigh-speed digital image correlation applications. Yet, these improvements come at the cost of lower spatial resolution and comparatively larger exposure times. We have concluded, after extensive testing of state-of-the-art technology, that the benefit of lower noise level outweighs the drawback of lower resolution when it comes to analyzing images with correlation algorithms, where noise minimization is paramount. Second, resolving the interfacial discontinuities, associated with rupture propagation, requires appropriate correlation algorithms and post-processing approaches. We have used an enhanced version of the commercial code Vic-2D, the development of which was motivated by our requirements and which is capable of resolving displacements discontinuities on the interface. This algorithm enables imaging many discontinuity features, as shown in this work. However, this approach is based on extrapolating displacements computed, by correlation, half a subset size away from the interface and does not guarantee continuity of tractions along the interface. Enforcing the expected symmetries in displacements for the specific experiments studied in this work removes some of the associated artifacts and highlights features that do not depend on the associated error. We are currently working on developing alternative approaches that would satisfy interfacial conditions without introducing any other theoretical interpretation in the measurements. Third, the complex rupture behavior results in the formation of features spanning a range of length scales, which cannot be captured by a unique set of experimental and analysis parameters. By employing both large and small fields of view, we have designed experiments tailored to visualize either larger-scale or localized features, and we have chosen analysis parameters in order to compromise between fine resolution of sharp rupture features and noise mitigation through spatial averaging. In summary, due to the complex interplay of all these factors, resolving the spatiotemporal features of dynamic ruptures is only possible by the meticulous selection of multiple experimental and analysis parameters. This study highlights the wealth of information that is possible to obtain about the behavior of shear ruptures by ultrahigh-speed photography combined with digital image correlation. The level of detail at which we can now study shear ruptures was, until recently, only attainable with numerical simulations. Using ultrahigh-speed DIC, we have recently been able to visualize new phenomena, such as the formation of pressure and shear shock fronts associated with supersonic ruptures [11], and to more accurately measure quantities that were only possible to obtain in an average sense, such as evolving dynamic friction histories inferred by tracking individual, transiently growing, rupture events without invoking the assumption of uniform interfacial sliding [10]. Indeed, using ultrahigh-speed DIC tailored to study dynamic rupture has significantly expanded the horizon of observable quantities and phenomena. The application of this method can enable the study of the full-field behavior of many other key rupture characteristics, including crack-like vs. pulse-like behavior, attenuation of the particle motion away from the interface, and rupture interaction with the free surface, among other highly transient dynamic phenomena. where u 1 (i, j, k) and u 2 (i, j, k) are the interface-parallel and interface-normal displacement components, respectively, for pixel (i, j) and frame k, expressed in microns; 2h s defines the stencil size and s is the step size, both expressed in pixels; p is the pixel size, expressed in microns. Here, we take h s = 1 pixel. Note that the minus sign in ε 22 (i, j, k) and in the first term of ε 12 (i, j, k), is due to the fact the that the row index i increases in the opposite direction of the x 2 axis. Close to the interface, we use the backward or forward difference scheme to compute strains above and below the interface, respectively. Below the interface, the forward difference approximation reads: The stress fields are computed from the strain fields employing the standard plane-stress linear elastic constitutive equations. Since Homalite is a strain-rate dependent material, we use the dynamic Young's modulus E d = 5.3 GPa [58], corresponding to high-strain rates, to compute the dynamic stress change [10,59], together with a Poisson's ratio of ν = 0.35. The displacement fields are computed using the loaded specimen configuration (before rupture) as reference, so the strains and stresses computed from these fields are changes over the reference configuration. To recover the actual level of stress, we add the computed levels of normal and shear pre-stresses (σ pre 22 x 1 ; x 2 ð Þ¼σ 0 and σ pre 12 x 1 ; x 2 ð Þ¼τ 0 , respectively) to the DIC measured stresses as: where σ 22 (x 1 , x 2 , t) and σ 12 (x 1 , x 2 , t) are the total levels of normal and shear stresses, respectively, and Δσ DIC 22 x 1 ; x 2 ; t ð Þand Δσ DIC 12 x 1 ; x 2 ; t ð Þare the normal and shear stress changes obtained from DIC. The total normal stresses, as well as the prestresses, are positive in compression, following the geomechanics sign convention, while the stress changes obtained from DIC are expressed using the solid mechanics sign convention; this is why the second term at the right hand side of Eq. (7) is negative. For left lateral ruptures, such as those analyzed in this study, the shear stress on the positive face of a material element (e.g. with normal x 1 ) is positive when acting in the negative axis direction (of the x 2 axis); this explains the minus sign in Eq. (8), since this sign convention is opposite to the one adopted in solid mechanics. For a right lateral rupture, the shear stress on the positive face of a material element is positive when acting in the positive axis direction and therefore the second term at the right hand side of Eq. (8) should have a positive sign. This convention allows expressing stresses most commonly encountered as positive quantities. For example, the constant level of shear stress τ 0 = 9.8 MPa before rupture arrival in Fig. 3(a) is positive in this sign convention, and the stress change due to the rupture propagation results in a shear stress drop, in analogy to the study of natural earthquake ruptures. The procedure of adding the computed levels of pre-stress to the measured ones to obtain the total level of stress is justified by the fact that the resolved levels of shear and normal pre-stress are nearly uniform along the interface [10]. The uniformity of pre-stresses in our experiments is supported by earlier studies providing direct evidence using photoelastic images of preloaded specimens and indirect evidence through repeatability of ruptures in different experiments performed under the same far-field experimental loading [5,6,46]. The uniform distribution of normal and shear pre-stresses is also consistent with the nearsteady rupture propagation through our observation window.
are p = f − 1 independent groups of measurements, where f = 128 is the number of frames, with the size of each group being q k = m x n = 250 x 400 pixels 2 , and the total number of measurements being N ¼ ∑ p k¼1 q k . Since q k = q is the same for all frames, the right-hand side of Eq. (9) can be recast to express the combined variance s 2 c as: The sequence of u k and s k vs. frame number is shown in Fig. 16(a), for the interface-parallel displacement component captured with the large FOVand analyzed with a subset size of 41 × 41 pixels 2 . These values are used to produce Gaussian functions g k u ð Þ ¼ 1= s k ffiffiffiffiffi ffi 2π p À Á •exp − u− u k ð Þ 2 = 2s 2 k À Á for each frame, while the overall mean and standard deviation produce a probability distribution function for the entire data set ( Fig. 16(b)).
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.