1 Basics, Initial and Boundary Value problems, Steady state, Phase portrait, Chemical kinetics
Contents
1 Basics, Initial and Boundary Value problems, Steady state, Phase portrait, Chemical kinetics#
# import all python add-ons etc that will be needed later on
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
init_printing() # allows printing of SymPy results in typeset maths format
plt.rcParams.update({'font.size': 14}) # set font size for plots
1.1 Types of equations#
Differential equations are characterized into ordinary and partial ones. Ordinary differential equations only involve one independent variable and so have ordinary (total) derivatives; for example, the equation for a damped harmonic oscillator is
with time \(t\) as the independent variable, \(m\), \(a\), and \(b\) are constants, and \(y\) is the displacement at time \(t\). Partial differential equations, such as the diffusion equation
have two or more independent derivatives, \(t\) and \(x\).
Differential equations are solved by integration and as with any integration, constants of integration are produced. These constants determine the exact solution of a differential equation and are determined by the initial values and/or the boundary values where \(y\) is given at some \(x\) and \(t\) values. Thus it essential to determine the initial / boundary conditions when solving any differential equation. These initial and boundary conditions are determined by the physics/chemistry of the problem being studied, for example initial concentrations in a reaction or initial position of a pendulum.
There are different notations used in differential equation \(dy/dx \equiv y'\) or \(y'(x)\). The second derivative is written as \(y''\) and so on. An alternative notation is \(\dot y\) or \( \ddot y \) etc., especially if the derivative is with respect to time.
1.2 Initial Value Problems (IVP)#
In these problems, the starting conditions only are specified. For example, with the equation \(dy/dt = ay\), the initial condition (or initial value) could be \(y(t_0)\) = 2 where \(t_0\) is the initial time, which is often zero. There is one initial condition because only a single integration step is necessary to solve the equation; the constant of integration is satisfied if \(y\) is specified at one value of \(t\). The initial value of \(y\) and \(t_0\) depends entirely on the problem being studied and is the equivalent of specifying the constant of integration in a normal integral.
If the equation is of second or higher order, two, three, or more initial conditions apply. One condition is needed for each stage of integration. Suppose the equation
has the particular solution
i.e we know from the physics/chemistry we are studying that these are the initial values. Our aim is to follow what happens at later times by solving the differential equation.
The second condition describes the gradient evaluated at time \(t_0\) and what this gradient is depends on the problem being analysed. You can appreciate these initial conditions if you play football, rugby, or golf; to get the ball to where you want it go to, the initial direction and speed (\(dy/dt\)) with which it is hit clearly matter.
The first stage of integration leaves us with one constant of integration and an equation in \(dy/dt\) and therefore \(dy/dt\) is needed at \(t_0\) to find this constant. Integrating this equation produces a second constant; hence, two initial conditions are needed. In any initial value problem once the calculation is started, there is no telling what value \(y\) will have since only the initial value has been fixed. Because many potential initial conditions could apply, all trajectories could start at \(t_0\), and have different \(y(t_0)\) and gradients or start at different \(t_0\) with the same gradient and so forth. These different starting conditions generate a field or swarm of solutions, or trajectories, Fig. 1. Solutions are sometimes represented on a graph by sets of arrows to indicate trends at different places and these are drawn as well as a particular solution.
Fig. 1 Different initial conditions lead to different trajectories from the same equation.
1.3 Boundary Value Problems (BVP)#
The boundary value problem defines conditions at two places that have to be satisfied throughout the calculation, thereby constraining the solutions. The most familiar example is probably the Schroedinger equation describing the quantum mechanical particle in an infinitely high box. This is solved with boundary conditions that the wavefunction must always be zero at both sides of the box, because if this were not the case an infinite amount of energy would be needed.
The Schroedinger equation equates the kinetic and potential energy, \(V(x)\), with the total energy but as the potential is zero inside the box only the kinetic energy remains and
The boundary conditions are \(\psi(x_0) = 0\) and \(\psi(x_L) = 0\) where \(x_0\) and \(x_L\) define the sides of the box.
1.4 Separable variables. Chemical kinetics 1.#
The rate of a chemical reaction can always be written down because this is proportional to the rate of loss or gain of a molecule’s population. For example, in a first-order process \(A \to B\), such as a cis-trans isomerisation or decay of a radioactive nucleus or of an electronically excited state, the rate of change of species A is proportional to the amount of A left unreacted or
where \(A_0\) is the amount initially present. The constant of proportionality \(k_1\) is the rate constant; more properly called the rate coefficient because in a chemical reaction, it depends on temperature. The negative sign is present because A decays into something else, and \([A]\) is the amount of A unreacted at time \(t\). It is important to remember that [\(A\)] is changing with time, although this is never explicitly stated in the equations; i.e. we do not normally write \([A]_t\). If the equation were written
with a positive sign, then it could, for instance, describe the rate of growth of the number of bacteria A in the presence of an unlimited food supply. ( Note that, as the notation \([A]\) or \([B]\) can be rather cumbersome, \(A\) and \(B\) will often be used instead to represent concentrations.)
First-order equations are derived by considering the difference in the number of molecules \(N\) present at a time \(t\) and the number reacted during a small time interval \(\delta t\). This can be expressed as, number unreacted at time
which can be written as
from which
and in the limit when \(\delta t\) becomes small, this equation becomes
This becomes equation (1) if the number of molecules is converted into a concentration and can be solved provided the amount of \([A]\) present initially is defined; the initial condition is \([A]_0 = A_0\) where \(A_0\) is a constant. Integration from \(t =0 \to t \) produces
as shown below.
If the reaction scheme is more complex, A \(\to\) B \(\to\) C, then the rate of decay of A is the same as just described, but B changes also and does so as
which is the rate with which A converts to B, less the rate that B reacts to form C. To solve this scheme, the amount of \([A]\) at time \(t\) is inserted into this equation. Reactions with many steps are dealt with in a similar manner, but eventually the scheme can become so complicated that only a numerical solution to the equations is possible.
In passing recall that if a reaction has the balanced equation of the general form,
which does not tell us what the mechanism actually is, then the rate of change of each species is
The first-order equation is found to describe many physical phenomena and is also representative of many types of equations where the variables are separable. This means that the equation can be written with \(y\) on one side and \(x\) on the other, or, for a first-order chemical reaction A on one side and \(t\) on the other. The general form of the separable equation is
which integrates to
which is the general solution since it contains an arbitrary constant \(c\). The particular solution is that obtained when either the initial or the boundary conditions are used.
In the first-order reaction of species A, the steps in the integration are
and the result should be familiar as the integration is a standard one. If the initial condition is that \([A]_0 = A_0\) at \(t = 0\) then the solution is
and indicates more clearly how the concentration varies with time. If the initial amount is \(A_0\) at time \(t = t_0\), then the equation is
As species A decays it must form another B, but if this does not decay then B can be obtained from the initial conditions because the total number of molecules must be constant; therefore
and B rises at the same rate as A falls. If, however, B decays to C, the scheme being A \(\to\) B \(\to\) C, then there is another rate equation
and we can substitute A into this to obtain
but now the terms are not separable, and another method is needed, which is described in Section 4. Sets of equations for sequential and parallel reactions such as \(dA/dt = \cdots\), and \(dB/dt = \cdots\) can alternatively be solved by eigenvalue methods described in Chapter 7 provided no terms with a product or ratio of concentrations is present, i.e. no A \(\times\) B terms, unless pseudo first-order conditions apply.
2 Steady State vs. Equilibrium#
In many chemical schemes, after the reaction has started it enters a period where intermediate or transient species can be identified, and their rate of change is effectively zero. This does not mean, however, that their concentration has to be small but if the rate of change is zero, then most rate equations can be solved relatively easily. This is of great utility because the complete solution can be very complex and often only numerical solutions are available. An example of this situation is the scheme \(A\to B\to C\) where the rate of change species B may be set to steady state but only once the reaction has started.
In other situations an equilibrium can be established between two species from which a second and slow reaction may proceed reaction, an example of this is
where the rate constant \(B\to C\) is small compared to those establishing the equilibrium.
In the scheme
the set of rate equations could be integrated, but if the rate of change of B is small then we can use the steady state approach,
making \(\displaystyle B = \frac{k_1A}{k_{-1}+k_2}\), and
If the equilibrium had been established first, with \(k_{-1} \gg k_2\) writing the equilibrium constant as
gives
which means that both the equilibrium is established before the rate determining step and that \(k_{-1}\gg k_2\).
2.1 Iodine recombination, I + M reaction#
The scheme
describes how iodine atoms are recombined to form I\(_2\) in the presence of an inert buffer gas \(M\). The transient collision complex \(IM\) can be supposed to exist at steady state once the reaction has started then,
making \(\displaystyle [IM]=\frac{k_1[I][M]}{k_{-1}+k_2[I]}\).
At steady state the rate of formation of I\(_2\) molecules is
The recombination of iodine atoms has be studied experimentally and the presence of IM species is confirmed by using flash photolysis, now called the pump-probe method. See G. Porter & J. Smith. Proc. Roy. Soc. (Lond.) v26,p142, 1961.
2.2 Decomposition of \(\mathrm{BH_3CO}\)#
In the gas phase decomposition of borane carbonyl the equilibrium
is rapidly established and the subsequent reaction
is so slow that the equilibrium is effectively maintained. Steady state can be applied to the \(\mathrm{BH_3}\) and then the fact that \(k_2\gg k_3\) is used to simplify the equation. For clarity we let \(BH_3\equiv X,BH_3CO\equiv XCO\). The rate equation for \(\mathrm{BH_3}\) is
At steady state \(d[X]/dt=0\) and therefore \(\displaystyle [X]=\frac{k_1[XCO]}{k_3[XCO]+k_2[CO]}\).
As
substituting for \([X]\) produces, after some re-arranging and simplifying,
which, by the stoichiometry, is also \(-d[XCO]/dt\). We are given that \(k_3\ll k_2\) and therefore
where \(K_e=k_1/k_2\) is the equilibrium constant. This equation can be solved by recognising that, by stoichiometry, if an amount \(x\) of CO is produced then \(a_0-x\) of \(\mathrm{BH_3CO}\) remains when \(a_0\) is the initial amount;
and by substituting \((a_0-x)=z\), \(-dx=dz\) and integrating produces the solution
where \(x \le a_0\). It is not obvious from this equation how the concentration of CO increases with time because \(x\) cannot be isolated on one side of the equation. However, a plot can be made by writing the equation as \(t=\cdots\) and put \(t\) on the \(x\)-axis. The \(\mathrm{[CO]}\) increases initially with a very large gradient until \(x/a_0 \sim 0.8\), in a plot of \(x/a_0\) vs. \(t/2k_3K_e\), and then the gradient rapidly reduces towards zero, i.e. the plot becomes almost horizontal, such that at long times \(x/a_0 \to 1\). At very early times expanding both terms in \(x/a_0\) gives \(x/a_0\sim t^2\).
2.3 Fluorescence yield#
The fluorescence yield of a molecule \(\varphi\) is the ratio of molecules excited to those that emit a photon and is defined as
If G is the ground state, S\(_1\) the excited singlet and T the triplet state, the scheme is
In this scheme \(k_a\) is the rate constant for absorption of a photon by the ground state, \(k_f\) the rate constant for fluorescence and \(k_s\) the intersystem crossing rate constant forming the triplet state. The S\(_1\) rate equation is
where \(k_aG\) is the rate of absorption. At steady state \(dS_1/dt = 0\), making \(\displaystyle S_1=\frac{k_aG}{k_f+k_S}\). The fluorescence yield is therefore
and is the fraction of molecules that fluoresce. A molecule fluoresces with a rate constant that is the sum of all process destroying the excited state; therefore, if \(k = k_f + k_S\) and as the fluorescence lifetime is the reciprocal of \(k\) or \(\tau = 1/k\), the yield becomes \(\varphi = k_f \tau\).
2.4 Quenching an excited state. Electron spin can affect the outcome#
The presence of paramagnetic species, such as O\(_2\), or heavy atoms or ions, Xe, I\(^-\) dissolved in solution will quench the excited state. The quenching process are energy transfer and spin orbit coupling. Quenching here means shortening the lifetime of the excited states, both singlet and triplet. Using nanosecond pump-probe (flash photolysis), which first became available shortly after the discovery of Q-switching in lasers in the 1960’s, it was discovered that the singlet excited state of several types of aromatic molecules was quenched at diffusion controlled rate constant but that their triplet excited states were quenched at approximately one ninth of this. The values for anthracene quenched with molecular oxygen are, for the singlet, \(k_q = 3.1\cdot 10 ^{10}\) and for the triplet, \(3.4\cdot 10^9\;\mathrm{dm^3 mol^{-1} s^{-1}}\).
To understand this behaviour a model of the quenching is made by supposing that a short lived collision complex is formed between the excited molecule \(M^*\) and the quencher \(Q\), i.e. \((M^*\cdots Q)\). The scheme is quite simple and is
where \(k_d\) is the diffusion controlled rate constant for the approach of the molecules and \(k_{-d}\) and \(k_R\) the first-order rate constants for the break up of the complex to reactants and products respectively. The rate of change of the excited complex is
We could include the decay of \(M^*\) by fluorescence or phosphorescence and solve the equation but it is necessary only to consider this at steady state in which case
The overall rate of reaction is by definition, \(k_q[M^*][Q]\), where \(k_q\) is the experimentally measured rate constant, thus
and therefore by substituting
(i) Reaction or Activation Control#
When the reaction is activation controlled this means that the reaction rate constant \(k_R\) is small thus the reaction’s activation energy is high and many, possibly millions of collisions occur before reaction takes place \(k_R \ll k_{-d}\) and
(ii) Diffusion Control#
When the reaction is diffusion controlled \(k_q=k_d\) and means that the reaction rate constant \(k_R\) is so large that reaction occurs essentially on first contact of the two species and then \(k_R\gg k_{-d}\).
The diffusion controlled rate constant has been calculated by Smoluchowski and combined with the Stokes-Einstein equation is
where \(\eta\) is the viscosity of the solution. This rate constant can vary widely, for example in ethanol \(k_d = 5.0\cdot 10^9\; \mathrm{dm^3 mol^{-1} s^{-1}}\) where the viscosity \(\eta = 1.2\) cP (cP = centiPoise = \(10^{-3}\) Pa s), but in glycerol \(k_d = 4.4\cdot 10^6 \mathrm{dm^3 mol^{-1} s^{-1}}\) because \(\eta = 14900\) cP.
(iii) Spin complexes#
The difference in quenching rate constant with molecular oxygen between singlet and triplet excited states can only be understood if the spin properties of quencher and excited state are examined. The oxygen molecule is paramagnetic with a triplet ground state, the term symbol is \(^3\Sigma_g^-\). There are two low lying excited states \(^1\Delta\) at \(8000\;\mathrm{cm^{-1}}\) (\(\approx 1250\) nm) and \(^1\Sigma_g^+\) at \(13000\;\mathrm{cm^{-1}}\) (\(\approx 770\) nm). Both the singlet and triplet excited state energies of the excited aromatic molecule are at least \(10000\;\mathrm{cm^{-1}}\) above that of the oxygen \(^1\Sigma\).
The spin angular momentum has to be conserved in any transition thus there are three pathways for quenching of the two triplet states. The initial states are the excited triplet state \(T\) and O\(_2\) ground state. The energy of the aromatic’s excited state must also be greater than any of the quencher states. (In the scheme \(^1M_0\) is the aromatic molecule’s ground state).
The superscript number \(1,3, 5\) is the spin multiplicity of each type of complex. There are four electrons in the complex so the total spin is \(2\) but according to addition rules for quantum spin the other allowed values values are \(2,1,0\) via a Clebsh-Gordon series. The total number of types of complexes is given by the total multiplicity which is \(9\). If \(S\) is the total spin quantum number which is \(0,1,2\) for singlet, triplet and quartet complexes respectively the multiplicity of each is \(2S+1\) making a total of \(9\) types.
Because the quintet complex can form no product as the triplet is an aromatic, the quenching rate constant is limited to \(4/9\) of its potential maximum, and has to be the sum of that due to energy transfer and intersystem crossing and this occurs via spin orbit coupling. The rate constant is
Experimentally as \(k_q=k_d/9\) this implies that \(k_e\gg k_{-d} \gg k_{isc}\), and so the experimental data is understood by considering the spin quantum numbers. Recall that the spin quantum number of an electron is \(1/2\) with projection or azimuthal or magnetic quantum numbers \(\pm 1/2\).
Notice that this mechanism indicates that triplet-triplet annihilation is a spin allowed process provided the energy is sufficient, i.e. there must be enough energy to produce the singlet product, which means that energy conservation must always apply. Note also that when a singlet excited state is quenched by a triplet there is only one type of complex formed because the singlet has total spin of zero and multiplicity of \(1\) which means that quenching always occurs at its full value.
2.5 The Phase Portrait#
The steady state is found when the gradient is zero, i.e. \(dy/dx = 0\). In non-linear differential equations which have higher powers of \(x\), such as \(dy/dx = x + ax^2 + bx^3\), there are many ‘steady states’, which are also called critical, fixed, or equilibrium points, and it is now useful to see where these are by plotting \(dy/dx\) vs. \(x\). This graph is called the phase portrait.
The new feature is that some equilibrium points are stable and others not. If an initial value of \(x\) is chosen close to a stable equilibrium point, conditions will change as \(x\) changes so that this point will be reached; if it is an unstable point, conditions now change so that this point is avoided. This may or may not mean that another stable point will be found,see Figure 2, which shows a schematic phase portrait. Some steady state points therefore can also be called equilibrium points because once the system reaches one of these it will remain there.
In Fig. 2 the arrows indicate the direction \(x\) will follow depending on its initial value; for example if 3 \(\le\) x \(\lt\) 2, then \(x\) will move towards, and end at, the stable point \(x\) = 3. The rule of thumb is that if the curve is above the \(x\)-axis the arrows move to the right, and they move to the left if the curve is below the line. A point starting in the range \(0\to 1\) or \(1\to 2\) moves to the stable point 1, and point 2 is unstable. Starting in the range \(\gt 2\), moves to stable point \(3\), but starting at zero or any negative value is unstable. A negative value of \(d^2y/dx^2\) at a steady state point also indicates stability.
Fig. 2 The stable ( filled circles ) and unstable equilibrium points in a schematic of a phase portrait.
If a problem is described by two coupled non-linear differential equations \(dy/dx = f (x,y)\) and \(dz/dx = g(x,y)\) then the phase portrait is not usually plotted but instead the phase plane in which is \(y\) vs \(z\) is plotted. Chapter 11 gives examples of using the phase plane. Strogatz (1994) discusses phase portraits and phase planes in detail and illustrates these with many interesting examples.
3 Phase planes and solving equations by separating variables#
3.1 Laser Gain#
A laser has two basic forms, an amplifier and an oscillator. A laser amplifier is usually a single or double pass device in which an input laser is amplified by spontaneous emission in the gain material in which there is a population inversion. If the gain material contains atoms, ions, or molecules this inversion contains many more excited states than ground states and an incoming photon of the correct energy will stimulate an excited state to enter a lower state by releasing its energy as a new photon. The population of this lower state has to be kept close to zero if the laser is to work well, because the transition rate is proportional to the population difference. The initial inversion is produced by an external source, for example, by using another laser, flash lamps, or an electric current.
In a laser oscillator, the gain material is placed between mirrors so that feedback of photons between these is possible. The lasing builds up to a steady state because a little of the spontaneous emission from the population inversion is captured by the mirrors and repeatedly passes back and forth through the gain material causing each photon to stimulate another every time it does so. The process is therefore non-linear. If you have tried to align a laser, you will have noticed that as soon as the alignment is correct, and the losses are therefore drastically reduced, the laser instantly lights up. A steady state is reached because there are losses in the cavity as well as gain and these balance one another.
Fig 3 Gain threshold and stability in a basic model of a laser when \(a \gt\) 0. \(n\) is the number of photons.
A basic model of a laser measures the rate of change of the number of photons \(dn/dt\) as the difference between the number due to gain and loss and, clearly, if the laser is to work the gain must initially exceed the loss. The threshold to lasing occurs when gain equals loss, (see Haken 1978, p. 127; and Svelto 1982). The increase or gain in the number of photons is due to stimulated emission. The losses are caused by reflections from surfaces in the laser cavity, such as those on the laser rod or the quartz of a dye cell, and from transmission through the output-coupling mirror. Some gain materials also absorb at the same wavelength as the laser operates, and this is a cause of loss. The gain minus loss in the number of photons is therefore
where \(g\) represents the gain coefficient, \(N\) the number of excited states in the gain medium, which assumes that the lower state has zero population. The rate constant for decay of photons out of the cavity \(k\) accounts for all the losses. This equation is non-linear although it appears not to be. The non-linearity is introduced because the number of excited states \(N\) is not expected to be constant since photons are stimulated out of it. Suppose, therefore, that by continuous external excitation, the number of excited states is kept constant at \(N_0\), then \(N = N_0 - \gamma n\) where \(\gamma\) is a constant. The rate equation now becomes non-linear;
which can be simplified to
with the constant \(a = gN_0 - k \) and \(b =g\gamma\).
If the excitation producing the excited states is weak, then \(a\) is negative and the number of photons does not increase. The lasing threshold condition occurs when \(a = 0\) giving the minimum \(N_0\) as \(k/g\). Thus, if the laser gain coefficient \(g\) is high, the lasing threshold is small. If the rate of loss of photons in the cavity \(k\) is small, which corresponds to a long cavity lifetime, the threshold is also small. If excitation is strong, making \(N_0\) large, then \(a\) is positive and the laser operates.
The phase portrait \(dn/dt\) vs. \(n\) is now an inverted curve as \(a \gt 0\) and has steady states at \(n_{ss} = 0\) and at \(n_{ss} = \alpha/\beta\). Since the rate of change of the number of photons has to be positive for the laser to work, the gradient at threshold must be negative; see Fig. 3. The second derivative at threshold is \(-a\), which is negative, and confirms the stability as \(a\) is positive. The presence of the steady state means that the laser intensity will rise to a constant value, given by the pump intensity. When this is achieved, the gain must balance the loss because the rate of change of the number of photons is zero. Note that loss includes the number of photons that form the laser beam itself.
The rate equation is solvable by separating variables and integrating;
The integral can be separated into simpler terms using partial fractions, Chapter 4.2.13, or SymPy used. The result is
Hence
where \(c\) is determined by the initial conditions. This function is plotted in fig 4 and the solution obtained using SymPy is also given.
Fig. 4 Relative population \(n(t)/n_0\) vs. at for a laser with \(a/(bn_0) = 2\) and \(1/2\).
The calculation is given below
# laser gain
a, b, t, n0 ,C1 = symbols('a, b, t, n0, C1', positive = True)
n = Function('n')(t)
f01 = diff( n,t) -a*n + b*n**2
print('equation is ')
f01
equation is
ans = dsolve(f01,n)
ans
Solving for \(e^{aC_1}\) with the initial conditions \(n(0)= n_0 ; t=0\) produces
If \(n(t)\) is plotted vs time, the curve rises or falls to reach a steady state number of photons, \(n_{ss} = a/b\). This is reached because the exponential terms tend to zero as \(t\) increases. This steady state can be greater or smaller than \(n_0\) because \(n_{ss} = (gn_0 - k)/g\gamma\) and when the gain \(g\) is large and \(k\) the rate of loss of photons is small, then \(n_{ss} \gt n_0\) and the number of photons increases from \(n_0\) to the steady state value. The threshold condition is \(k/g\) and this needs to be small for lasing to occur. The expression for the steady state also shows that the number of photons increases linearly with pump power, which is proportional to \(n_0\). Thus, by examining the phase portrait, the same conclusions are arrived at as by solving the equations and then looking for the long time limit.
The growth or decrease in relative population \(n(t)/n_0\) is shown in Fig. 4 plotted vs. reduced time \(\tau= at\) and with dimensionless parameter \(a' = a/(bn_0)\). Making these changes reduced the equation to a simpler form;
3.2 Autocatalytic reactions#
The simplest autocatalytic reaction is
and this will lead to the characteristic ‘S’ shaped curve of the Verhulst or Logistic equation. The form of the rate equations describing the laser is mathematically the same as for this autocatalytic reaction; this is not so very surprising because in a laser one photon stimulates another. An autocatalytic reaction is also described in question Q18.
The rate equations are
and the method of solution given in the answer to Q18. To achieve this we use the initial conditions, \(A=A_0,\,B=B_0\) at \(t=0\) and because
Letting \(x=A\) at time \(t\) substituting for \(B\) gives
Integrating and rearranging gives
which shows that the amount of \(A\) is initially \(A_0\) but tends to zero at very long times. A graph of this function is in the answer to Q18 and is one form of the Logistic equation. The solution for \(B\) has the familiar sigmoid typical of the Logistic equation, starting with a small value but becoming constant at long times as shown in figure 4a below. This data shows the catalytic conversion of trypsinogen into trypsin, Kunitz & Northrop, J. Gen. Physiol. v19, p991, 1936.
Fig 4a. The data points show the conversion of trypsinogen, A in our model, into trypsin B. The solid line is a fit to \(B\) using the last equation above knowing the concentrations \(A_0, B_0\) and by varying the rate constant \(k\). The values used were \(A_0 = 72, B_0 = 0.1, k = 0.0175\) /hour/concentration. (Data points from Kunitz & Northrop, J. Gen. Physiol. v19, p991, 1936)
3.3 Population dynamics#
How a population of animals or humans changes with time was first examined mathematically in 1798 by Thomas Malthus. He argues that if in a small time interval \(\delta t\) the chance of an offspring being born is \(b\delta t\), and in the same interval the chance of death is \(d\delta t\) and assuming that the population is large so that the number of individuals can be represented as a continuous variable then
subtracting \(N\) from both sides and dividing by \(\delta t\) then in the limit make \(\delta t\to dt\) the differential equation becomes
where \(N_0\) are the number at time zero. Clearly when births exceed deaths the population grows without limit, a ‘Malthusian Catastrophe’. The quantities \(b,d\) in chemical parlance are first order rate constants so have units of 1/time. The ratio \(R_0=b/d\) is called the reproductive ratio, as clearly when \(R_0>1\) the population will increase.
This model is clearly a good start but not that good as it allows an semi-infinite number of individuals to exist without regard for vital things such as availability of food and water. Vanhulst made a modification to the rate equation to account for limited resources that all species have to confront. The change was to add an extra term to the rate equation with a parameter that constrains the population because of environmental effect, limited food, excess heat or cold etc., the equation which is quite empirical, becomes
where \(r=b-d\) is called the nett per capita growth and \(K\) is the environment carrying capacity. An important point to notice is that the growth rate decreases as the population increases an effect that is called compensatory. This equation has the form of that discovered for an autocatalytic reaction and for laser gain. It turns out that this is a poor model for estimating human populations but really quite good for predicting bacterial ones. The population is given by
3.4 Harmonic oscillator#
Simple harmonic motion is described by the generic equation
where \(\omega\) is the oscillator frequency, \(x\) position, and \(t\) time. As the phase portrait is a plot of \(dx/dt\) vs \(x\), the equation has to be integrated once to get it into this form. This can be done by defining \(v = dx/dt\), then
Next using the chain rule
and the equation becomes,
Integrating by separating variables gives
where \(c\) is a constant. The trajectories in the phase plane are ellipses whose exact values are determined by their initial conditions. However, all of the trajectories take the same time to complete, as only one frequency is associated with the oscillator.
The origin of the phase portrait is a position of stable equilibrium, and corresponds to the oscillator being stationary. If a pendulum is observed, with \(x\) as the initially positive angle and with an initial velocity of zero, then the starting point is on the positive x-axis. The phase plane is followed clockwise, as the pendulum increases in (negative) angular velocity to the left and the angle decreases towards zero, at which point the pendulum is vertically downwards and the angular velocity has its maximum negative value. The pendulum continues to move until the velocity is again zero at the opposite angle to the starting point. It now changes direction, the velocity reversing (becoming positive), and returns to the starting point, travelling over the top of the phase plane.
Fig. 5 Phase portrait for the simple harmonic oscillator. The position of a pendulum at various points is shown. Motion to the left is a negative velocity.
As the velocity has been calculated if this is integrated the time to move from one position to another can be now be found using
and then
where \(c\) is proportional to the energy. Integrating produces
Suppose that the oscillator is a diatomic molecule with frequency \(\omega\) (in radians/s) where
with \(k\) as the force constant and \(\mu\) the reduced mass. The energy is
with quantum number \(n=0,1,2\cdots\). Substituting the energy for \(c\) gives
and to find the period of oscillation the end points of the oscillation must be found. These occur when the bond is most compressed and when most stretched and at these points the velocity is zero, hence at these points
which gives the connection between bond extension \(x\) and \(n\) quantum number. Substituting for both limits with
and multiplying by two to make a round trip gives the round trip time as
which shows that the period \(t\) is independent of the quantum number, just as expected for the harmonic oscillator.
3.5 The Morse anharmonic oscillator#
The Morse potential is often used as a model of a diatomic molecule as an anharmonic oscillator because not only does the potential allows the molecule to dissociate at high energy, unlike the harmonic oscillator so is inherently more realistic but also the Schroedinger can be solved analytically with this potential. The energy (eigenvalues) are
where \(\omega_e=\sqrt{k/\mu}\) is the fundamental oscillator frequency and \(x_e\) the (dimensionless) anharmonicity. If the bond extension about the equilibrium point ( the bottom of the potential well, \(r_e\) ) is \(q = r-r_e\) and the dissociation energy \(D_e\) the potential has the form
where \(a=\sqrt{k/2D_e}\) and has dimensions of 1/length. The force constant is \(k\) and is that value that would pertain just at the minimum of the potential well where the motion is effectively harmonic. Figure 5a shows the potential and some energy levels. That the energy levels become closer together as the bond stretches can be understood by recalling that the same effect occurs in a ‘particle in a box’ as the box is made larger and is a general effect in quantum systems.
Figure 5a. Example of a Morse potential with several energy levels. Parameters used are not from any particular molecule and are \(D_e = 1,\; a=1,\; \mu = 30\).
Following the example of the harmonic oscillator (ii) and assuming the motion to be classical the equation of motion becomes
which can be integrated by separating variables to find the velocity
and from this equation the velocity vs. position can be plotted for different values of the constant \(c\). Below the dissociation limit \(D_e\) these values are given by \(E_n\) for quantum number \(n\) and some values of the momentum \(p=\mu v\) vs. position \(q\) are shown in figure 5b. The dot at (0,0) shows the minimum of the well and two quantum numbers \(n = 0,6\) are labelled. When the energy is below \(D_e\) the motion is periodic, i.e. the molecule is vibrating. At \(q = 0\) molecule is at the middle of the potential well and its momentum (or velocity) is at a maximum, whereas at \(p = 0\) the molecule is at the inner or outer turning point and is instantaneously stationary. Above \(D_e\) the motion is no longer periodic, i.e the molecule is dissociating and a set of equally spaced values of \(c\) are shown.
Figure 5b. Momentum vs. position for the Morse potential. The coloured lines are below the dissociation energy, the blue one at energy \(D_e\) (the separatrix ) and the grey ones above this energy.
When there are very many levels in the potential, such as when \(m\) is large then it is possible to excite the molecule to very close to the dissociation limit. In this case molecules can have a huge size, many hundreds of nanometres if not larger. These can be called ‘Rydberg’ molecules.
3.6 Hamilton’s Equations of Motion#
The calculation of the period as was done with the harmonic oscillator is, in principle, just the same: the equation for the velocity is integrated by separating variables. However, in this case a very hard integral is produced and another and far more powerful approach can be used. Thus uses Hamilton’s formulation of mechanics and is in itself a very powerful method but is huge topic in itself and so only an brief example will be given here.
The basis of the method we use is based on two equations connecting the total energy (in this example) with momentum and position. The total energy is called the Hamiltonian, \(H\) and position and momentum are conventionally labelled \(q,\;p\) respectively.
Hamilton’s ‘canonical’ equations of motion are;
where
The total energy of the oscillator is the sum of the kinetic and potential parts, thus with mass \(m\) and kinetic energy \(p^2/2m\)
Using this equation it is clear that a contour plot of \(p\) vs \(q\) can be immediately produced, where the values of the contour can, for example, be chosen to be the energy of the quantum oscillator e.g.
which is only valid up to the dissociation energy.
The period of the motion at each energy \(E_n\) can be found by integrating \(dq/dt\) and this equation is
which means finding an equation for the momentum. However, we already have this as it is the Hamiltonian, all that is needed is to equate it with an energy and we can choose \(E_n\). Thus
and rearranging
integrating between positions \(q_1,q_2\) gives
The integral is easier using the substitution \(u=e^{-ax}\) which produces
When expanded the left-hand side has the standard form
with \(a\to -D, b\to 2D,c\to E_n-D\) giving
The limits, which are the turning points have to be found. As these occur when the momentum is zero (bond fully stretched or compressed) these are found when
which produce limits in terms of parameter \(u\) of
The conversion
is also needed.
The term \(\sqrt{E_n -D_e(1-u)^2 }\) in the denominator becomes zero with these limits, making \(\tan^{-1}(\infty)=\pi/2\), and the solution is greatly simplified yielding a period of
As a check work out the dimension of this result; \(a\) has dimensions of m\(^{-1}\), \(D\) of energy (kg (m/s)\(^2\)), thus \(t\) has units of time. As expected for an anharmonic oscillator, the period become bigger as the energy increases. The excursion the atoms make is also increased as may be seem in figure 5b. When \(E_n=D\) the period is infinity as here the molecule has dissociated. This trajectory is the blue line in fig 5b and is called a separatrix or the homoclinic orbit .
Figure 5c. The period of an anharmonic oscillator vs energy. The parameters are as used in figure 5b. The dissociation energy \(D=1\). The circles show the energy \(E_n\) for \(n=0\cdots 9\) the only bound energy levels in this particular potential.
4 Chemical Kinetics 2#
Here we describe a few examples of chemical kinetics including those reaching equilibrium, reaction with flow, and reaction and diffusion. In sections 9.2.iv,v and vi of this chapter and Chapter 11.4.10, 11.7 and 11.8 (Numerical methods) examples of more complex reactions with diffusion or multiple steps, which can generally only be solved numerically, are examined.
4.1 Bimolecular reactions#
In a bimolecular reaction, such as A + B \(\rightarrow\) C, the rate of reaction of A would normally be written as
and as it is written it is of little help in solving the equation. If the initial amounts of A and B are \([A]_0 = a_0\), \([B]_0 = b_0\), and if \(x\) moles have reacted at time \(t\), then the amount of A at any time is \(a_0-x\) and the rate equation can be solved as
The variables are separable, and to integrate the right-hand side this can be split into partial fractions,
integrating and combining terms produces
and the constant \(c\) can be evaluated because \(x = 0\) at \(t = 0\), therefore
which can easily be solved for \(x\), the amount reacted, at time \(t\) by first exponentiating the rearranging.
In the case when \(2\mathrm{A}\rightarrow \mathrm{B + C}\) the rate equation is
Integrating with \(a=[A],\; a_0=[A]_0\) at \(t=0\) and \(b_0=[B]_0\),
and so the product follows
and sometimes the factor of \(2\) is incorporated into the rate constant.
The calculation can also be done by letting \(x\) be the amount that reacts, i.e. amount of product formed, then
Integrating gives \(\displaystyle \frac{x}{(a_0-x)}=2a_0k_2t\) which on extracting \(x\) gives the equation above for \(B\).
4.2 Reversible or opposing reactions \(A \overset{k_1}{\underset{k_2} \rightleftharpoons }B+C\)#
Suppose that the reaction
reaches equilibrium there being an amount \(\mathrm{[A]_0}\) initially and \(\mathrm{[B]_0=[C]_0}\). The simplest way to integrate the equation is to let an amount \(x\) react then
The equation is
letting \(b_0=c_0=0\) substituting for \(x\), the amount reacted at time \(t\), gives
This equation has a standard form
where \(\alpha = k_1a_0,\; \beta=-k_1,\; \gamma=-k_2\) and the condition is that \(\beta^2-4\alpha \gamma>0\) which is true. Substituting produces
and also evaluating the constant, gives
and from this we can see that the reaction comes to equilibrium with a first order rate constant of
where \(\tau\) is the relaxation time or simply the lifetime. It would be a mistake to think that the decay was exponential, however, as this is true only at longer times, see fig. 5d.
Making substitutions \(\displaystyle w=\frac{k_1+q}{2k_2},\;v=\frac{k_1-q}{2k_2}\) and exponentiating gives
with limits at \(t=0,\; x\to 0\) and at long times \(t\to \infty,\; x\to -v\) and \(v\) is negative because \(q>k_1\).
Figure 5d. Comparison of the true decay (circles) vs. that by making the approximation that \(x^2\) is small (solid line). The equilibrium (long time) amount has been subtracted from the true decay (circles) eqn. 5c, to make the form of the decay clearer. The straight line is eqn. 5d and is the same as produced when \(x^2\) is small as described in section 4.3 where a temperature jump is described.
This form of the equation is clearly rather messy. It is simplified if the equilibrium amounts are incorporated. At equilibrium the rate of change of A becomes zero and so
where \(x_e\) is the equilibrium amount and so
The equilibrium constant is by definition
Substituting for \(k_2\) and after a great deal of simplifying the result is
and as \(x,x_e,a_0\) are all concentrations this is dimensionally correct.
4.3 Temperature Jump; a perturbation from the equilibrium \(\mathrm{H_2O}\overset{k_1}{\underset{k_2} \rightleftharpoons }\mathrm{H^+}+\mathrm{OH^-} \)#
Once a reversible reaction such as
or specifically
comes to equilibrium a small perturbation, such as a temperature jump, can be used to measure the return to equilibrium. This may take only a few microseconds, and from this the rate constants \(k_1,k_2\) can be obtained. After a perturbation the change in absorption in the visible or IR, or fluoresce is commonly measured, but sometimes optical rotation or light scattering can used to instead. Often the perturbation is in the form of rapid heating of the solvent by a degrees or so, called T-jump and by using this method numerous enzyme equilibria have been studied. In the past heating was achieved by passing an electric current but today this can be effected most easily by using a nanosecond duration laser tuned, for example, to water’s \(1.9\) micron absorption in the infra red. This can be achieved via the Raman effect by passing a nanosecond duration \(1064\) nm YAG laser through a gas, such as hydrogen, at very high pressure, \(20 - 60\) bar. Alternatively a dye with a low fluorescence yield, such as crystal violet, is added to the solution. This dye will absorb visible laser light which will indirectly heat the solution as its excited state decays.
Consider the general equilibrium
and let \(a_e, b_e, c_e\) be the equilibrium amounts of A, B and C respectively. An amount \(x\) is transferred to products as a result of heating the solution and the concentration of \(a\) is \(a = a_0-x\). The rate equation is
and at equilibrium
where \(x_e\) is the equilibrium amount of A. A change in amount \(x\) due to the temperature jump is \(\Delta x=x-x_e\) and substituting for this gives
This equation has the form \(\displaystyle \frac{d\Delta x}{dt}=a+b\Delta x+ c\Delta x^2\) where \(a,\;b,\;c\) are constants which can be integrated using a standard form, see the previous section (4.2) or chapter 4-2.14, however, as \(\Delta x\) is small the term in \(\Delta x^2\) is very small compared to the others and may be ignored giving a far simpler equation which fig 5d shows is only valid at longer times,
which can be simplified using the equilibrium condition
to
which on integrating gives
when \(t=0, \Delta x=\Delta x_0\) so that
which has a relaxation lifetime of \(1/(k_1+2k_2x_e)\). If the equilibrium constant is also known then both rate constants can be evaluated because \(K_e=k_1/k_2\).
In water at pH = \(7\) and \(298\) K the relaxation time is measured to be 37 microseconds. The equilibrium constant is
where pure water is \(55.55\) molar. In using the equilibrium constant it is understood that each concentration is divided by 1 molar so that it is dimensionless. At this temperature \(K_w = 10^{-14}\), thus \(x_e=10^{-7}\) and in the denominator can be ignored giving \(K_e=1.80\cdot 10^{-16}\).
The equation is \(\displaystyle \frac{1}{37\cdot 10^{-6}}= k_2(K_e+2x_e)= k_2(1.80\cdot10^{-16}+2\cdot 10^{-7})\) giving \(k_2=1.35\cdot 10^{11} \mathrm{dm^3\;mol^{-1}\;s^{-1}}\), and therefore, \(k_1=2.43\cdot 10^{-5}\mathrm{s^{-1}}\).
The large value for \(k_2\) suggests that it is diffusion controlled. This may be estimated using the Stokes-Einstein equation
which gives a value of \(6.6\cdot 10^9 \;\mathrm{ dm^3\,mol^{-1}\,s^{-1}}\) in water with a viscosity \(\eta=0.001\) Pa S. This is clearly way-off the measured value, but the formula was derived for neutral species and therefore ignores the fact that ions are involved which in this case are attracted to one another. Adding a correction (Wilkinson, 1980) increases \(k_2\) by a factor of around twenty (depending on ion separation) bringing the value towards \(10^{11}\;\mathrm{ dm^3\,mol^{-1}\,s^{-1}}\) much closer to the experimental measurement.
4.4 Reversible reactions \(A + B \overset{k_1}{\underset{k_2} \rightleftharpoons } C + D\)#
The gaseous decomposition of HCl follows follows the stoichiometry \(\mathrm{2HCl \rightleftharpoons H_2 + Cl_2}\) although the reaction is more complex involving radicals. This is a particular case of a reversible reaction \(A+B=C+D\). The scheme general is, using the same notation as above,
and again the integral
has the standard form of eqn 5a with \(\alpha = k_1-k_2,\; b= -2a_0k_1x,\; c=k_1a_0^2\)
which shows that equilibrium is reached with a rate first order constant \(2a_0\sqrt{k_1k_2}\) because \(k_1,k_2\) are both second order rate constants.
Again, using the equilibrium condition, which in this case is
the equation can be simplified to
A reaction such as \(\mathrm{2HCl \rightleftharpoons H_2 + Cl_2}\) can be converted into \(A+A\rightleftharpoons C+D\) quite simply by using the form of the rate equation for \(A+B\rightleftharpoons C+D\) given above but letting \(k_2'\) be the new back-reaction rate constant, i.e.
This is has the form of \(A+B\rightleftharpoons C+D\) and the equations become the same if we make \(k_2'=4k_2\).
4.5 Hydrogen iodide equilibrium#
One of the first reaction to be studied in the late 1800’s was the dissociation of hygrogen iodide vapour at high temperatures \(\approx 450\) K. It is interesting to derive the kinetic equations for
from scratch. The rate equations are
Suppose that initially only hydrogen and iodine are present in amounts \([H_2]_0, [I_2]_0\) then
thus by substituting so that the expression is only in HI, the other concentrations being constant,
If initially only HI is present then \( [HI]_0-[HI]=2[H_2]=2[I_2]\) then similarly calculating
which leads to the same rate equation when the numbers of H and I atoms are equal, i.e. \([H_2]_0=[I_2]_0\) and \([H_2]_0+[I_2]_0=[HI]_0\) and so the same equilibrium is reached whether the reaction stats as pure HI or as hydrogen and iodine.
At equilibrium the forwards and backwards rates are equal and then the equilibrium condition is
and solving for the equilibrium constant \(\displaystyle K_e=\frac{k_1}{k_2}\) gives
To solve eqn 5f a change in notation clarifies the equation, let
making the solution the integral,
which has a similar solution to those given above but is more compact if expressed in terms of \(\tanh^{-1}\) solution where
The solution is
Using the identity
and isolating \(x\) gives
and starting with iodine and hydrogen only and integrating from \(1\to x\), gives
The approach to equilibrium both starting with HI and with Hydrogen and iodine is shown in figure 5e. The data was recorded in 1899 by Bodenstein and are replotted below. The temperature was \(448^\text{o}\) C. The lines are fitted from the integrated equations. The quantity labelled as \(x=[HI]/[HI]_0\)
Figure 5e. \(\mathrm{2HI \rightleftharpoons H_2 + I_2}\) reaching equilibrium starting at zero HI, rising points, and eqn 5g, and only HI, falling points, and eqn 5h. The equilibrium value \(x_e=0.786\) at \(448^\text{o}\) C. The curves show that the same equilibrium is reached starting either with products or reactants.
4.6 Reversible reaction when only the equilibrium constant is known \(A \overset{k_1}{\underset{k_{-1}} \rightleftharpoons } B\)#
Sometimes we may want the rate expression in terms of the equilibrium constant as this may be known when neither rate constant has been measured. Examples are the cis - trans isomerisation of some hydrocarbons and of \(\alpha -,\beta-\) glucoses and the ortho to para conversion of H\(_2\) with a catalyst such as Nickel on \(\mathrm{Al_2O_3}\). This reaction has the form
Letting \(a\equiv o-H_2\) and \(b\equiv p-H_2\) concentrations produces the rate equation
The initial amounts are \(a_0,b_0\) and as the total amount of ortho- and para-hydrogen is constant at all times
where subscript \(e\) denotes the equilibrium. The equilibrium constant is
and by rearranging and substituting for \(K\), the rate expression becomes
Using the total amount
and substituting again gives
which is now a standard integral with \(a\) being the only variable
producing
or
which illustrates an exponentially decaying population of \(a\) (ortho hydrogen) from \(a_0\) to reach a constant value of \(a_e\) at equilibrium with a measured rate constant of
Both \(k_1,k_{-1}\) can be found if \(K\) is known, for example from thermodynamic data via the standard Gibbs energy change,
Notice that the initial amount of species \(b\) is irrelevant to the decay of \(a\).
4.7 Kinetics of ligand exchange using isotopes#
A number of reactions involving electron and ligand exchange can be followed using isotopes ( H. McKay J. Am. Chem. Soc. 1943, v 65, p 702). Examples are
where the \(*\) represents the radioactive tracer. In both reactions \(\Delta H=0\) but \(\Delta S \gt 0\) because of the scrambling of isotopes between the products. As \(\Delta H=0\) this does not means that there is not an activation energy for the process. The energy barrier in the second example is caused by the different bond lengths when an electron is transfered between \(\mathrm{Fe^{II}-OH_2}\) and \(\mathrm{Fe^{III}-OH_2}\) the former being \(0.22\) nm the latter \(0.205\) nm. One bond has to stretch the other contract and this takes energy. There will also be a contribution from re-arranging water molecules surrounding the ions as the charge changes. The reaction is followed by precipitating out some \(\mathrm{Fe^{3+}}\) to stop further reaction, and measuring its radioactivity. (We shall assume that any decay during the reaction is insignificant, although if it is not, it can be accounted for knowing the half life).
The general equations are
and let
At \(t=0\), \([A^*X]=0\) and at infinite time we measure \(x_{\infty},y_{\infty}\) and as mass is conserved
At the end of the reaction the tracer is distributed in proportion the the number of moles of each that are present, thus
The reaction scheme only gives the stoichiometry and does not indicate any mechanism, first or second order for example, so instead of using rate constants we use a proportionality constant \(R\) which is the gross rate of exchange of ligands, which is constant because the chemical composition is constant, and has units of concentration/time. If the reaction is in fact bimolecular then \(R=k_2ab\).
The rate of production of \([A^*X]\) is
therefore
We shall need to eliminate \(y\) before integration. To do this start by rearranging,
where \(y/b\) is the fraction of \([B^*X]\) and \(x/a\) for \([A^*X]\), and substitute for \(y\) using \(x+y=x_{\infty}+y_{\infty}\), which leads to
Separating variables and integrating gives
which is
The measurements produce \(x/x_{\infty}\) at different times so a log plot has a gradient of \(-R(a+b)/(ab)\) from which \(R\) is found as \(a,b\) are known. The form of \(R\), for example \(k_2ab\) or \(k_1a\) is found by repeating the experiment with different values of \(a\) and \(b\).
4.8 Kinetics and flow#
Water entering and leaving a tank, reaction vessel, or even a lake can be modeled by calculating the difference in the amount material flowing in and out, viz.,
in an analogous way to a chemical reaction. Terms for chemical reactions can be added if needed. It has to be assumed not only that the rate of inflow is the same as outflow, otherwise the lake will fill or empty, but also that material entering is fully mixed before it leaves (a well-stirred reactor), otherwise a concentration gradient would exist and this will complicate the calculation. If the concentration is \(x_0\) before any flow starts, the material entering has a constant concentration \(x_{in}\), the in-flow rate is \(f \,\mathrm{dm^3\,s^{-1}}\), the volume of the lake is \(V\) and if \(x\) is the concentration (mole dm\(^{-3}\)) of the tank/lake at any time \(t\) then,
and notice that \(f/V\) has dimensions of a first-order rate constant or s\(^{-1}\). Solving this initial value problem by separating variables gives
which has the solution
as \(x=x_0\) at \(t = 0\) then
Working in reduced values makes it simpler to determine the range of behaviour of the results, see Fig. 6. If the initial concentration \(x_0 \ne 0\) the reduced equation can be written as
where
so that now only the ratios need be considered. If \(r_{in} \gt 1\), the concentration ratio increases with time and the lake becomes polluted. At values of \(r_{in} \lt 1\), which means that \(x_0 \gt x_{in}\), the lake becomes more dilute and if pure water is added, \(r_{in} = 0\), the polluted lake would be cleaned after a time larger than approximately \(t = 5V/f\).
If the initial amount could be zero, i.e. the lake is clean but becomes polluted then we can write
making
which always reaches a constant value of \(1\) but with initial value depending on \(r_0\).
Barnes & Fulford (2002) discuss several similar models.
Fig. 6. Left. The change in the ratio \(x/x_0\) vs. reduced flow \(ft/V\) for different \(r_{in} = x_{in} /x_0\) ratios, showing the range of curves obtainable with different starting conditions. Right. The same equation but rewritten with ratio \(r_0=x_0/x_{in}\). The ratios in both plots are \(0, 0.5, 1\) and \(1.5\).
4.9 Dissolution kinetics#
(i) Dissolving#
When a solid solute is dissolved in a solvent, the rate equation is found by considering the change in the amount dissolved in solution during a short time period. In dissolving a solid, a saturated solution will eventually be formed and this limits how much solid will dissolve. If \(x_0\) is the initial amount of solid to be dissolved in \(m\) grams of solvent, \(s_x /m\) the saturated solution concentration, \(k\) the rate of dissolution (mass s\(^{-1}\)) and \(x\) the number of grams of solid remaining at time \(t\), then
\(\qquad\qquad\) Amount of \(x\) dissolved in time \(t + \delta t\) = amount \(x\) at \(t - \) amount dissolved in \(\delta t\).
The amount dissolved in time \(\delta t\) is proportional to the amount of solid undissolved \(x\), multiplied by the difference in concentration compared to that of a saturated solution, and this product is
The term \((x_0 - x)/m\) is the concentration dissolved, making the rate equation,
To solve, separate variables, and letting \(a = x_s - x_0\) gives
which is a standard integral as \(a\) is a constant. If this is not recognized, then converting to partial fractions gives,
and the integral is
If \(x=x_0\) at \(t=0\) then
and therefore
which can be solved for \(x\) if the concentration profile with time is needed. This result shows that the amount of solid \(x\) decreases exponentially to a constant value determined by \(x_s - x_0\) from an initial value \(x_0\) provided \(x_s \lt x_0\), otherwise the amount of solid material becomes zero.
(ii) Spreading of disease#
A rate equation with the form of equation (6) also appears in a completely different context, which is that of spreading a disease among \(N\) individuals if the rate of spread is proportional both to the number infected \(x\) and the number who are not infected. The equation then has the form
where \(k\) is the rate of spreading the infection in units of number of individuals time\(^{-1}\). Integrating gives
4.10 Mean free path and probability of obtaining a path of a given length#
If \(\bar c\) is the average velocity and there are \(v\) molecules in unit volume then the total distance covered is \(\bar c v\). Let \(\gamma\) be the number of collisions that occur per unit volume per second and as each collision terminates two free paths, there are \(2v\) free paths in the same time. As the average distance travelled is \(\bar c v\) the average length of the free path \(\lambda\) is
from which the number of collisions can be calculated. The mean free path can be determined by thermal diffusion or the viscosity of a gas, and was one of the early successes of the kinetic theory of gases.
The chance of a collision in a sufficiently small length \(dL\) is proportional to that length and is independent of the path \(L\) already travelled, thus we can make this chance \(\beta dL\) where \(\beta\) will be determined later. After a collision let the probability equal \(p(L)\) that a path of length \(L\) occurs. At this distance the chance of a collision in the small increment \(dL\) is \(\beta dL\) and the chance that no collision will occur is therefore \((1-\beta dL)\). The chance that a path of length \(L\) will occur and then no collision in the next increment \(dL\) is the product of the two probabilities because they are independent of one another and is
The quantity \(p(L+dL)\) has to be expressed in terms of \(p(L)\) to be of any use. In the diagram below \(dL\) is an infinitesimal. The gradient is
and so the increment in \(p\) is \(\displaystyle \Delta\) so that
therefore
giving,
which has the solution \(\displaystyle p(L)=Ae^{-\beta L}\). When \(L=0\) the constant \(A\) is unity so that \(\displaystyle p(L)=e^{-\beta L}\).
Figure 6a. Calculating \(p(L+dL)\).
To determine the constant \(\beta\) we need to average over all paths. The chance of having a path from \(L\to L+dL\) is the product of the chance of a path of length \(L\) and of a collision in \(dL\) because these are independent events. This product is \(\beta p(L)dL\). The mean free path is the average value of \(L\) but weighted with the number of free paths of different lengths, i.e. the distribution of path lengths, or
thus
and shows how rapidly the chance of mean free path \(L\) decreases at fixed \(\gamma/v\) and also as \(\gamma/v\) increases at fixed \(L\).
5 Diffusion of heat and molecules#
At steady state, the constant quantity of heat, or the heat flux, \(Q\) in watts ( J s\(^{-1}\) ) flowing through an area \(A\), is proportional to the spatial rate of heat loss or
where \(\theta\) is the temperature, \(x\) the thickness, and \(k\) the thermal conductivity coefficient. This equation is called Fourier’s law of heat conduction. The heat flux \(Q\) is also \(dH/dt\) = constant, if \(dH\) is the constant amount of heat transferred in time \(dt\). In general, the rate of heat transferred from one body to another will depend on the shape of the body and its composition, as well as the temperature difference and whether these are held constant. Therefore, the flux and \(d\theta/dx\) may change across an object. This detail, however, adds nothing new, per se, to the problem. The corresponding steady state equation for the mass diffusion of fluids such as gases and liquids, is Fick’s first law, which is
where \(J\) is the flux, which is the amount of material diffused/unit area/unit time and this is a constant quantity. The diffusion coefficient is \(D\) (\(\mathrm{m^2\,s^{-1}}\)), \(c\) is the concentration, and \(x\) the distance over which diffusion occurs. If these last two ‘diffusion’ equations were to be rearranged, they would have the mathematical form of zero-order rate equations because the rate of change does not depend on \(\theta\) or \(c\), but is a constant quantity for example,
An example from chemical kinetics of a zero-order rate equation is
5.1 Heat loss#
The heat and diffusion equations can be used to solve a variety of different problems. For example, a fridge has a wall that is \(5\) cm thick, and is kept \(20^\mathrm{o}\) cooler than the room. The thermal conductivity coefficient \(k = 0.1 \mathrm{J\, s^{-1}\, m^{-1}\, K^{-1}}\) ( or \(0.1\) watt/metre/ kelvin ), which is typical of good insulating materials, then the steady heat flow into the fridge is
for each m\(^2\) of surface area. Because the temperature change is fixed as is the wall thickness, then
5.2 Ice forming#
The following example describes ice forming on a still lake. The ice layer increases in thickness with the square root of time when the water temperature is \(0^\mathrm{o}\) C and the air temperature is lower and constant. The ice acts to insulate the water from the colder air above. A quantity of heat, \(dH\), is taken from the water to freeze it in time \(dt\), and the rate of heat transfer is \(dH/dt\). By Fourier’s law, equation (8) this is proportional to the temperature gradient \(\Delta \theta/x \) across the ice, therefore \(dH/dt \propto \Delta\theta/x\). As the ice thickens, the temperature gradient decreases and so heat transfer is reduced per unit time, which must be proportional to the thickness \(dH \propto dx\), hence \(dx/dt = \alpha\Delta\theta/x\) where \(\alpha\) is a constant of proportionality. When integrated, this equation shows that the ice thickness increases as the square root of time; \(x \propto \sqrt{t}\). This, and the fact that ice is less dense that water means that lakes do not generally freeze solid; there are usually not enough cold days.
An interesting consequence concerns the freezing of water droplets in the atmosphere. Non-linear freezing means that the outside of the liquid droplet freezes rapidly, forming a shell. This may become strained as it thickens, causing the ice to crack thus releasing a spurt of liquid water from the interior in the form of micro-droplets, which rapidly freeze in the cold air. Such a process may be important in cloud formation.
5.3 Diffusion through a cell wall#
A third example concerns diffusion. A thin porous pipe of outer diameter \(b\) has gas at pressure \(p_0\) flowing inside it. The inside of the pipe has diameter of \(a\) and the concentration of gas outside the pipe is zero. If the diffusion coefficient through the pipe to the outside is \(D\), the rate of gas loss can be calculated using Fick’s first law. The ‘pipe’ could be, for example, a capillary in the lungs containing O\(_2\) and CO\(_2\).
The flux (mass flow) per unit area/second through the pipe is a fixed quantity, which is given by
where \(r\) is the radial distance from the centre of the pipe and \(c\) the concentration of the gas. The thickness of the pipe’s wall is \(b - a\) and the surface area is \(2\pi r\) per unit length, hence
per unit length. Integrating through the wall of the pipe, from \(a\) to \(b\) gives
with the result
where \(q\) is the integration constant. With the initial conditions \(c(b) = 0,\, c(a) = c_0\), the constants \(J\) and \(q\) can be found and eliminated using
which results in the concentration at position \(r\) as,
To find the mass flow \(J\), differentiate this equation with respect to \(r\) and substitute into Fick’s law. The flow rate is therefore,
Does this make sense? Common sense would suggest that the flux of material is going to be large if the diffusion coefficient and initial concentration is large, and if the wall is thin. If the wall is very thin then \(a \to b\) and the log becomes small (\(\ln(1) = 0\)) and the flux increases rapidly as the wall becomes thinner. The effect is quite dramatic; if initially the pipe has \(b/a = 2\) and then this ratio is reduced to \(1.25\), a reduction by a factor of \(1.6\) times, but the flux is \(\approx 6.5\) times greater. If two gases are in the pipe, the ratio of the amount of each transported outside is in the ratio of their diffusion coefficients and their concentrations. This is important for O\(_2\) /CO\(_2\) diffusion in the lungs.
5.4 Reaction and diffusion#
In the presence of a chemical reaction where the molecule is also diffusing, the diffusion equation is
where \(R\) is the rate of reaction, which, for simplicity, is assumed to be a constant, independent of \(c\). This can be the case for O\(_2\) diffusing into a cell, or the reaction of ATP molecules when their concentration is high enough. ATP is used by the molecular motor protein ATPase. In this protein, proton flow across a membrane containing the ATPase drives a rotor that then forces closed one part of the protein and thereby splits ATP to ADP and releases phosphate. The reverse reaction can occur depending on the prevailing conditions. ATP is also used to drive similar molecular motors to ATPase that operate the flagellum of a micro-organism and so enable it to move. In this case, the ATP has to diffuse down the flagellum from its base where it has a constant concentration. The possibility exists that diffusion cannot supply enough ATP molecules, and hence energy, to move the flagellum and using the diffusion equation this condition is sought. The general model is therefore of a fixed concentration \(c_0\) of a species outside a boundary, which diffuses (in one dimension in this model) into the bulk where it also reacts. We want to know to how far into the bulk reaction can occur.
Integrating at steady state when \(\partial c/\partial t = 0\) produces
and if this is zero at the end of the flagellum where \(x = L\), the integration constant \(a\) can be found. Integrating again produces,
if at \(x = 0,\, c = c_0\). At the end of the flagellum, the concentration is
which obviously cannot be negative; the limiting condition is, therefore, that
As a figure of merit, it is found that if \(\displaystyle \frac{RL^2}{DC_0} \lt 1\) diffusion is sufficient to supply enough ATP to the molecular motors or to supply O\(_2\) into cells. Using typical values, reaction can only take place within \(10\) microns at most from the boundary. More sophisticated reactions with, say, spherical geometry lead to similar conclusions. Thus, it is understandable why small insects, for example, have tracheal tubes to increase the O\(_2\) available to diffuse into their muscles and so maintain the high metabolic rate necessary for flight.
5.5 Cooling and heating. Trommsdorff - Norrish effect#
The time taken to cool a cup of coffee or tea follows Newton’s law of cooling. This law states that the rate of cooling depends on the difference in temperature compared to the surroundings and the form that the rate of heat transferred from one body is the same as that gained by another. This means that the rate is proportional to the temperature difference or
where \(t\) is time, \(\theta_s\) is the temperature of the surrounding air, and \(k\) is the constant, measuring the convective transfer of the heat from one medium into the other. When divided by the surface area, this is sometimes called Newton’s cooling coefficient with units of watts \(\mathrm{m^{-2}\, K^{-1}}\). The cooling equation assumes that the cooling is convective and the object is in a slight draught of air so that the air temperature next to the object \(\theta_s\) is constant. If the draught becomes a forced flow of air, then the cooling coefficient becomes dependent on the air speed and a ‘wind chill’ effect occurs. This is soon felt if you are outside on a cool day and move into a strong wind from a sheltered spot.
Note that Newton’s law of cooling has the same form as a first-order rate equation. The rate of change of temperature is proportional to the temperature of the body.
Mammals cool themselves by evaporating water from their bodies, by sweating or panting. The heat balance equation must, therefore, have extra positive terms due to the chemical reactions keeping them alive (their metabolism) and negative terms due to evaporation. This last term could be incorporated into the constant \(k\). Heat flow is also important in simple chemical reactions and exothermic reactions are usually cooled. If the heat loss from the reaction vessel is not sufficient, the heat generated by an exothermic reaction, such as fermentation, can result in spontaneous combustion as used to happen with hayricks.
In radical polymerization reactions, as the reaction proceeds towards completion, the mixture must become very viscous. This means that the termination steps cannot occur and the rate of reaction increases because this is inversely proportional to the rate of termination. Secondly, the increased viscosity means that stirring can become ineffective and the only cooling mechanism is thus thermal diffusion to the vessel walls, which is slow. Consequently, if heat gain is too great, any gases or vapours trapped in the polymer may cause it to explode. This is called the Trommsdorff - Norrish effect and is an example of an auto-acceleration process or one with positive feedback. The rate of temperature change has an extra term for the heat generated and should have the generic form
at temperature \(\theta\). The activation energy is \(E_a\), and \(k_r\) is a constant proportional to the pre-exponential term from the Arrhenius equation and the heat capacity of the reaction mixture.
6 The Centrifuge: Forced separation#
If a cylinder of radius \(r\) and width \(h\) is filled with a fluid, whose several components are to be separated, is spun about its axis at a speed of \(\omega\, \mathrm{rad\, s^{-1}}\), the heavier components are forced to come to equilibrium further towards the outside of the cylinder than the lighter ones do. The balance is between centripetal forces acting towards the axis of rotation and the centrifugal forces acting in the opposite direction.
A small cylindrical shell of fluid in the region \(r\) to \(r + dr\) from the rotation axis, experiences a net force inwards of \(2\pi hr \cdot dP\) where \(P\) is pressure (Margenau & Murphy 1943). This is equal to the inwards centripetal force due to the rotational motion which is \(2m\omega r\). As the mass of the cylindrical shell is \(2\pi rh\rho \,dr \), with \(\rho\) being the density, equating the forces produces
or
In a liquid, the density \(\rho\) is a constant, and integrating produces the pressure
This also explains why, if a liquid is spun in a beaker about its vertical axis, the surface obtains a parabolic profile. This can be seen clearly if two immiscible liquids are spun and one of them is coloured with some dye.
In a centrifuge, imagine strips taken out of the rotating cylinder and each replaced with a tube. Several of these are spun about the axis to balance the rotor. The distance \(r\) is therefore the radial distance from the axis, and \(h\) the height of the tube when it is spinning in a horizontal plane, which is normally its width.
In an ideal gas, the pressure is proportional to the density since
Solving the equation
produces
if the pressure is \(P_0\) at \(r = 0\).
Fig. 7 Horizontally spinning disc; the fluid level is shown by the dashed line.
7 Absorption of photons and scattering of X-rays and electrons#
Photons in the ultraviolet and visible part of the spectrum excite electrons from the HOMO to LUMO orbitals of molecules giving rise to the colours we see around us. Vibrational and rotational transitions are excited by infrared light. When photons are absorbed, if the cross section for absorption at a certain wavelength is \(\sigma\), then \(\sigma ndL\) photons per second are absorbed by \(n\) molecules in a layer of thickness \(dL\). The fractional decrease in intensity of the photons is then
If the initial intensity is \(I_0\) then the transmitted intensity is
which is called the Beer–Lambert law. The absorbed light intensity is
because the total number of photons is those absorbed plus those transmitted. In many instances, the quantity \(\sigma nL\), the optical density, is replaced with \(\epsilon [c]L\) where \(\epsilon\) is the extinction coefficient, which depends on the wavelength and has units of \(\mathrm{dm^3\,mol^{-1}\, cm^{-1}}\), \([c]\) is the solution concentration and \(L\) the path length.
A beam of X-rays can be scattered by the electrons in a molecule’s atoms, rather than being absorbed. In a single crystal, parallel planes of atoms occur at regular intervals. The scattered intensity from the many planes similarly oriented to the X-ray beam can add in phase or cancel out, and this forms the basis of X-ray crystallography. The beam of X-rays has intensity \(I\), and if a sample contains \(n\) electrons per unit volume, the first layer of thickness \(dL\) intercepts \(ndL\) photons, each electron presenting an area of \(4\pi d^2\) where \(d = e^2/mc^2\) and is approximately the closest approach distance of two electrons. An equation similar to the Beer-Lambert equation describes scattering but if this is small, and as \(nL\) is the number of scattering electrons in the sample, the fractional scattering of a single electron becomes \(I_s/I_0 =4\pi d^2=4\pi(e^2/mc^2)^2\), see Q5.
8 Velocity and Acceleration#
A number of problems involving the motion of bodies can be solved by using Newton’s laws to equate forces, and separating the variables of the equation produced. Usually, the force acting downwards is gravity, and it is often convenient to choose the positive direction as downwards. By Newton’s second law of motion, the net force \(f\) acting on a body, is equal to the rate of change of momentum, or
where \(m\) is the (constant) mass and \(v\) the velocity at time \(t\). This equation can also be written as, force equals mass times acceleration or
\(x\) being the position at time \(t\). Both of these equations are of fundamental importance in solving problems in dynamics.
8.1 Skydiving#
Consider now a skydiver who jumps from an aeroplane and falls under gravity, but with air resistance that is linearly proportional to velocity. Naturally, we shall want to work out how velocity increases with time, and if it will reach a constant, limiting value. After a certain time, the parachute must open, and the resistance to falling will increase so that it is now proportional to the square of velocity. At any time the forces due to gravity and air resistance are balanced, and therefore act in opposite directions and this is used to solve the problem. The force of gravity is \(mg, \,g\) being the acceleration due to gravity and \(m\) the mass in kg. The air resistance (drag) is \(av\) where \(a\) is a constant with units of kg s\(^{-1}\) and \(v_0\) is the initial downwards velocity which could be zero. The equation of motion in terms of force is
where the initial condition of velocity being \(v_0\) at zero time is included. Note that \(v\) is a function of time, although this is not shown explicitly in the equation. To calculate the result it is easier to change the force equation to
where \(h = a/gm\). Separating variables and integrating gives
and using initial values
this produces
After substituting for \(h\) and rearranging, the velocity at any time \(t\) is,
To find the terminal velocity, the limit at long time is taken, which produces \(v_{term} = mg/ a\).
The distance from the ground at any time is found by integrating the velocity because \(v = dx/dt\),
If the initial height is \(x_0\) the result is
The time to reach the ground is found when \(x = 0\) and \(x_0\) is the height above ground taken as a negative number. The resulting transcendental expression has to be solved numerically for \(t\) which can be done by using the Newton-Raphson method to find the root of the equation. The figure below (8a) shows the velocity and distance plots vs time. The mass was taken to be \(75\) kg, approx an adult’s weight, and the height \(3000\) m (\(\approx 10000\) ft) typical of a skydiver’s jump height. The friction constant \(a = 13\) kg/s. The terminal velocity (left figure) is \(56.5\) m/s (or \(126.5\) mph) and the jump lasts for just under a minute before hitting the deck (right figure).
Fig 8a. Velocity and distance vs. time for a skydiver.
# Calculate time to hit the ground
v = lambda t : m*g/a - (m*g/a-v0)*np.exp(-a*t/m) # velocity
x = lambda t : m*g/a*t + (m/a**2)*(m*g - a*v0)*(np.exp(-a*t/m) -1) + x0 # distance
m = 75.0 # kg
g = 9.8 # acceleration due to gravity m/s/s
a = 13.0 # air resistance kg/s
v0 = 0.0 # initial vertical velocity
x0 = -3000.0 # initial height
vmax = m*g/a # terminal velpocity m/s
print('{:s}{:6.2f}{:s}{:6.2f}{:s}'.format('terminal velocity ',vmax,' m/s or ', vmax*2.237, ' mph'))
tm = 10.0 # initial guess
for i in range(3): # Newton-Raphson. V rapid convergence.
tm = tm - x(tm)/v(tm)
print('iteration ',i,tm)
pass
print('{:s}{:6.2f}{:s}'.format('time to reach ground ', tm, ' seconds'))
terminal velocity 56.54 m/s or 126.48 mph
iteration 0 68.07207783233369
iteration 1 58.83034251371377
iteration 2 58.83024022254945
time to reach ground 58.83 seconds
8.2 Rockets#
A rocket is fired vertically into the air with an initial velocity \(v_0\). To find the speed at any height where it has travelled a distance \(r\) the equation to use is,
where \(g\) is the acceleration due to gravity and \(R\) the radius of the earth and \(r\ge R\) since the rocket starts at the earth’s surface. Separating variables and integrating gives
The shape of curve of velocity \(v\) vs distance \(r\) shows that only when \(v_0^2 \ge 2gR\) does the rocket not return to earth. This condition defines the escape velocity \(v_e\) as does the calculation
where \(M\) is the mass of the earth and \(G\) the gravitational constant. The acceleration due to gravity \(g=GM/R^2\) and on earth \(v_e = 11.18\) km/s.
To find the time to travel a distance \(r\) the velocity equation is rearranged to isolate \(r\) and then \(dv/dt\) integrated. This is most easily done using the method shown in section 5.2. The result of this calculation is also shown in fig 8b.
Fig 8b. Rocket moving vertically and without air resistance. Left: velocity vs distance from earth’s surface (in units of earth’s mean radius \(R = 6372\) km). The initial rocket velocity on km/s is shown in the legend. The escape velocity curve is labelled with \(v_e\). Right. Time taken to cover distance \(r\) when the initial velocity is the escape velocity (black) and at a greater speed of \(12.5\) km/s (orange).