View source for Mason–Weaver equation

The '''Mason–Weaver equation''' (named after [[Max Mason]] and [[Warren Weaver]]) describes the [[sedimentation]] and [[diffusion]] of solutes under a uniform [[force]], usually a [[gravitation]]al field.<ref name="mason_1924" >{{cite journal | last = Mason | first = M |author2=Weaver W | year = 1924 | title = The Settling of Small Particles in a Fluid | journal = [[Physical Review]] | volume = 23 | issue = 3 | pages = 412–426 | doi = 10.1103/PhysRev.23.412 | bibcode=1924PhRv...23..412M}}</ref> Assuming that the [[gravitation]]al field is aligned in the ''z'' direction (Fig. 1), the Mason–Weaver equation may be written

:<math>
\frac{\partial c}{\partial t} = 
D \frac{\partial^{2}c}{\partial z^{2}} + 
sg \frac{\partial c}{\partial z}
</math>

where ''t'' is the time, ''c'' is the [[solute]] [[concentration]] (moles per unit length in the ''z''-direction), and the parameters ''D'', ''s'', and ''g'' represent the [[solute]] [[diffusion constant]], [[sedimentation coefficient]] and the (presumed constant) [[acceleration]] of [[gravitation|gravity]], respectively.

The Mason–Weaver equation is complemented by the [[boundary conditions]] 
:<math>
D \frac{\partial c}{\partial z} + s g c = 0
</math>
at the top and bottom of the cell, denoted as <math>z_a</math> and <math>z_b</math>, respectively (Fig. 1).  These [[boundary conditions]] correspond to the physical requirement that no [[solute]] pass through the top and bottom of the cell, i.e., that the [[flux]] there be zero.  The cell is assumed to be rectangular and aligned with 
the [[Cartesian coordinate system|Cartesian axes]] (Fig. 1), so that the net [[flux]] through the side walls is likewise 
zero.  Hence, the total amount of [[solute]] in the cell
:<math>
N_\text{tot} = \int_{z_b}^{z_a} \, dz \ c(z, t)
</math>
is conserved, i.e., <math>dN_\text{tot}/dt = 0</math>.


==Derivation of the Mason–Weaver equation==
[[Image:Mason Weaver cell.png|frame|left|Figure 1: Diagram of Mason–Weaver cell and Forces on Solute]]
A typical particle of [[mass]] ''m'' moving with vertical [[velocity]] ''v'' is acted upon by three [[force]]s (Fig. 1): the 
[[drag (physics)|drag force]] <math>f v</math>, the force of [[gravitation|gravity]] <math>m g</math> and the [[buoyancy|buoyant force]] <math>\rho V g</math>, where ''g'' is the [[acceleration]] of [[gravitation|gravity]], ''V'' is the [[solute]] particle volume and <math>\rho</math> is the [[solvent]] [[density]].  At [[mechanical equilibrium|equilibrium]] (typically reached in roughly 10 ns for [[molecule|molecular]] [[solutes]]), the 
particle attains a [[terminal velocity]] <math>v_\text{term}</math> where the three [[force]]s are balanced.  Since ''V'' equals the particle [[mass]] ''m'' times its [[partial specific volume]] <math>\bar{\nu}</math>, the [[mechanical equilibrium|equilibrium]] condition may be written as

:<math>
f v_\text{term} = m (1 - \bar{\nu} \rho) g \ \stackrel{\mathrm{def}}{=}\  m_b g
</math>

where <math>m_b</math> is the [[buoyant mass]].

We define the Mason–Weaver [[sedimentation coefficient]] <math>s \ \stackrel{\mathrm{def}}{=}\  m_b / f = v_\text{term}/g</math>.  Since the [[drag coefficient]] ''f'' is related to the [[diffusion constant]] ''D'' by the [[Einstein relation (kinetic theory)|Einstein relation]]

:<math>
D = \frac{k_B T} f
</math>,

the ratio of ''s'' and ''D'' equals

:<math>
\frac{s}{D} = \frac{m_b}{k_B T}
</math>

where <math>k_B</math> is the [[Boltzmann constant]] and ''T'' is the [[temperature]] in [[kelvin]]s.

The [[flux]] ''J'' at any point is given by

:<math>
J = -D \frac{\partial c}{\partial z} - v_\text{term} c 
  = -D \frac{\partial c}{\partial z} - s g c.
</math>

The first term describes the [[flux]] due to [[diffusion]] down a [[concentration]] gradient, whereas the second term 
describes the [[convective flux]] due to the average velocity <math>v_\text{term}</math> of the particles.  A positive net [[flux]] out of a small volume produces a negative change in the local [[concentration]] within that volume

:<math>
\frac{\partial c}{\partial t} = -\frac{\partial J}{\partial z}.
</math>

Substituting the equation for the [[flux]] ''J'' produces the Mason–Weaver equation

:<math>
\frac{\partial c}{\partial t} = 
D \frac{\partial^{2}c}{\partial z^{2}} + 
sg \frac{\partial c}{\partial z}.
</math>

==The dimensionless Mason–Weaver equation==

The parameters ''D'', ''s'' and ''g'' determine a length scale <math>z_0</math>

:<math>
z_0 \ \stackrel{\mathrm{def}}{=}\  \frac{D}{sg}
</math>

and a time scale <math>t_0</math>

:<math>
t_0 \ \stackrel{\mathrm{def}}{=}\  \frac{D}{s^{2}g^{2}}
</math>

Defining the [[dimensionless]] variables <math>\zeta \ \stackrel{\mathrm{def}}{=}\  z/z_0</math> and <math>\tau \ \stackrel{\mathrm{def}}{=}\  t/t_0</math>, the Mason–Weaver equation becomes

:<math>
\frac{\partial c}{\partial \tau} =
\frac{\partial^{2} c}{\partial \zeta^{2}} + 
\frac{\partial c}{\partial \zeta} 
</math>

subject to the [[boundary conditions]]

:<math>
\frac{\partial c}{\partial \zeta} + c = 0
</math>
at the top and bottom of the cell, <math>\zeta_a</math> and 
<math>\zeta_b</math>, respectively.

==Solution of the Mason–Weaver equation==

This partial differential equation may be solved by [[separation of variables]].  Defining <math>c(\zeta,\tau) \ \stackrel{\mathrm{def}}{=}\  e^{-\zeta/2} T(\tau) P(\zeta)</math>, we obtain two ordinary differential equations coupled by a constant <math>\beta</math>

:<math>
\frac{dT}{d \tau} + \beta T = 0
</math>

:<math>
\frac{d^2 P}{d \zeta^2} + 
\left[ \beta - \frac 1 4 \right] P = 0
</math>

where acceptable values of <math>\beta</math> are defined by the [[boundary conditions]]

:<math>
\frac{dP}{d\zeta} + \frac{1}{2} P = 0
</math>

at the upper and lower boundaries, <math>\zeta_a</math> and <math>\zeta_b</math>, respectively.  Since the ''T'' equation has the solution <math>T(\tau) = T_0 e^{-\beta \tau}</math>, where <math>T_0</math> is a constant, the Mason–Weaver equation is reduced to solving for the function <math>P(\zeta)</math>.

The [[ordinary differential equation]] for ''P'' and its [[boundary conditions]] satisfy the criteria
for a [[Sturm–Liouville theory|Sturm–Liouville problem]], from which several conclusions follow.  '''First''', there is a discrete set of [[orthonormal]] [[eigenfunction]]s 
<math>P_k(\zeta)</math>  that satisfy the [[ordinary differential equation]] and [[boundary conditions]].  '''Second''', the corresponding [[eigenvalue]]s <math>\beta_k</math> are real, bounded below by a lowest 
[[eigenvalue]] <math>\beta_0</math> and grow asymptotically like <math>k^{2}</math> where the nonnegative integer ''k'' is the rank of the [[eigenvalue]].  (In our case, the lowest eigenvalue is zero, corresponding to the equilibrium solution.)  '''Third''', the [[eigenfunction]]s form a complete set; any solution for <math>c(\zeta, \tau)</math> can be expressed as a weighted sum of the [[eigenfunction]]s

:<math>
c(\zeta, \tau) = 
\sum_{k=0}^\infty c_k P_k(\zeta) e^{-\beta_k\tau}
</math>

where <math>c_k</math> are constant coefficients determined from the initial distribution <math>c(\zeta, \tau=0)</math>

:<math>
c_k = 
\int_{\zeta_a}^{\zeta_b} d\zeta \ 
c(\zeta, \tau=0) e^{\zeta/2} P_k(\zeta)
</math>

At equilibrium, <math>\beta=0</math> (by definition) and the equilibrium concentration distribution is

:<math>
e^{-\zeta/2} P_0(\zeta) = B e^{-\zeta} = B e^{-m_b gz/k_B T}
</math>

which agrees with the [[Boltzmann distribution]].  The <math>P_0(\zeta)</math> function satisfies the [[ordinary differential equation]] and [[boundary conditions]] at all values of <math>\zeta</math> (as may be verified by substitution), and the constant ''B'' may be determined from the total amount of [[solute]]

:<math>
B = N_\text{tot} \left( \frac{sg} D \right)
\left( \frac 1 {e^{-\zeta_b} - e^{-\zeta_a}} \right)
</math>

To find the non-equilibrium values of the [[eigenvalue]]s <math>\beta_k</math>, we proceed as follows.  The P equation has the form of a simple [[harmonic oscillator]] with solutions <math>P(\zeta) = e^{i\omega_k\zeta}</math> where

:<math>
\omega_k = \pm \sqrt{\beta_k - \frac 1 4}
</math>

Depending on the value of <math>\beta_k</math>, <math>\omega_k</math> is either purely real (<math>\beta_k\geq\frac 1 4</math>) or purely imaginary (<math>\beta_k < \frac 1 4</math>).  Only one purely imaginary solution can satisfy the [[boundary conditions]], namely, the equilibrium solution.  Hence, the non-equilibrium [[eigenfunctions]] can be written as

:<math>
P(\zeta) = A \cos{\omega_k \zeta} + B \sin{\omega_k \zeta}
</math>

where ''A'' and ''B'' are constants and <math>\omega</math> is real and strictly positive.

By introducing the oscillator [[amplitude]] <math>\rho</math> and [[phase (waves)|phase]] <math>\varphi</math> as new variables,

:<math>
u \ \stackrel{\mathrm{def}}{=}\  \rho \sin(\varphi) \ \stackrel{\mathrm{def}}{=}\  P
</math>

:<math>
v \ \stackrel{\mathrm{def}}{=}\  \rho \cos(\varphi) \ \stackrel{\mathrm{def}}{=}\  - \frac 1 \omega
\left( \frac{dP}{d\zeta} \right)
</math>

:<math>
\rho \ \stackrel{\mathrm{def}}{=}\  u^2 + v^2
</math>

:<math>
\tan(\varphi) \ \stackrel{\mathrm{def}}{=}\  v / u
</math>

the second-order equation for ''P'' is factored into two simple first-order equations

:<math>
\frac{d\rho}{d\zeta} = 0
</math>

:<math>
\frac{d\varphi}{d\zeta} = \omega
</math>

Remarkably, the transformed [[boundary conditions]] are independent of <math>\rho</math> and the endpoints <math>\zeta_a</math> and <math>\zeta_b</math>

:<math>
\tan(\varphi_a) = 
\tan(\varphi_b) = \frac 1 {2\omega_k}
</math>

Therefore, we obtain an equation

:<math>
\varphi_a - \varphi_b + k\pi = k\pi = 
\int_{\zeta_b}^{\zeta_a} d\zeta \ \frac{d\varphi}{d\zeta} = 
\omega_k (\zeta_a - \zeta_b) 
</math>

giving an exact solution for the frequencies <math>\omega_k</math>

:<math>
\omega_k = \frac{k\pi}{\zeta_a - \zeta_b}
</math>

The eigenfrequencies <math>\omega_k</math> are positive as required, since <math>\zeta_a > \zeta_b</math>, and comprise the set of [[harmonic]]s of the [[fundamental frequency]] <math>\omega_1 \ \stackrel{\mathrm{def}}{=}\  \pi/(\zeta_a - \zeta_b)</math>.   Finally, the [[eigenvalue]]s <math>\beta_k</math> can be derived from <math>\omega_k</math>

:<math>
\beta_k = \omega_k^2 + \frac 1 4
</math>

Taken together, the non-equilibrium components of the solution correspond to a [[Fourier series]] decomposition of the initial concentration distribution <math>c(\zeta, \tau=0)</math>
multiplied by the [[weight function|weighting function]] <math>e^{\zeta/2}</math>.  Each Fourier component decays independently as <math>e^{-\beta_k\tau}</math>, where <math>\beta_k</math> is given above in terms of the [[Fourier series]] frequencies <math>\omega_k</math>.

==See also==
* [[Lamm equation]]
* The Archibald approach, and a simpler presentation of the basic physics of the Mason–Weaver equation than the original.<ref>{{cite journal | last=Archibald | first=William J. | title=The Process of Diffusion in a Centrifugal Field of Force | journal=Physical Review | publisher=American Physical Society (APS) | volume=53 | issue=9 | date=1938-05-01 | issn=0031-899X | doi=10.1103/physrev.53.746 | pages=746–752| bibcode=1938PhRv...53..746A }}</ref>

==References==

{{reflist|1}}

{{DEFAULTSORT:Mason-Weaver equation}}
[[Category:Laboratory techniques]]
[[Category:Partial differential equations]]