Impact Statement

The goal of this work is to provide a modern look at the Trinity test using tools that were unavailable to Taylor and other scientists at the time, and to reconduct the same experiments that led Taylor to his famous calculation of the energy yield. We also aim to provide a comprehensive background on the multiple methods that scientists, both modern and past, have used to calculate the energy yield of the Trinity test and provide a comparison between values.

2. Background

The Trinity explosion, shown in Figure 1, was the name for the first test of the atomic bomb that occurred on July 16, 1945 in New Mexico, just a month before the bombs were dropped into Hiroshima and Nagasaki. The test bomb released a huge amount of energy, causing radioactive fallout in the area and a blast that could be seen from space, reaching 40,000 feet in only 7 minutes (Szasz, Reference Szasz1995). The atomic bomb would go on to end the war and have devastating consequences for the civilians of Japan, raising profound questions on the very nature of our existence and the use of weapons of such cataclysmic destruction. Those involved in the Manhattan projectFootnote ¹, that is, the name of the team that developed the atomic bomb, as well as others such as Taylor, spent much time actualizing the idea of a fission bomb and understanding the effects of such a device.

Figure 1. Picture of the Trinity explosion (nps.gov, 2023).

There are multiple approaches one can take when calculating the energy of a blast, such as one from an atomic bomb. During the time of the Trinity test and World War II, there were three main minds working on the problem. Taylor, the focus of this paper, is often cited with the calculation as he was the only one who gave proof of the validity of a point source model, as well as calculating the energy within a range of error. Two others, Hungarian-American John von Neumann and Russian Leonid Sedov, also performed similar calculations at the time, and actually achieved better accuracy than Taylor. Taylor’s approach was to use partial differential forms of the equations for motion, continuity, and state as well as integral forms of energy to arrive at an approximate formula for the energy of the blast wave. Sedov and von Neumann followed similar approaches that began with dimensional analysis. The two then determined integral forms for energy and computed them analytically, although neither directly calculated energy. A more in-depth exploration of their calculations is provided in their own papers (Sedov, Reference Sedov1946; Neumann, Reference von Neumann1958) and in the work of Deakin (Reference Deakin2011b). All three did not find exact values as the radiative energy was not taken into account (Deakin, Reference Deakin2011b).

Table 1 shows the energy estimates of various authors, including Taylor’s estimate. Taylor’s, Deakin’s, Sedov’s, and von Neumann’s calculations are all under the assumption that the specific heat ratio $ \gamma $ is constant and equal to 1.4. Hanson’s calculations—done more recently—use physical samples from the explosion to calculate the amount of plutonium and other elemental constituents of the bomb, from which the yield is determined (Hanson et al., Reference Hanson, Pollington, Waidmann, Kinman, Wende, Miller, Berger, Oldham and Selby2016). President Truman and General Groves, who were involved in the Manhattan project, released the values presented in the table once the information was declassified (Deakin, Reference Deakin2011a,Reference Deakinb).Reference Deakin2011).

Table 1. Summary of Trinity test calculations

2.1. Notation

Table 2 contains the variables used both in this paper and by TaylorI (Reference TaylorI1950), TaylorII (Reference TaylorII1950), and will be used throughout this brief paper.

Table 2. Nomenclature

2.2. Dimensional analysis approach

Prior to detailing Taylor’s approach, it will be useful to touch upon dimensional analysis with the Buckingham Pi theorem. This theorem states that with $ N $ variables and $ K $ fundamental units there are $ P=\left(N-K\right) $ dimensionless quantities (Buckingham, Reference Buckingham1914). By using this theorem, we can identify the unit-less quantities that will enable us to calculate the energy of the explosion. To complete such an analysis, the three scientists, Taylor, Sedov, and von Neumann, reduced the explosion to starting from a point source and expanding outward in a spherical shock wave, rapidly releasing large amounts of energy. The use of a shock wave suggests that the state of the air can be expressed in two categories: internal shock wave properties and external conditions of the undisturbed atmosphere. Therefore, we can devise seven variables: $ E $ the energy, $ R $ the radius, $ t $ the time, $ p $ the static pressure of the wave, $ {p}_0 $ the pressure of the undisturbed atmospheric pressure, $ \rho $ the static density of the shock wave, and $ {\rho}_0 $ the undisturbed density of the atmospheric density (see Figure 2). All variable units, when broken down into their base units, contain three fundamental units: mass, length, and time. In terms of the Buckingham Pi theorem, there are $ N=7 $ variables and $ K=3 $ fundamental quantities; therefore, this problem contains $ P=4 $ dimensionless quantities (Deakin, Reference Deakin2011b).

Figure 2. Flow properties inside and outside a spherical shock front.

We wish to focus on a quantity that involves the target variable of energy, as well as easily measurable quantities, which in this case are radius, time, and the density of undisturbed air. To find the exact form of the dimensionless quantity, we can use dimensional analysis to solve for energy. In SI units, energy is recorded in Joules $ (J) $ which is equivalent to a Newton-meter $ \left(N\cdot m\right) $ or kilogram-meter-squared-per-second-squared $ \left( kg\cdot {m}^2/{s}^2\right) $ . Therefore, an equation for energy must contain units of mass $ (kg) $ , length $ (m) $ , and time $ (s) $ . Mass is not directly used within the context of the problem, instead density (units of $ kg/{m}^3 $ ) takes its place. In Taylor’s paper, he is interested in the time since the blast, the density of the undisturbed atmosphere, and the radius of the shock wave. Thus, the equation for energy has the form

(1) $$ E=K{\rho}_0^a{R}^b{t}^c $$

where $ a,b,c, $ and $ K $ are the constants (Díaz, Reference Díaz2021).

There are a few important points to note regarding the above. First, energy is per second squared, and since neither density nor radius contain units of time, the equation for energy must contain time, leading to the negative second power, or $ c=-2 $ . Second, density is the only variable that contains units of mass, and since both energy and density contain mass only to the first power, then $ a $ must equal 1. However, the inclusion of density adds an additional $ {m}^{-3} $ to the equation that must be corrected in order to get $ {m}^2 $ in the final form. Therefore, to fix this, we need a factor of $ {m}^5 $ , which can be provided by the radius, meaning $ b=5 $ and we have that

(2) $$ E=K{\rho}_0{R}^5{t}^{-2}, $$

where $ K $ is a constant and proven to be a function of $ \gamma $ by TaylorI (Reference TaylorI1950). This dimensional analysis is taken from Díaz (Reference Díaz2021).

2.3. Taylor’s first principles approach

While many believe that Taylor used the above approach (Deakin, Reference Deakin2011b), he instead used first principles and similarity assumptions to arrive at an equation for energy. He eventually arrived at an expression of the same form as equation (2), identifying the constant $ K $ . Taylor’s process and assumptions are captured in Figures 2–4; a detailed derivation is provided in the Data Availability Statement section.

Taylor provides the similarity assumptions in Figure 3 (1a-c) based on his earlier works in 1946 and 1950 on spherical detonation and expansion; he comes to the conclusion these are the appropriate forms for the expansion of a blast, proving that they obey continuity and motion as outlined above. In Figure 3, Taylor uses Navier–Stokes under the assumption of incompressible flow and no gravitational or viscous effects, resulting in the incompressible Euler equations. This process lead him to three equations for the derivatives of the non-dimensionless ratios $ f $ , $ \psi $ , and $ \phi $ of pressure, density, and velocity, respectively, where $ {f}^{\prime }= df/ d\eta $ , $ {\phi}^{\prime }= d\phi / d\eta $ , and $ {\psi}^{\prime }= d\psi / d\eta $ . These are expressed as

(3) $$ {f}^{\prime}\left({\left(\eta -\phi \right)}^2-\frac{f}{\gamma}\right)=f\left(-3\eta +\phi \left(3+\frac{1}{2}\gamma \right)-2\frac{{\gamma \phi}^2}{\eta}\right) $$

(4) $$ {\phi}^{\prime}\left(\eta -\phi \right)=\frac{1}{\gamma}\frac{f^{\prime }}{\psi }-\frac{3}{2}\phi . $$

(5) $$ \frac{\psi^{\prime }}{\psi }=\frac{\phi^{\prime }+\frac{2\phi }{\eta }}{\eta -\phi } $$

Figure 3. Introduction of nondimensional pressure, density, velocity, and first principles analysis.

Using equations (3)–(5), Taylor performs a discretized step (first-order) calculation for the values of $ f $ , $ \phi $ , and $ \psi $ using

(6) $$ {\displaystyle \begin{array}{l}f\left({\eta}_i\right)\approx f\left({\eta}_{i-1}\right)+{f}^{\prime}\left({\eta}_{i-1}\right)\cdot \left({\eta}_i-{\eta}_{i-1}\right),\\ {}\phi \left({\eta}_i\right)\approx \phi \left({\eta}_{i-1}\right)+{\phi}^{\prime}\left({\eta}_{i-1}\right)\cdot \left({\eta}_i-{\eta}_{i-1}\right),\\ {}\psi \left({\eta}_i\right)\approx f\left({\eta}_{i-1}\right)+{\psi}^{\prime}\left({\eta}_{i-1}\right)\cdot \left({\eta}_i-{\eta}_{i-1}\right),\end{array}} $$

where $ {\eta}_i-{\eta}_{i-1} $ represents a small decrement in the nondimensional radius. Taylor performs his step calculation over the nondimensional radius range $ \eta \in \left[\mathrm{0.5,1}\right] $ with intervals of $ {\eta}_i-{\eta}_{i-1}=-0.02 $ for the values of $ f $ , $ \phi $ , and $ \psi $ . He sets the initial conditions in equation (6) to be at $ \eta =1 $ , found by the Rankine–Hugoniot solutions and calculated in box (7) of Figure 4. In this case, $ f $ , $ \psi $ and $ \phi $ can be expressed solely as functions of $ \gamma $ (see top half of Figure 4). From this calculation and the processes outlined in Figures 3 and 4, Taylor finds approximate forms of $ f $ , $ \phi $ , and $ \psi $ . Note that while his exact procedure is difficult to follow, we performed a similar analysis, plotted as a black line (labeled step calculations) in Figure 5.

Figure 4. Approximate forms of nondimensional parameters.

Figure 5. f, ϕ, and ψ versus η_.

To solve for analytical forms of the three nondimensional quantities, Taylor substitutes equation (5) into (3), to arrive at

(7) $$ \frac{f^{\prime }}{f}\left(\eta -\phi \right)=\gamma {\phi}^{\prime }-3+\frac{2\gamma \phi}{\eta }. $$

However, this equation cannot be directly solved unless $ {f}^{\prime } $ is assumed small, such that the $ {f}^{\prime } $ term can be ignored. This assumption can be made at early values of $ \eta $ where $ \phi $ is approximately linear as $ f $ minimally increases in a lower $ \eta $ range, as shown in Figure 5, at $ \eta $ values of <0.6 and $ {f}^{\prime }<0.035 $ . Using this assumption, he evaluates the initial approximate form as $ \phi =\frac{\eta }{\gamma } $ . He notes this form, shown as the blue dashed line (labeled initial approximation) in Figure 5, does not agree well with the step calculations above $ \eta \approx 0.7 $ , where the $ {f}^{\prime } $ small assumption breaks down. He suggests the approximate form would not be linear—containing a term of order $ n $ . This form is shown in box (8) of Figure 4. Taylor calculates the values of the constants $ \alpha $ and $ n $ using the initial conditions set forth by the Rankine–Hugoniot solutions at $ \eta =1 $ . From the corrected approximate form of $ \phi $ , he used equations (3)–(5) to arrive at approximate forms for $ \psi $ and $ f $ as well (TaylorI, Reference TaylorI1950), as shown in step (8) of Figure 4. These corrected approximation forms are shown as red dashed lines (labeled approximate form) in Figure 5. One can observe that the step calculations (black) have relatively good agreement with the corrected approximate forms (red), with almost complete agreement for $ \psi $ . Once Taylor found these approximate forms, he was able to calculate the constant $ K $ in his energy equation.

In Figure 6, we show his process for formulating the equation for energy. This marks the end of what he accomplished in his first paper with regard to calculating the energy of the blast (TaylorI, Reference TaylorI1950), as values for radius and time are still needed. Once the Trinity test was declassified, he obtained those values, and in his second paper (TaylorII, Reference TaylorII1950), he linearly interpolates data on the log scale plot of $ {R}^{5/2} $ versus $ t $ and uses the intercept value for $ {R}^5{t}^{-2} $ in his energy calculations.

Figure 6. Energy calculations.

3. Revisiting Taylor’s analysis

In this paper, we take the same approach as Taylor but make use of modern techniques to improve upon his calculation. We started by reviewing his work, as documented in the previous section, to get a full understanding of his approach and assumptions.

To emulate his radius measurements, we made use of the open source Computer Visualization Annotation Tool (Sekachev et al., Reference Sekachev, Manovich, Zhiltsov, Zhavoronkov, Kalinin, Hoff, TOsmanov, Kruchinin, Zankevich, DmitriySidnev, Markelov, Johannes, Chenuet, a-andre telenachos, Melnikov, Kim, Ilouz, Glazov, Priya, Tehrani, Jeong, Skubriev, Yonekura, truong, Zliang, lizhming and Truong2020) to annotate the sets of photos by Mack, the U. S. Atomic Energy Commission (Reference Mack1946), and others presented by Taylor in his second paper (TaylorII, Reference TaylorII1950), though we could not find an image for $ t=3.26 $ ms. With this tool, we made an outline of the blast for earlier times as well as annotating two circles that represent the smallest enclosing circle and the largest enclosed circle for all times. Later times do not include a blast outline as the entire blast is not visible within the photographs. For later time stamps, an arrow is provided by Mack and U. S. Atomic Energy Commission (Reference Mack1946) that designates the shock front, and this is used as the maximum circle for these times. The minimum enclosing circle was checked using the Welzl method (Welzl, Reference Welzl and Maurer2005) and the maximum enclosed circle by testing circles centered at Voronoi points (Voronoi, Reference Voronoi1908) of the polygon created by the blast outline. An example of this method along with the radii that Taylor provided for $ t=0.36 $ and $ t=3.53 $ ms is shown in Figure 7a,b. The larger of the two circles is closer to the approximation that Taylor made, and is the more likely radius as the shock wave is the radius we are interested in, which occurs in front of the fire blast.

Figure 7. Estimating the blast radius from images taken at different times and plotting the data on a logarithmic scale.

The assumption that Taylor makes based on his first principles analysis, shown in Figure 3 step (3a), is that the radius and time are related by $ At=\frac{2}{5}{R}^{-5/2} $ where A is a constant. Therefore, the radius and time of the blast should correspond to:

(8) $$ \frac{5}{2}\log (R)=\log (t)+\log (n) $$

where $ n $ is a constant and $ {n}^2 $ corresponds to the value of $ {R}^5{t}^{-2} $ . Using the radii found in the previous paragraph, we plot the logarithmic-scale graph of $ {R}^{5/2} $ versus $ t $ as Taylor did in his second paper (TaylorII, Reference TaylorII1950); this is shown in Figure 7c. This figure contains three lines with a slope of one whose y-intercepts correspond to values of $ {R}^5{t}^{-2} $ for the minimum, maximum, and average radii based on equation (8). These values are calculated using standard linear least-squares interpolation to fit a linear approximation. Figure 7 also displays a zoomed-in view of the calculations from $ t\in \left[\mathrm{3.53,4.61}\right] $ ms as the difference between our values and Taylor’s lead to a large discrepancy in the comparison of energy values shown in the next section.

To calculate $ K $ in the energy equation shown in Figure 6, we used the equations for $ \phi $ , $ \psi $ , and $ f $ put forth by Taylor in his first paper (TaylorI, Reference TaylorI1950) and shown in Figure 3 step (8). These integrals do not have an analytical solution, so the open-source code equadratures (Seshadri et al., Reference Seshadri, Wong, Scillitoe, Ubald, Hill, Virdis and Ghisu2021) was used to computationally estimate these integrals, that is,

(9) $$ K\approx \frac{4{\omega}_j}{25}\left(2\pi \sum \limits_{j=1}^{J=20}\psi \left({\eta}_j\right){\phi}^2\left({\eta}_j\right){\eta}_j^2+\frac{4\pi }{\gamma \left(\gamma -1\right)}\sum \limits_{j=1}^{J=20}f\left({\eta}_j\right){\eta}_j^2\right) $$

where $ {\left\{{\eta}_j,{\omega}_j\right\}}_{j=1}^{J=21} $ are Gauss–Legendre quadrature points and weights of order 20 across the support $ \eta \in \left[0,1\right] $ . This order was found to be sufficiently high for numerical precision (to nine decimals), a convergence analysis can be found in the GitHub link. We calculated $ K $ using (9) for 10 values of $ \gamma $ in the range $ \left[\mathrm{1.2,1.667}\right] $ . The values calculated using this method show minor differences from those Taylor calculated; see Table 3. It is difficult to ascertain what approach Taylor took in calculating the integrals, as he simply states in TaylorI (Reference TaylorI1950) that he evaluates the integrals and uses a step calculation. Based on the differences in the K calculation alone (not considering differences in radius/time), this would put our energy ranges within a 1.89% difference from Taylor.

Table 3. K Values computed using (9) contrasted with those reported by Taylor

Using the values for $ {R}^5{t}^{-2} $ and K calculated in (8) and (9), respectively, the total energy was estimated using the same equation as Taylor, shown in Figure 6 step (9). This was done for both the minimum and maximum radii. These values are contrasted to those that Taylor presents in his paper (TaylorII, Reference TaylorII1950).

6. Lessons learned

• • While dimensional analysis provides a simple process for finding the structure of a driving equation, it cannot account for constant values as Taylor’s first principles approach can.

• • While not exact, using reliable assumptions (like those pursued by Taylor) can provide one with an easier process to calculate a value within a margin of error.

• • Revisiting the analysis with modern tools assisted in providing a value range, more accurately representing the photographic data, and quantifying the uncertainty, an aspect Taylor failed to include.

• • This is very useful in situations where a back-of-the-envelope-type calculation is required, and where perhaps expensive simulations are not feasible or necessary.

• • Attempting a recreation of historical works and accessing the same references is difficult with the differences in technology and an inability to contact the author.