Dissociation Fit Test

  The procedure used to fit the computed dissociation data is very sensitive to the value of A_eq so we should do a more rigorous test by adding some random errors to the data. The dissociation products of the diatomic molecule do not have to be the same since the dissociation produces two molecules and the resulting equations end up the same.

An estimate of the equilibrium constant K'_eq can be found from the ratio of the products after equilibrium has been reached and averaging eliminates most of the measurement error.

As mentioned before the first reaction rate constant dominates at the start of the dissociation since there are no products. We can estimate it by doing a linear fit of a few of the initial data points. The function F(x) is arbitrary and was chosen to give a good fit for the data. With K'_eq and k'₀ known we can compute the second reaction rate constant k'₁.

Like K'_eq, and estimate of A'_eq can be found by averaging.

We can then compute D' and λ' with the formulas above from the solution of the first order differential equation and solve for κ' using the first two data points and the formula for A.

The result gives a fairly good fit for the data.

There is an unused relation between A₀, A_eq, K_eq and D so the one might be able to adjust the values of K' , A'_eq and D' for a slightly better fit by assuming that A₀ is correct.

Supplemental (Oct 31): Since the value of the first data point, A₀, is most likely well known one gets a better estimate of κ' by solving
2(A₀ - A'_eq) = D'(1 - tanh(κ')) for κ'.

Temperature Dependence of the Dissociation Rates and the 'Time of Expectation'

 In 1885 J. J. Hood published a paper On the Influence of Heat on the Rate of Chemical Change in Philosophical Magazine in which he gave the empirical fit by an exponential function for some data on the reaction rate for the oxidation of ferrous sulfate (FeSO₄) by potassic chlorate (KClO₃) at various temperatures. The justification for the use of an exponential function is that the data appears to be almost linear in a semi-log plot.

Later in 1889 Svante Arrhenius published Über die Reaktionsgeschwindigkeit bei der Inversion of Rohrzucker durch Säuren (On the Reaction Rate in the Inversion of Cane Sugar by Acids) in Zeitschrift für Physicalische Chemie which gave a better fit for Hood's data based on some thermodynamics considerations. In the table on the page of the link "beob." is the observed rate of the reaction, "ber.₁" is the calculated rate using formula (1) in the paper and "ber.₂" gives the calculated rates for Hood's empirical formula for comparison.

Formula (1) in this paper is equivalent to the statement that the reaction rate ρ divided by e^-A/T is equal to a constant which is essentially the Arrhenius equation. This is the exponential that is found in Schrödinger's equation for the 'time of expectation'. The Arrhenius constant A is now written as activation energy divided by the gas constant R, the constant from the ideal gas law, PV = nRT. The k in Schrödinger's equation is the corresponding constant for a molecule. His 'time of expectation' is inversely proportional to the reaction rate.

Curve Fit for the Dissociation Rate Constants

  If one looks at the equation for the amount of the dissociated molecule as a function of time in the last blog one one can see that it can be rewritten as a linear function of A equal to the hyperbolic tangent of a linear function on the time t. Not all the constant coefficients in this equation are independent but we can find the first two as a function of the third by evaluating the equation at t = 0 and t = ∞ since the initial amount A₀ and equilibrium amount A_eq are known. Substituting the functions for μ and ν and separating out the coefficient of tanh(κ) we get two new functions of A, α and β, the sum of which is 1.

So we are left with two constants to be determined since α, β and t can be obtained from the experimental data. If κ is known we can compute λ so we can estimate the function λ(κ). One then has to search for the value of κ for which the square of the difference between the two sides of the equation is a minimum.

I did this for some of the data for A in the last blog and found a minimum for the variance V(κ) when κ = 0.458. With κ known we can go back and compute μ, ν and λ and consequently the equilibrium constant K and the dissociation rate constants k₀ and k₁.

The computed values for the fit are almost identical to the values chosen in the last blog since we used values for A that were nearly exact and so there was negligible error. The round-off errors and the estimate for λ(κ) contributed some error to the results. As we can see the variance changes very slowly near the optimum value for κ.

edit (Oct 26): Corrected the formula for λ(κ) in the derivation by removing the 2nd atanh that crept in with a paste. The correct formula was used for the fit so fortunately the error in the derivation didn't affect the results. The number of data points determines the accuracy of the fit but one gets a reasonably fit with just 20 data points. The calculation above used 100 data points and 200 came closest to the original values. Measurement errors also would affect the fit and consequently an error in κ will propagate through to the dissociation constants.

The Experimental Determination of Reaction Rates

  By the end of the 19th century chemists had succeeded in formulating rules for determining the rates of reactions and their dependence on temperature. Arrhenius discussed the dissociation of a substance into its component atoms. As an example of how one can determine the rates of a reaction and their dependence on temperature we can study the simple case where a diatomic molecule is split up into its constituent atoms.

The first two highlighted equations give the rates for the formation of the molecule A and its atoms B. There are two ways that we can change the amount of A. The first term in the first equation is the rate at which it spontaneously breaks up into two atoms of B. It is negative and proportional to A because every molecule of A has an equal chance undergoing the dissociation. The second term takes gives the rate at which two atoms of B come together to form A. The rate is positive and proportional to B² since each atom of B can potentially interact with all the other atoms of B. The two equations are difficult to solve simultaneously but we can use the fact that the total number of atoms is a constant. The result is a nonlinear rate law for A which can be solved using a standard integral formula. To simplify the calculation of A we can define some additional constants. The system converges to an equilibrium value A_eq and an equilibrium constant K can be defined that relates the amounts of A and B. The formula for A_eq is found by letting
t → ∞.

The solid curves in the figure below show the changes in the amount of A and B with time. Since initially B = 0 and A follows the rate law for simple decay with k₀ playing the role of the decay constant. By determining the initial slope of the curve for A and the equilibrium constant we can experimentally determine the reaction rate constants, k_i.

By repeating this experiment at different temperatures we can determine the dependence of the rate constants on temperature.

Supplemental (Oct 24): One can find papers by Deville, Wilhelmy, van't Hoff, Guldberg & Waage, Arrhenius & Oswald in Leicester, A Source Book in Chemistry, 1400 - 1900. Ludwig Wilhelmy studied the dissociation of sucrose into its simpler sugars and gave the rate equation and the formula for exponential decay in the amount of sucrose. Arrhenius won the Nobel Prize in 1903 for his work on conduction in solutions and is known for studying the influence of atmospheric CO₂ on ground temperature.

The "Life Expectancy" of a Molecule

  Schrödinger's "time of expectation" in What is Life? could use a little clarification. It is the expected amount of time for some change or transition to take place in a molecule at a given temperature. Suppose the transition in question is the time it takes to inactivate a virus. Using probability theory we can assign a probability q to the chance that the molecule will make the transition in a given time interval Δt and a probability p = 1 - q that it does not. The transition state diagram is shown in the figure below.

The molecule starts in the upper left state and it moves one step to the down or right depending on whether on not it makes the transition. The expected number of intervals that the molecule will survive is determined as follows from the transition matrix.

The powers of the transition matrix can be determined by doing repeated multiplications and noting that the right column doesn't change and that the top left component is just p^k. To get the last component we use the fact that the sum of each column is 1 for a probability matrix. The expected value of k is just the sum of k times the probability that the molecule has survived k times. The sum can be shown to be equal to p/q². The expected time of survival is the time interval multiplied by the expected number of survivals. We can check this result by comparing the two functions above for the expectation.

   The result is not a linear function of the probability p since q = 1 - p. For very low rates of survival however the survival time is approximately a linear function of the probability of surviving an interval of time.

  We are now in a position to modify Schrödinger's formula for the time of expectation. First we need to rewrite p as an exponential function of time by noting k = t/Δt to get the decay rate λ. This gives the same curve that we got for the finite steps. Substituting for p and q gives the correction for the time of expectation.

For chemical reactions the decay rate λ is temperature dependent as Schrödinger indicated.

Supplemental (Oct 22):  see also:  Arrhenius equation   van 't Hoff equation

A Precessing Gytoscope

  The motion of a gyroscope can be described by a double rotation. One can get one to precess about its center of mass by placing it on a glass tabletop so that the bottom of the spin axis skates across the surface. This gyroscope is spun up by a piece of string while holding the outer frame. It might help to do this at an angle and give a little twist as one sets it down to get it to precess about the center.

An Animation of a Permutation Rotation

  To find a rotation that would take the axes x → y, y → z and z → x two generators are required and the calculation proceeds as follows.

To find out what happens if we let the change occur incrementally can allow the two angles to change from 0° to 360 ° at a steady rate and compute a set of X, Y and Z values that indicate the position of the tips of the axes unit vectors.

We can view the results as an animated video. The red, green and blue colors indicate the x-, y- and z-axes respectively. Initially the rotation is difficult to follow but if one focuses on the blue z-axis one sees that it moves in a circle in the xz-plane while in the rotating x'y'-plane the x-axis always moves towards the y-axis.

Supplemental (Oct 6): I was able to upload a higher resolution video clip to YouTube and then add it to this post. Here's the link to the YouTube video. It is best viewed there in Theater mode. The color code is (R, G, B) represents (x, y, z). The same double rotation is generated if the second rotation is replaced by z → x.

The Equivalent of a Series of Rotations

  One can subject the coordinate axes to a series of plane rotations and then try to find an equivalent rotation to produce the same orientation of the axes. Suppose one first moves the x-axis in the xy-plane leaving the z-axis fixed then moves the position of the new x-axis towards the z-axis leaving the direction of the new y-axis fixed. We can follow this with another rotation about the fixed x-axis. The rotation generator formula allows us to find the required rotation matrices for each step and multiplying these in the correct order gives the equivalent rotation. The generator G(a,b) allows us to define a vector function Γ(a,b) which produces a set of matrix coefficients N, P and G for the plane rotations. These are saved in the matrix X and multiplied by the trig functions to give the R_k(θ). The order of the plane rotations in the formula for the equivalent rotation is from right to left since they act on vectors and matrices on their right side.

The matrix coefficients don't just depend on the unit vectors that define them. As this example shows the plane of rotation can depend on a previous rotation so x'=x'(θ₀).

nD Rotation Summary

  This is a good point to collect some of the results. Note that the unit vectors a and b can be thought of as a radial unit vector and its tangent that define a point in the plane and the direction of its rotation. I've added the effects of operating on an arbitrary vector v with matrices in the set {N, P, G}. The vector n is assumed to be normal to the plane. Pv projects v onto x in the plane and Gv projects the same magnitude onto y which is also in the plane and perpendicular to x in the direction of rotation.

Multiplicative Inverses

  A matrix of the form Z = aN + bP + cG has an inverse Z' if Z satisfies certain conditions.

The inverse of a rotation R = N + α P + β G is its transpose R^T.

edit (Oct 2 1am): Z must satisfy conditions for inverse

Using Multiplication Tables to Compute Powers of R

  If one looks at the multiplication tables for the sets of matrices {I, P, G} and {N, P, G} one sees that the second set has more zeros which suggests that it be better for calculating repeated rotations of R. This is like finding the powers of a complex number with G playing the role of i. The multiplication tables are independent of G and consequently so are the formulas which use the variables α and β, the direction cosines in the plane.

Supplemental (Oct 1): The multiplication table for the set {I, P, G} indicates that the algebraic structure is a ring rather that a group since there is no multiplicative inverse. It acts a lot like numbers of the form a + b√q where q is some integer.

Simple Rotations are Matrix Exponential Functions

  One can rewrite the formula for the rotation R in terms of an exponential function. First let us recall that  for vectors in the plane the equation for this formula is,

R(θ) = I cos(θ) + G sin(θ).

We can rewrite any vector v in the entire space as the sum of a vector in the plane and another perpendicular to it. So v = Pv + (I - P)v where P is a matrix that projects v onto the ab-plane. If we calculate the square of G = ba^T - ab^T we get

G² = (ba^T - ab^T)(ba^T - ab^T) = b0a^T - b1b^T - a1a^T + a0b^T = -(aa^T + bb^T)

The part in brackets at the end is a matrix that will project any vector onto the ab-plane so we can let P = -G² and the rotation formula for any vector will be,

R(θ) = I - P + [I cos(θ) + G sin(θ)]P = I - P + P cos(θ) + G sin(θ)

since GP = G. If we go to the Wikipedia article on matrix exponentials we will find that the expression on the right is just e^Gθ. So we conclude that for simple rotations,

R(θ) = e^Gθ =  I - P + P cos(θ) + G sin(θ).

If we note that,

e^GθP =  P cos(θ) + G sin(θ) then,

R(θ) = I - P +e^GθP = I - (I - e^Gθ)P,

the formula in the last blog, which works fairly well but isn't as elegant as R(θ) = e^Gθ.