<figure><img src="image_1.png" /></figure>

HAL
open science

# How cracks are hot and cool: a burning issue for paper

Renaud Toussaint, Olivier Lengliné, Stéphane Santucci, Tom Vincent-Dospital, Muriel Naert-Guillot, Knut Jørgen Måløy

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Renaud Toussaint, Olivier Lengliné, Stéphane Santucci, Tom Vincent-Dospital, Muriel Naert-Guillot, et al.. How cracks are hot and cool: a burning issue for paper. Soft Matter, 2016, 12, pp.5563-5571. 10.1039/C6SM00615A. hal-01340886

HAL Id: hal-01340886

https://hal.science/hal-01340886

Submitted on 5 Jul 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
How cracks are hot and cool: a burning issue for paper†

Renaud Toussaint,\*a,b Olivier Lengliné,\*a,b Stéphane Santucci,\*c,b Tom Vincent-Dospital,\*a Muriel Naert-Guillot,\*a and Knut Jørgen Måløy\*d,b

Abstract

Material failure is accompanied by important heat exchange, with extremely high temperature – thousands of degrees – reached at crack tips. Such temperature may subsequently alter the mechanical properties of stressed solids, and finally facilitate their rupture. Thermal runaway weakening processes could indeed explain stick-slip motions and even be responsible for deep earthquakes. Therefore, to better understand catastrophic rupture events, it appears crucial to establish an accurate energy budget of fracture propagation from a clear measure of the various energy dissipation sources. In this work, combining analytical calculations and numerical simulations, we directly relate the temperature field around a moving crack tip to the part $\alpha$ of mechanical energy converted into heat. Monitoring the slow crack growth in paper sheets with an infrared camera, we measure a significant fraction $\alpha = 12\% \pm 4\%$. Besides, we show that (self-generated) heat accumulation could weaken our samples with microfibers combustion, and lead to a fast crack/dynamic failure/ regime.

<footnote>a Institut de Physique du Globe de Strasbourg, CNRS, EOST-University of Strasbourg, 5 rue Descartes, 67084 Strasbourg Cedex, France. Fax: +33 3 68 85 01 25; Tel: +33 3 68 85 03 37; E-mail: renaud.toussaint@unistra.fr</footnote>

<footnote>b Centre for Advanced Study at The Norwegian Academy of Science and Letters, Drammensveien 78, 0271 N-Oslo, Norway.</footnote>

<footnote>c University of Lyon, ENS de Lyon, University Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France.</footnote>

<footnote>d Department of Physics, University of Oslo, PB 1048 Blindern, NO-0316 Oslo, Norway.</footnote>
Computing the change in energy in a linear elastic solid due to the creation of a fracture surface, Griffith proposed in his pioneering work [1] a criterion for the progress of fracture in brittle solids. Irwin [2] completed this description taking into account the presence of a plastic process zone around the crack tip, where energy is dissipated in processes which are not pure surface energy. No thermal effects were considered in these theories. Later on, it was recognized that local temperature elevation could play a role in fracture propagation, and Rice and Levy [3] were the first to derive the temperature field around a crack tip. Assuming that all work done was converted into Joule heating – the other extreme case compared to the previous assumptions, they derived a criterion for crack propagation in glass and PMMA based on the local temperature rise. Experimentally, fast rupture in PMMA was shown to lead indeed to high temperatures around the crack tip [4]. Recently, the spectral characteristics of photonic emissions observed during fast rupture in glass allowed to determine a local temperature elevation of several 1000°C around the emitting crack tip [5]. Similar temperature rises were observed during shear fracture in PMMA [6–8].

Moreover, temperature has also been shown to determine the slow dynamics of rupture, controlling the lifetime of various materials submitted to a constant load (creep tests) or to rising load via an Arrhenius law [9–14]. Since then, thermodynamics has slowly emerged as a framework to describe slow fracturing processes. In the case of brittle materials, several statistical models have been proposed to predict the lifetime [15–18], the average dynamics [11], as well as the burst size distribution of a slowly growing crack [19]. In these models [12, 13], the thermal noise inside the material induces stress fluctuations locally at the crack tip, triggering micro-cracks nucleation as soon as the stress exceeds the local rupture threshold of the material.

Therefore, the self-heating associated with crack propagation is expected to play a crucial role in material failure. Indeed, it has been estimated that in elastomers [20] and polymers [8] such increase of temperature could modify the fracture energy (leading to its non-monotonic behavior with crack growth velocity) and finally cause a “Stick-Slip” rupture instability. The temperature rise has also been shown to significantly decrease the yield stress in steel [21], metallic glass [22] or quasicrystals [23], where it is associated to the formation of dimples at the fracture surface [23, 24]. In Earth sciences, based on the observation of pseudotachylite (melted crushed rock by friction during faults dynamics), a recent theoretical approach proposed that self-localizing thermal runaway processes [25] could be responsible for deep earthquakes [26]. The reactions triggered by heating also play a major role on rock mechanics deformation during Earthquake faulting [27, 28], or in general during metamorphism and the localization of deformation [29].

Other features of fracture seem rather related to non thermal, but rather quenched disor-
der, i.e. heterogeneities of the material properties. The distribution of disorder influences the
localization pattern and mode and is responsible for the intermittency and avalanches [30–35].
In general, this intermittency is responsible for crackling noise [36] associated to fracture in
disordered materials, as e.g., in paper during crumpling [37] or fracturing [38].

To understand the physical mechanisms of material failure and the possible instabilities due to the feedback of the temperature on the fracture propagation, it appears crucial to be able to establish an accurate energy budget of crack growth, from a clear identification and measure of the various dissipation sources.

In general, from the total change in mechanical energy during propagation, a fracture can create new surfaces and defects in a process zone (dislocations or irreversible internal energy changes). It can also emit acoustic emissions that dissipate further away [38, 50], it can radiate light [5, 50], and generate heat during friction and Joule effect [45]. The conversion of mechanical energy into heat is sometimes assumed to be total in models [3], though it has been shown experimentally that only a fraction of it is converted [21, 22]. In the materials probed in these experiments in metals and metallic glasses, this conversion into Joule heating through plastic work is still close to complete. We will show here that for fracture in paper, the conversion
efficiency into heat is much more inefficient, and most of the energy goes into defect formation and creation of new crack surfaces.

Here, we report experimental measurements of surface temperature obtained from an infrared (IR) camera during slow tensile fracture in paper. We implement a simple physical model for the heat flow, taking into account in plane diffusion, losses towards the surrounding air, and Joule heating in a process zone around the crack tip. Performing a numerical integration, we notice that a simplified quasi one dimensional model with lateral losses for heat propagation renders for the thermal evolution observed experimentally. This allows to formulate an invariant, based on integrals of the temperature profile perpendicularly to the crack propagation, that leads to an estimate of the energy brought by Joule heating during the process. Comparing the predicted shapes of the temperature fields to the measured ones, one can evaluate the fraction of the mechanical energy converted into heating during the fracture process, found to be around $\alpha \simeq 12 \pm 4\%$ in our experiments. We also show that this fraction $\alpha$ rises significantly during fast propagation stages, increasing with velocity, from $12\% \pm 4\%$ to about 40% as the crack speed jumps from 1 mm/s to 18 mm/s – while no significant change of the energy release rate is recorded.

Finally, we discuss the temperature reached at the hottest points, and show from the dimensions of the fibers, that it can reach above 200°C, compatible with a thermal decomposition of a few paper fibers, reported to auto-ignite around Fahrenheit 451 (Bradbury, 1953 [39]). We show that the heat accumulation allows a runoff of the fracture at a speed of a few centimeters per second set by heat diffusion in the process zone, in accordance with the measured speeds.

# 1 Method

## 1.1 Experiments

We have monitored the slow growth of a single crack in paper sheets submitted to a pure tensile loading. Actually, the fracture of such material has attracted a lot of attention in the statistical physics [11–13, 19, 38, 40, 41] and material science [42–45] communities.

Samples of regular laser printing paper of dimensions $L_w \times L_l = 10 \times 10$ cm$^2$ and thickness $h = 112 \pm 20$ µm (measured by optical techniques) are fixed to an external frame consisting of two cylindrical rods.

These samples are pre-cut on one side with an initial defect of length $l_{p.c.}$ between 1 and 2 cm, in order to have a better control of the rupture process (localized at the crack tip). The width of the intact paper sample, after subtracting this precut length, is $L'_w = L_w - l_{p.c.}$. They are submitted to elongation at a constant speed of 40µm/s, while the longitudinal elongation is recorded, as well as the total force exerted across the paper. An infrared camera, initially calibrated for temperature measurements, records the paper surface temperature, with 320x80 pixels of linear size 0.3 mm, at a rate of 150 frames per second. Simultaneously, a fast camera (Photron SA4) records optical images of the same field at 250 frames per second, with a resolution of 1024×1024 pixels of linear size 0.096 mm per pixel. The longest direction of the fibers is roughly perpendicular to the elongation direction.

After a few mm of elastic loading, a crack starts to propagate slowly, the temperature rising around the tip. A snapshot of the moving crack, with the IR image superimposed on the corresponding optical image, is shown in Fig. 2. A video recorded with the IR camera is provided as supplementary data. In Fig. 3, four instantaneous temperature maps are displayed: after a peak in temperature, preceding the apparently open section of the crack, a trail of elevated temperature is observed. The dynamics of the crack tip, identified as the hottest spot in the image, is displaying a varying longitudinal velocity, up to a few cm/s, as illustrated in Fig. 3.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Scheme of the experimental setup, paper sample and control device: side view (top) and front view (bottom). The paper sheet is precut (black line) over 1 to 2 cm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Snapshot of the crack propagation, with a colored image of the IR camera superimposed with the corresponding optical impage. The heated zone precedes an open crack. The black dots on the paper are patterns that were printed for deformation measurement and image analysis purposes.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top: Successive configurations of the IR image of the crack propagation. (a) and (b) correspond to a stable propagation of the crack tip around 1 mm/s, (c) and (d) to a final fast jump around 10 to 20 mm/s. The characteristic width of the hot zone, along the Y direction, is a few mm, and the characteristic length, along the X direction, is a few mm in (a,b), and a few cm in (c,d). Movie in supplementary material. (e) Longitudinal position and speed of the crack tip (hottest point) as function of time. (f) Tensile force versus extension of the sheet during the experiment.</figcaption></figure>
The force measurement during the intermittent slow crack growth allows to compute the total work W done, integrating the force displacement curve, shown in Fig. 3. Dividing this work by the cross section of the paper, $S = L'_w \times h = 8 \text{ cm} \times 112 \mu\text{m}$, the work per nominal surface area is found to be $W/S = G = 17 \text{ kJ m}^{-2} \pm 5 \text{ kJ m}^{-2}$ for the various experiments performed (it displays a non negligible variability from one paper sample to another). Two characteristics of the paper samples are independently measured during different calibration experiments, detailed in the supplementary material: the in-plane heat diffusion coefficient, $D \simeq 4.4 \cdot 10^{-8} \text{ m}^2 \text{ s}^{-1}$ and the decay rate due to loss into the surrounding air of the out-of-plane thermal flux, formulated as $(T - T_{air})/\tau$, with $\tau \simeq 5.2 \text{ s}$, where T and $T_{air}$ are respectively the temperature averaged over the paper thickness, and the surrounding local air temperature.

## 1.2 Theory

During the fracturing process, the energy brought by the external device to the sample is transformed into reversible elastic energy, and some irreversible energy losses. The difference between the work done $dW$ and the stored elastic energy $dE$ is proportional to the nominal new fracture surface created $dS$, with a factor of proportionality called the energy release rate $G$: $dW - dE = GdS$. The irreversible energy losses can be divided in:

1. Energy irreversibly stored into new surface and defects created in the process zone (such as dislocations).

2. Energy radiated under mechanical waves emitted both in the solid and in the surrounding atmosphere as sound waves dissipated far from the crack [38, 41]. It is usually a very small fraction of the mechanical energy, between $10^{-5}$ and $10^{-2}$ for fracture in glass [46–48] or earthquake sources [49]. We can thus neglect this term in the energy budget.

3. Energy emitted in the form of electromagnetic waves, such as light emissions [50]. This third form of energy losses, estimated to a fraction around $10^{-5}$ [50] can be also neglected in the global energy budget.

4. Energy dissipated during friction of fibers sliding against each other. This small scale vibrational energy corresponds to Joule heating, and contributes to a rise in temperature $dT$, with an internal energy density increase $\rho cdT$, where c is the specific heat of paper (in J.K$^{-1}$kg$^{-1}$), and $\rho$ its mass density. This Joule heating energy density is also equal to $\alpha GdS$, where $\alpha \in [0, 1]$ is the fraction of the irreversible energy contributing to heating.

Heat due to the frictional process is created locally around the crack tip, in the process zone. The form of this heat source is thus, per unit volume, $dw = \alpha GdS f(x, y, t)$, where $(x, y)$ denote the coordinates along the average fracture propagation direction and along the imposed elongation, and $t$ the time. $f(x, y, t)$ is the support function of the process zone, of normalized spatial integral (i.e. of units m$^{-3}$). It is centered around the crack tip at $(x_p(t) = x_0 - vt, y_p(t) = y_0)$, assumed for simplicity to be moving at constant speed $v$ in the $x$ direction, and of linear extension $l$, corresponding to the radius of the process zone (considered circular for simplicity). The open fracture surface during time $dt$ is $dS = hvdt$, with $h$ the paper thickness, and $v$ the crack tip velocity. Eventually, the heat diffuses according to Fourier law, both in-plane, and out-of-plane with losses into the surrounding atmosphere. One can write the three dimensional (3D) heat flux in the paper bulk as $\mathbf{j} = -\lambda \nabla T$, where $\lambda$ is the thermal conductivity. Along the paper surfaces, assuming a linear process for the coupling to the atmosphere [51], one can write the perpendicular heat flux to the atmosphere as $j_{\perp} = -\kappa(T - T_{air})$, where $\kappa$ is a heat transfer coefficient [52]. The 3D conservation equation of energy $\rho cdT = -\nabla \cdot \mathbf{j} dt + dw$, leads to the following two dimensional form after integration across the thickness of the paper:

$$ \partial_t \Delta T = D \nabla^2 \Delta T - \Delta T / \tau + \alpha G f(x, y, t) hv / (\rho c) \quad (1) $$
where $\Delta T(x, y) = (T - T_{air})$ is the temperature elevation through the paper at position $(x, y)$ (variations through the plane, along $z$, are negligible), $\nabla^2$ in the inplane Laplacian, $D = \lambda/(\rho c)$ is the thermal diffusivity, and $\tau = \rho ch/(2\kappa)$ is the characteristic damping time due to lateral thermal losses in the surrounding air. This is a diffusion equation, with a damping term, and a source coming from Joule heating during the crack propagation.

# 2 Results: evaluation of the energy fraction spent in Joule heating

## 2.1 Theory: numerical observations and analytical approximations

We solve equation 1 using an Alternate Direction Implicit algorithm [53]. The results are checked on some examples by comparison to the explicit direct numerical integration of the Green function of diffusion with analytically calculated loss term (see supplementary information):

$$ \Delta T(x, y, t) = \int_{0}^{t} ds \iint d\xi d\eta \frac{\alpha G h v f(\xi, \eta, s)}{\rho c \cdot 4\pi D(t-s)} \exp \left[ -\frac{(x-\xi)^2 + (y-\eta)^2}{4D(t-s)} - \frac{t-s}{\tau} \right] \quad (2) $$

After a short transient stage, the simulated temperature field is stationary (in the reference frame of the moving process zone), and a typical simulation result is displayed in Fig.4. This transient stage vanishes as soon as $\sqrt{Dt}$ exceeds a couple of times the source dimensions, i.e. for $t > 2l^2/D$, which corresponds to microseconds up to a fraction of second depending on the size $l$ in the simulations. The temperature field at scales above 0.4 mm is not affected by the process zone size $l$, if $l < 100 \ \mu m$ (See e.g. supplementary material comparing different simulations). This size only affects the temperature in the vicinity of the tip, at smaller scales.

A characteristic dimension corresponds to the skin depth of in plane thermal diffusion before the temperature reaches a characteristic background level due to lateral temperature losses, i.e. to out of plane diffusion with a characteristic time $\tau$. This is given by $l_{hot\perp} = \sqrt{4D\tau} \approx 0.9$ mm in the case of paper. This corresponds indeed to the observed extent normal to the crack propagation of significantly heated region around the crack tip, visible on Fig. 3. Another size, along the crack propagation, can be computed: $l_{hot\parallel} = v\tau \approx 5$ mm to 5 cm, for cracks propagating around 1 mm/s to 1 cm/s. This is order of the extent observed in Fig. 3 (a,b) in the slow propagation regime, or (c,d) in the fast one. These sizes, sometimes referred to as those of a plastic deformation zone in IR observations of cracked paper [42], correspond in fact to the significantly heated region, but is not particularly associated to plastic deformations: it exceeds by far the process zone size $l$, where the significant heat source $f$ is distributed.

Apart from the temperature in the tip vicinity, for process zones of size $l = 100 \ \mu m$ or below, the size of the process zone does not affect the far temperature field, at the measurement scale of the IR camera (pixels of 0.4 mm).

One notices that the temperature gradient behind the crack tip (in the “crack tail”) is almost perpendicular to the crack propagation. The temperature profiles along the y-direction are also found to be Gaussian up to a few % in this tail (see supplementary material). We observe that $\partial_{xx}T \ll \partial_{yy}T$ for zones behind the tip, $x > x_p$ where $x_p$ is the abscissa of the crack tip, centre of the process zone, moving towards low $x$ direction in the convention chosen. Consequently, Eq. (1) can be simplified there as:

$$ \partial_t \Delta T = D \partial_{yy} \Delta T - \Delta T / \tau + \alpha G f(x, y, t) h v / (\rho c) \quad (3) $$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Simulated temperature around a crack propagating at 1 cm/s, with $\alpha G = 2kJ/m^2$, in accordance with the parameters inferred from the experiments. The size of the process zone here is below the resolution size, and has not impact on the figure - i.e. it corresponds to a source function $f$ corresponding to a spatial Dirac function, with a size $l \to 0$. The explicit numerical integration was used here.</figcaption></figure>

Since the dimensions of the process zone are not felt at the mm scale and above, the geometrical source term can be simplified, conserving its weight and location, as $f(x, y, t) = \delta(x - x_p(t))\delta(y)/h$, where $x_p(t) = x_0 - vt$ and $\delta$ is a Dirac distribution. Hence, Eq. 3 is a one dimensional diffusion equation, independent of $x$, apart implicitly to set the time $t_p(x) = (x_0 - x)/v$ when the crack tip passes at $x = x_p(t_p)$. The crack deposits an energy $\alpha G$ per unit area at the moment when it passes. The rest of the time, no source term is present. The solution without the loss term of this one dimensional diffusion equation would conserve the energy along a y-profile. The loss term modifies this, introducing an exponential decay of the energy deposited, but this can be corrected for, so that the energy $\alpha G$ deposited by the crack in Joule heating, that diffuses along any y-profile and can be determined by spatial integration. Technically, the constant multiplication method amounts to introduce the variable: $\Delta T' = \Delta T e^{(t-t_p)/\tau}$ and consider its time derivative (see supplementary information). This quantity satisfies

$$ \partial_t \Delta T' = D \partial_{yy} \Delta T' + [\alpha G h v / (\rho c)][\delta(t - t_p(x))/v]\delta(y)/h, \quad (4) $$

with initial conditions $\Delta T'(x, y, t < t_p) = 0$. Its solution is obtained from the Green function of one dimension diffusion $G_0(y, t')$, as $\Delta T'(x, y, t) = \alpha G / (\rho c) G_0(y, t - t_p)$. The temperature elevation can thus be approximated, noting $\Theta$ the Heaviside function, as

$$ \Delta T = \frac{\alpha G}{\rho c \sqrt{4\pi D(t - t_p)}} e^{-\frac{y^2}{4D(t - t_p(x))}} e^{-\frac{t - t_p(x)}{\tau}} \Theta(t - t_p(x)), \quad (5) $$

Integrating this along a y-profile behind the crack tip (i.e. at any fixed $x$ so that $t > t_p(x)$, and correcting for the exponential decay of temperature to the surrounding, one obtains an invariant corresponding to the Joule energy deposited by the process:

$$ I(x, t)e^{(t-t_p(x))/\tau} = \frac{\alpha G}{\rho c}, \quad (6) $$

$$ \text{where } \int_{-\infty}^{\infty} dy \Delta T(x, y, t) = I(x, t) \quad (7) $$

## 2.2 Evaluation of the Joule heating term in the experiments

In the experiments, the crack speed fluctuates, and the crack tip motion is not perfectly straight. However, the analytical solution, at a position $(x, y)$ and time $t$ in Eq. (5), does not depend on the speed, but only on the time elapsed since the crack tip crossed, $t - t_p(x)$, and $y$ represents $y - y_p(x)$, the distance along the y-direction to the point where the crack tip crossed, $y_p(x)$ (which was assumed zero in the analytical solution for simplicity). To avoid fluctuations and
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Space time diagram of the experimentally observed temperature increase behind the tip, as function of the y-coordinate, with respect to the position of the tip when it passed along this profile, and the time elapsed since the tip crossed the profile.</figcaption></figure>

obtain better statistics, we transform the instantaneous spatial map of temperature into a space time map of temperature increase, i.e. for each $x$, for all $t > t_p(x)$ and all $y$, we represent $\Delta T$ as function of $(y - y_p(x), t - t_p(x))$.

This change of reference frame allows to account for the process intermittency. After averaging these space time maps over all $x$, we obtain the map displayed in Fig. 5.

Cutting at constant $t - t_p$ (vertical lines in Fig. 5) leads to the temperature profiles shown in Fig. 6, top. These correspond well to Gaussian profiles of growing width and decaying prefactor, as predicted by Eq. (5). The integral $I$ of these excess temperature profiles over $y$, on each available positive $t - t_p$, leads to the measures shown on the top Fig. 6 in semilogarithmic scale. The straight line is predicted by Eq. (7), i.e. $I = (\alpha G)/(\rho c)e^{-(t-t_p(x))/\tau}$, with a slope $1/\tau$ corresponding to the expected $\tau = 5.2$s from the calibrations. The prefactor of this law, at $t - t_p = 0$, allows to measure the deposited energy in the form of Joule heating: it corresponds, in this experiment, to $\alpha G/\rho c = I_1 = 2.5 \cdot 10^{-3} \text{ K} \cdot \text{m}$, for the slow moving part (dark line), and to $\alpha G/\rho c = I_2 = 8.7 \cdot 10^{-3} \text{ K} \cdot \text{m}$ for the fast stage (cf Fig. 6, bottom, inset and top points).

From the characteristics of paper, $c \simeq 1000\text{J} \cdot \text{kg}^{-1} \cdot \text{K}^{-1}$, and $\rho \simeq 800\text{kg} \cdot \text{m}^{-3}$, and thus its volumic heat capacity is $\rho c \sim 8 \cdot 10^5 \text{J} \cdot \text{m}^{-3} \cdot \text{K}^{-1}$, one can express $\alpha G = I_1 \rho c = 2000\text{J}/\text{m}^2$ during most of the experiment. With the determined value, $G = 17\text{kJ} \cdot \text{m}^{-2}$, this allows to estimate the ratio sought for, $\alpha = I_1 \rho c/G = 0.12$. Over 7 experiments where temperature, force and displacement were monitored, this Joule heating efficiency was found to be $\alpha = 0.12 \pm 0.04$. During fast crack growth events, the energy deposited in Joule heating jumps to $\alpha G = I_2 \rho c = 6900\text{J}/\text{m}^2$, and thus to a ratio $\alpha \simeq 0.4$. The occurrence of this fast stage is a common feature of the experiments performed, as is this increase of $\alpha G$. This fast stage does not systematically happen at the end of the experiment, in some of them the velocity drops down again to the previous low values before the paper is entirely cracked. Such an example is analyzed in details in Supplementary material, Fig. 13.

## 3 Discussion and considerations on the crack tip temperature

As shown above, the whole energy released during the fracture does not contribute to Joule heating, but some significant fraction of it does. A priori, a weakening mechanism due to the feedback from temperature is therefore possible. However, in general fracture problems, caution is needed regarding the exact amount of energy available as a heating source. For example, for evaluations of heating in plastic material, models [3] consider that the whole energy is available as a Joule effect, i.e. $\alpha = 1$. This result can be modified in the case of materials presenting a much smaller heating efficiency, i.e. $\alpha \ll 1$. In the case of the paper studied here, a large fraction (88 %) of the mechanical energy goes to non heating damage (creation of defects and isolated dislocations in the process zone) and fracture energy (surface energy). Indeed, creating
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Experiments: Top: Profile of temperature along the y-axis (continuous lines), and Gaussian fits (dashed ones). Bottom: Integral of the temperature elevation along profiles as function of the time elapsed since the crack passed. The dotted line has a constant slope corresponding to the theoretical prediction of exponential decay. The regression to t=0 gives the total amount of Joule energy brought by the fracture process in a given profile. The long dark curve and points refer to the stable propagation stage, with $I(t = 0) = I_1 = \exp(-6) = 0.0025K \cdot m$, the points above (and the zoom in the inset) to the final second of the experiment, configurations (c) and (d) in Fig.3, where the velocity jumped above 1 cm/s – with a corresponding $I(t = 0) = I_1 = \exp(-4.75) = 0.0087K \cdot m$. The gray lines correspond to two subsets of the experiment in slow propagation periods with a tip velocity around 1 mm/s, and are representative of the variations in this slow stage across experiments.</figcaption></figure>

new surface, defects or dislocations is not generating heating (and a temperature change) per se, but consumes energy by changing the potential energy in the internal energy – a contribution corresponding to the strain energy for the dislocations, or to the surface energy of the new surface created for the fracture surface. A contrario, the phonons (disordered vibrations) that are emitted during this creation, associated to a change in internal energy, notably via some extra kinetic energy at small scale in the vibrational degrees of freedom, i.e. to the deviation of velocity from the mesoscopic average of velocity, correspond to the heating.

The maximum crack tip temperature is not accessible due to the limited IR camera spatial resolution. In the case of paper, two process zone sizes can be considered, associated to two different processes [40]: First, breaking fibers and heating at this scale, which corresponds to a size $l \simeq 1\mu m$ to $10\mu m$, the order of magnitude of the cross section of the fibers. Next, disentangling fibers around the tip, which generates friction over a size around $l \simeq 10$ to $500\mu m$, the length of the fibers. These scales are the main candidates for the process zone size – although it is unclear which one of the two processes, breaking or disentangling fibers frictionally, dominates the Joule heating source.

Another typical size can be introduced: On any point around the crack trajectory, the process zone is present for a time $T = l/v$, during which the heat diffuses out of the process zone over a characteristic skin depth

$$ \delta = \sqrt{4DT} = \sqrt{4Dl/v}. $$

For a typical velocity of 1 mm/s to 1 cm/s, this corresponds to sizes in the range of $\delta \simeq 4$ to $13\mu m$ if $l = 1\mu m$, $\delta \simeq 12$ to $40\mu m$ if $l = 10\mu m$, and to $\delta \simeq 40$ to $130\mu m$ if $l \simeq 100\mu m$. These orders of magnitude of $\delta$ are close to the process zone size $l$.

Different cases of temperature fields, simulated by solving the heat flux in the tip, are illustrated on Fig. 7, and in the supplementary material, Figs. 3-11, for process zone sizes ranging between 10 and $100\mu m$, and for crack speeds between $v = 0.0001$ to $0.1$ m/s.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Simulated temperature elevation (in K) during the stationary stage for a crack at $v = 1 \text{cm/s}$ with a circular process zone of size $l = 10 \mu\text{m}$, limited by the dash dots. On the back (left) of the process zone, an elevation around $200^\circ \text{K}$ is reached over a $\simeq 1 \mu\text{m}$ size zone. Added to a $20^\circ \text{C}$ background temperature, this corresponds to the top most point at corresponding speed in Fig. 8 - close to the autoignition temperature of paper.</figcaption></figure>

The process generates a heating energy $\alpha Glh$ over a process zone, assumed as roughly cylindrical around the crack tip, of volumic size around $r^2h$, during a time $l/v$. The positions along the trajectory of the tip (the fracture line) is heated for the longest time, $t = l/v$.

One can consider three limiting cases:

Either the crack propagation is slow enough, so that heat diffuses efficiently out of the process zone in all directions, i.e. $\delta \gg l$, or

$$ \sqrt{4Dl/v} \gg l, \quad \text{i.e.} \quad v < v_c = 4D/l. $$

In this case, the corresponding energy spreads approximately as a Gaussian around the center of the process zone (the tip), over a cylindrical zone of cross sectional area $\pi(2\delta)^2$. The corresponding temperature rise corresponds then to

$$ \begin{aligned} \Delta T_{tip} & \simeq (\alpha Glh/(4\pi\delta^2h))/(\rho c) \qquad (8) \\ &= \alpha Gv/(\pi\rho cD). \end{aligned} $$

Another possibility is that the skin depth of diffusion exceeds the lateral distance between the crack tip and the side of the zone $l/2$, but not the full length of the zone $l$. In this case, the heat diffuses laterally over a length $\delta$ on two sides, but the diffusion does not allow the heat to cross entirely the process zone in the direction of propagation of the crack. This amounts to distribute the energy $\alpha Glh$ over a zone of approximate size $\delta lh = \sqrt{(4Dl/v)lh}$. This is the case assumed for the propagation of crack in glass, in [3]. Such a case is shown in supplementary material, Fig. 8, for $v = 1 \text{mm/s}$, $l = 100 \mu\text{m}$. The corresponding temperature rise corresponds then to

$$ \begin{aligned} \Delta T_{tip} & \simeq (\alpha Glh/(\delta lh))/(\rho c) \qquad (9) \\ &= \alpha G/(\rho c\sqrt{4Dl/v}) = [\alpha G/(\rho c)]\sqrt{v/4Dl}. \end{aligned} $$

Or eventually, for a crack propagating fast enough, the heat does not diffuse significantly out of the process zone during the crack motion. This is characterized by $\delta < l$, i.e. $v > v_c = 4D/l$. For a process zone size around $l \sim 10 \mu\text{m}$, this corresponds to a criterion $v > v_c \simeq 18 \text{mm/s}$. In
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Maximum temperature as function of the crack speed, evaluated from the approximations of Eqs. ((8,9,10)) (curves, for process zone sizes $l = 1, 10, 20, 50$ and $100\mu m$). The symbols, that can be compared to these proxys, are the maximum recorded temperature over all positions, in simulated temperature fields at $l = 10, 20, 50, 100\mu m$, at different speeds. The dash-dotted horizontal line represents an elevation of 200 ° K, corresponding to an oxidization temperature (paper auto-ignition for macroscopic material) of 220 ° C for a room at 20 ° C. The temperature saturates above $v = v_c$. For $l = 10\mu m$, this is reached at a few cm/s.</figcaption></figure>

this case, the whole energy $\alpha Glh$ stays in a zone of dimensions comparable to the process zone size $l^2h$ while the crack passes, leading to a temperature rise of

$$ \Delta T_{tip} \approx \alpha Glh/(l^2h)/(\rho c) = \alpha G/(\rho c l). \qquad (10) $$

With possible process zone sizes $l = 1, 10$ or $100\mu m$, this corresponds to potential temperature rises respectively of 1300 to 130°C, or 13°C. Such temperature rises, for energy sources in the process zones at sizes close to fiber cross section, at 1 to $10\mu m$, could lead to a temperature exceeding the autoignition of paper in this micrometric zone, around 220°C, i.e. Fahrenheit 451 [39].

On Fig. 8, we show these estimates of the maximum temperature, corresponding to the Eqs. (8,9,10), as function of the crack speed, for possible characteristic process zone sizes, $l = 1, 10, 20, 50$ and $100\mu m$. The autoignition temperature is also shown, supposing a laboratory temperature at 20° C. The curves represent the estimates of the approximations from Eqs. (8,9,10) above. The symbols represent the hottest temperature recorded in the stationary stage in numerical simulations, shown in Fig. 7 and in Figs. 3-11 in the supplementary material. This hottest temperature is located around the center of the process zone at low speed (see Fig. 9-10 in supplementary), or at large speed around the back of the process zone (in the part of the process zone lagging behind with respect to the propagation direction). This is clearly seen in Fig. 7, where the hottest point is on the left of the process zone, marked as a dash dotted curve – and on Fig. 3 of the supplementary material. On Fig. 8, the agreement between the analytical approximations Eqs. ((8,9,10)) of the maximum temperature and the simulation results seems satisfactory.

The analysis of the experimental crack velocity typically shows arrests of the crack, followed by fast crack propagation maintained up to the final rupture. This is seen on Fig. 3, around $t = 84s$, where the crack jumps to a velocity around $18mm/s$ for close to $1s$, and tears the rest of the paper apart. We note that this velocity is of the same order as the critical velocity $v_c = 4D/l$, which corresponds to the velocity obtained by allowing diffusion of heat over a skin depth equal to the process zone size. This observation allows to formulate the following hypothesis: when the process zone is small enough and the speed large enough, the temperature reached on the back of the process zone, behind the crack tip, can reach autoignition temperature. In this case, some energy is produced by the oxidization reaction in this part of the process zone, in addition to the one produced by Joule heating converting mechanical energy everywhere in the zone. This additional energy of chemical origin can reach the head of the process zone after the heat
is transported by diffusion through the process zone. This allows to increase the temperature of the fibers ahead and weaken them, leading to the fracture propagation. The speed at which this happens is set by the time of heat diffusion from the back to the head of the process zone, and corresponds to the observed velocity $v_c = 4D/l$ during the fast stages. Indeed, at such speed, for a process zone around $l = 10\mu m$, a part of the process zone reaches the temperature of paper autoignition – see Fig.7.

It is thus possible that the temperature plays an active role (i.e. affects the velocity) during this last propagation stage for the crack propagation, and is not simply a "passive marker" of the mechanical process via Joule heating. We also observe during this last stage, (Fig. 3 (c) and (d)), that the temperature observed is significantly larger, and the estimate of Joule energy release rate $\alpha G$ is more than 3 times larger than its characteristic value during the stable slow propagation - as can be seen on Fig. 6, bottom and inset: $ln(I(t = 0)/I_0)$ for this stage (blue circles around time $t = 0$) lies around 1.2 above the corresponding value for the rest of the propagation (dark full line), i.e. $\alpha G \propto I(t = 0)$ increases by a factor $e^{1.2} \sim 3$ during this stage. The mechanical energy release rate, $G$, is not changing significantly during this stage (see supplementary material, section V and Fig. 14). This increase in Joule energy release rate can be attributed potentially to a mechanism increasing significantly the share of the mechanical energy release contributing to the Joule heating, $\alpha$, from 12% to roughly 40%. We observe in different experiments that this increase happens at fast velocities, but anywhere in the paper sample, not particularly close to the boundaries - see an example in supplementary material where a fast propagation stage and simultaneous increase of heating efficiency happens far from the boundaries (Figs. 13 and 14) – It should thus not be an effect due to the mechanical interaction between the tip and the boundaries, but rather an effect related to the velocity of propagation, that would change the ratio between the Joule heating and the damage creation energy, for a fixed mechanical energy released.

Another possibility for this increase of heating ratio $\alpha$, if indeed the temperature comes close to the autoignition temperature, is that this increase of heating results from an overall increase of energy released with an extra non mechanical energy source, related to the exothermic character of the reaction (cellulose oxidization, $C_6H_{10}O_5 + 6O_2 \rightarrow 6CO_2 + 5H_2O$). Paper combustion of the process zone, indeed, brings potentially an enthalpy of $h_b \approx 10^4$ kJ/kg of cellulose [54], i.e. an extra energy source of $G_b = h_b\rho l = 80\text{kJ/m}^2$, with $l = 10\mu m$, which is roughly 4 times larger than the mechanical energy release rate $G = 17\text{kJ/m}^2$. During the slow stage, the Joule heating energy release rate is $\alpha G \sim 2\text{kJ/m}^2$. Suppose for simplicity that this stays fixed in the fast stage, and denote $\alpha_0 G$ this part coming solely from converting mechanical energy into heat. Then, oxidizing effectively a ratio $\alpha' = 5\%$ of the mass of the fibers in the process zone could thus be sufficient to provide the observed increase by a factor 3 to 4 of the total heat source, to reach in total a heating of $G_{heating} = \alpha G = \alpha_0 G + \alpha'G_b = 2\text{kJ/m}^2 + 0.05 \times 80 = 6\text{kJ/m}^2$ – which corresponds to the observed heating $\alpha G = 0.4 \times 17\text{kJ/m}^2 = 6\text{kJ/m}^2$ observed during the fast stage, in contrast with $\alpha G = 2\text{kJ/m}^2$ during the slow stage. In Fig 7, indeed, only a small fraction of the area of the process zone reaches temperature rises around 180 °C, around its back. Oxidizing some fibers in this small region, roughly micrometric in size, corresponding to oxidizing 5 % of the process zone, could provide the required extra energy. This process can be triggered initially due to a smaller or harder fiber met by the fracture, with a resulting smaller process zone or acceleration of the crack which lagged for a while behind the equilibrium position, allowing the back of the zone to reach this oxidization temperature. Once this is triggered, the extra heat generated by this reaction at the back of the process zone is transported by heat diffusion. When it reaches the head of the process zone, it can fragilize the next intact fibers lying in front of the crack. This process should happen at a speed set by a diffusive limiting factor, i.e. by the time $4D/l^2$ required to transport the heat by diffusion through the whole size $l$, i.e. at the speed $v_c = 4D/l \sim 18\text{mm/s}$, which is indeed of the order of the observed speed during the fast crack propagation stages.
In the proposed mechanism, it is required that the process zone reduces to the thickness of the fibers size, rather than the length of the fibers, to reach the autoignition temperature locally, and cause the subsequent fiber softening and dramatic crack acceleration at the speed selected by diffusion. For this to happen, it is sufficient that a group of fibers gets sufficiently entangled to form a knot, and mobilize a friction high enough to block the normal slip mechanism, which is certainly rather fiber disantanglement (happening at sizes around the fibers length, 10 to 100 µm). In the presence of such knot acting as a pinning point, the velocity can presumably after breaking the fibers composing the knot, jump to a speed of the order of cm/s due to the overloading to break these fibers, and the successive reduction of the stress threshold. Thus, in the proposed scenario, the presence of such pinning fiber configurations allows to fluctuate to higher velocities and smaller process zone size, reaching values of (l, v) around (10µm, 1cm/s), or (l, v) around (1µm, 3mm/s), and obtain local temperatures at the back of the process zone large enough to trigger the oxidization of cellulose (see corresponding points and curves on Fig. 8). The later transition from slow to fast propagation regime, due to extra heat generated by oxidization of cellulose, is the consequence of both the velocity fluctuations and reduction of process zone size at the meeting between the crack tip and a fiber knot/strong fiber.

In conclusion, we have studied the temperature of paper during crack propagation under tensile failure driven at slow constant extensional speed. Using a simple theoretical description of the inplane and out of plane diffusion, we have shown that this temperature measurements allow to determine the amount of Joule heating. It was measured that 12±4% of the mechanical energy release rate is converted into heating, the rest being spent in damage and surface energy and radiation of mechanical waves. The temperature at the crack tip was determined, depending on the process zone size and the crack speed. We have shown that, for process zones under 10 µm, when the crack velocity comes close to $v_c = 4D/l$, temperature increases could locally reach the autoignition temperature in the process zone. We have also observed that the crack tip velocity can jump to a significantly larger speed, during which the Joule heating release rate increases by more than 200%. We have shown that our temperature measurements are compatible with the trigger of an additional non mechanical energy release, coming from the oxidation of a fraction of the fibers in the process zone. The observed velocity in this regime is also compatible with a mechanism where the heat diffusion plays an active role in transporting heat and softening paper.

The authors acknowledge fruitfull discussions with Signe Kjelstrup, Dick Bedeaux, Alex Hansen, Eirik G. Flekkøy, Mikko Alava, Elisabeth Bouchaud and Stefan Nielsen, and support from the Norwegian Research Council FRINATEK project 205486, the University of Strasbourg IDEX, the University of Lyon, the ENS Lyon, the CNRS, the regional alsatian project REALISE, the ITN FLOWTRANS european network, and the ANR LANDQUAKE.

## References

[1] Griffith, A. A., *Philosophical Transactions of the Royal Society of London*, 1921, **A 221**, 163–198.

[2] Irwin G., 1957, *Journal of Applied Mechanics*, **24**, 361–364.

[3] Rice J. R. & Levy N., in *The Physics of Strength and Plasticity*, edited by Argon A. S. (The M.I.T. Press, Cambridge, Mass.; London, England) pp. 277–293 (1969).

[4] Fuller, K.N.G., Fox, P. G. & Field, J. E., *Proc. R. Soc. Lond. A*, 1975, **341**, 537–557

[5] Pallares, G., Rountree, C. L., Douillard, L., Charra F. & Bouchaud, E., *Europhysics Letters*, 2012, **99**, 28003. doi: 10.1209/0295-5075/99/28003
[6] Ben-David, O. & Fineberg, J., *Tribology Letters*, 2010, **39**, 235-245.

[7] Ben-David, O., Rubinstein, S. M. & Fineberg, J., *Nature*, 2010, **463**, 76-79, doi:10.1038/nature08676,

[8] Williams, J.G., *Int. J. frac. Mech.*, 1972, **8**, 4, 393-401

[9] Brenner, S.S., *J. Appl. Phys.*, 1962, **33**, 33. doi: 10.1063/1.1728523

[10] Zhurkov, S.N., *Int. J. Fract. Mech.*, 1965, **1**, 311-323, reprinted in 1984, *Int. J. Fract.* 26, 295-307.

[11] Santucci, S., Cortet, P.-P., Deschanel, S., Vanel, L. & Ciliberto, S., *Europhys. Lett.*, 2006, **74**, 595.

[12] Santucci, S., Vanel, L. & Ciliberto, S., *Eur. Phys. J. ST*, 2007, 146.

[13] Vanel, L., Ciliberto, S., Cortet, P.-P. & Santucci, S., *J. Phys. D : Appl. Phys.*, 2009, **42**, 214007.

[14] Lengliné, O., Toussaint, R., Schmittbuhl, J., Elkhoury, J.E., Ampuero, J.-P., Tallakstad, K.T., Santucci, S. & Måløy, K.J., *Phys. Rev. E*, 2011, **84**, 3, 036104. doi: 10.1103/Phys-RevE.84.036104

[15] Pomeau, Y., *C.R. Acad. Sci. Paris II*, 1992, **314**, 553.

[16] Roux, S., *Phys. Rev. E*, 2000, **62**, 6164.

[17] Scorretti, R., Ciliberto, S. & Guarino, A., *Europhys. Lett.*, 2001, **55**, 626.

[18] Santucci, S., Vanel, L., Scorretti, R., Guarino, A. & Ciliberto, S., *Europhys. Lett.*, 2003, **62**, 320.

[19] Santucci, S., Vanel, L. & Ciliberto, S., *Phys. Rev. Lett.*, 2004, **93**, 095505.

[20] G. Carbone & Persson, B. N. J., *Phys. Rev. Lett.*, 2005, **95**, 114301.

[21] Zhender, A.T. & Rosakis, A.J., *J. Mech. Phys. Solids*, 1991, **39**, 3, 385-415.

[22] Flores, K.M. & Dauskardt, R.H., *J. Mat. Res.*, 1999, **14**, 3, 638-643

[23] Ponson et al. L. Ponson, L., Bonamy, D. & Barbier, L., *Phys. Rev. B*, 2006, **74**, 184205.

[24] Wang, G., Chan, K.C., Xu, X.H. & Wang, W.H., *Acta Mater.*, 2008, **56**, 5845-5860.

[25] Braeck S. & Podladchikov, Y.Y., *Phys. Rev. Lett.*, 2007, **98**, 095504.

[26] John, T., Medvedev, S., Rüpke, L.H., Andersen, T.B., Podladchikov Y.Y. & Austrheim, H., *Nature Geoscience*, 2009, **2**, 137-140, doi:10.1038/ngeo419

[27] Brantut, N., Schubnel, A., Corvisier, J., & Sarout, J., *J. Geophys. Res.*, 2010, **B05314**, B05314.

[28] Brantut, N., Sulem, J., & Schubnel, A. (2011). *J. Geoph. Res.* 2011, **116**, B05304.

[29] Jamtveit, B., Austrheim, H. & Putnis, A., *Earth-Science Reviews*, 2016, **154**, 1-13.

[30] Tallakstad, K.T., Toussaint, R., Santucci, S., Schmittbuhl, J. & Måløy, K.J., *Phys. Rev. E*, 2011, **83**, 4, 046108. doi: 10.1103/PhysRevE.83.046108
[31] Tallakstad, K.T., Toussaint, R., Santucci, S., & Måløy, K.J., *Phys. Rev. Lett.*, 2013, **110**, 145501. doi: 10.1103/PhysRevLett.110.145501

[32] Toussaint, R. & Pride, S.R., *Phys. Rev. E* **66**, 2002, **3**, 036135. doi: 10.1103/Phys-RevE.66.036135

[33] Toussaint, R. & Pride, S.R., *Phys. Rev. E* **66** 2002, **3**, 036136. doi: 10.1103/Phys-RevE.66.036136

[34] Toussaint, R. & S.R. Pride, S.R., *Phys. Rev. E* **66**, 2002, **3**, 036137. doi: 10.1103/Phys-RevE.66.036137

[35] Måløy, K.J., Santucci, S., Schmittbuhl J. & Toussaint, R., *Phys. Rev. Lett.* **96**, 2006, **4**, 045501 (2006). doi: 045501 10.1103/PhysRevLett.96.045501

[36] Sethna, J.P., Dahmen, K.A. & Myers, C.R., *Nature*, 2001, **410** (6825), pp. 242-250. (2001). doi: 10.1038/35065675

[37] Houle, P.A. & Sethna, J.P., *Phys. Rev. E*, 1996, **54**, 278.

[38] Salminen, L.I., Tolvanen, A.I. & Alava, M.J., *Phys. Rev. Lett.*, 2002, **89**, 185503.

[39] Bradbury, R. Fahrenheit 451. Ballantine Books (1953). ISBN 978-0-7432-4722-1.

[40] Alava, M. & Niskanen, K., *Reports on Progress in Physics*, 2006, **69**, 3, 669-723.

[41] Stojanova, M., Santucci, S., Vanel L. & Ramos, O., *Phys. Rev. Lett.*, 2014, **112**, 115502.

[42] Tanaka, A. & Yamauchi, T., *J. Pack Sci. & Tech. Jpn.*, 1997, **6**, 5, 268-276.

[43] Tanaka, A., Otsuka, Y. & Yamauchi, T., *Tappi J.*, 1997, **80**, 5, 222-226.

[44] Yamauchi, T., Okumura, S. & Noguchi, M., *J. Mater. Sci.*, 1993, **28**, 4549.

[45] Yamauchi, T., 2012, Application of IR Thermography for Studying Deformation and Fracture of Paper, Infrared Thermography, Dr. Raghu V Prakash (Ed.), ISBN: 978-953-51-0242-7, InTech, doi: 10.5772/29961.

[46] Boler, F.M. & Spetzler, H., *Pure and Applied Geophysics*, 1986, **124**, 4-5, 759-772.

[47] Boler, F. M., *J. Geophys. Res.*, 1990, **95**, B3, 25932607, doi:10.1029/JB095iB03p02593.

[48] Gross, S.P., Fineberg, J., Marder, M., McCormick, W.D. & Swinney, H.L., *Phys. Rev. Lett.*, 1993, **71**, 3162.

[49] Mc Garr, A., *J. Geophys. Res.* **104**, 1999, **B2**, 3003-3011.

[50] Tantot, A., Santucci, S., Ramos, O., Deschanel, S., Verdier, M.-A., Mony, E., Wei, Y., Ciliberto, S., Vanel, L. & Di Stefano, P.C.F., *Phys. Rev. Lett.*, 2013, **111**, 15.

[51] De Groot, S. R. & Mazur, P., Non-Equilibrium Thermodynamics, Dover, Mineola NY, 1984.

[52] Welty, J.R. et al. Fundamentals of Momentum, Heat and Mass transfer, 5th ed., John Wiley and Sons. (2007) ISBN 978-0470128688

[53] Press, W.H. et al., Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge University Press, 2nd ed. (1992).

[54] Jessup, R. S. & Prosen, E. J., *Journal of Research of the National Bureau of Standards*, 1950, **44**, RP2086.