You are currently browsing the category archive for the ‘Uncategorized’ category.

Random fluctuations and noise are ubiquitous within nature and will generally occur on some length scale. Consider long linear DNA strands suspended within a solution at some temperature T and salt concentration Or in-vivo, compactified within a living cell, within a biological environment with average ambient temperature T. Thermal fluctuations will cause the DNA ribbon to fluctuate weakly about its average geometry. In addition, the thermal fluctuations will also result in the spontaneous opening and closing of the bonds or bridges spanning the edges of the DNA ribbon. The DNA ribbon is therefore not a static structure but fluctuates randomly, essentially becoming a stochastic or random geometry.} The aim therefore, will be to describe the fluctuating DNA ribbon mathematically, in terms of random geometrical and topological terms. First, some statistical mechanical issues are considered. The following proposition follows the arguments in Huang. The following assumptions are made:

A DNA suspension within a liquid at temperature and some given salt concentration, or in=vivo within a cell. It can be considered within the canonical ensemble as a subsystem within a heat bath at fixed temperature. The DNA strand of length \ell and sugar-phosphate backbone edges  h,\bar{h}) has N bonds or weak hydrogen bridges across the ribbon binding the A-T and G-C nucleic acids. At temperature T the ribbon weakly fluctuates and the bridges open and close. A open bond can be assigned an energy E and an open bond has energy  \mathcal{E}=0. A link can be open only if all the other links along the ribbon are already open. The possible states are labeled by the number of open bridges or links  \alpha=0,1,2,...,N so that the energy with  \alpha open links is  E_{\alpha}=\alpha{E} The at low temperatures  E/T\gg 1, and so there are few open links in that  \langle\alpha\rangle\sim \exp(-E/T). At high temperatures  E/T\ll 1, and so most of the links are open in that  \langle\alpha\rangle\sim \alpha
The partition function is

 Z_{n}=\sum_{\alpha=0}^{n}\exp[-\alpha{E}/T]\equiv\frac{1-\exp[-\beta(n+1)E}{1-exp[-\beta E]}

where \beta=1/T is the inverse temperature. The expected number of open links at temperature T is then given by the thermodynamic relation

 \langle\alpha\rangle=-\frac{1}{E}\frac{\partial \ln[Z_{n}]}{\partial\beta} so that    \langle\alpha\rangle = \frac{\sum_{\alpha=0}^{N}\alpha n\exp[-nE/T]}{\sum_{\alpha=0}^{n}\alpha \exp[-nE/T]}=\frac{\exp[-E/T]}{1 - \exp[-E/T]}-\frac{(n+1)\exp[-(n+1)E/T}{1-\exp[-(n+1)E/T]}

Then at low temperatures  E/T\gg 1, and so there are few open links in that  \langle\alpha\rangle\sim \exp(-{E}/T). At high temperatures E/T\ll 1, and so most of the links are open in that \langle\alpha\rangle\sim \alpha

This suggests that DNA will exist and function only within certain temperature ranges. The molecule dissociates at high temperatures. (Which also implies that DNA-based life cannot arise and evolve in high-temperature conditions in the universe.) Superhelicity also plays a role in that the free energy  \delta G associated with supercoiling is of the quadratic form  \Delta G \sim K(Lk-Lk_{0})^{2}, where K is a proporitonality constant. This quadratic relation has been experimentally shown (refs.) DNA molecules that differ only in \mathbf{LK} are  topoisomers. The concentration C[Lk,T] of a given topoisomer with linking number  Lk at temperature T, is the Boltzmann distribution

 C[\mathcal{LK},T]=\frac{1}{Z}\exp[[-\Delta G/RT]\sim \frac{1}{Z}\exp[-K(\mathcal{Lk}-\mathcal{L}_{r})^{2}/R  which is a normal or Gaussian distribution centred on \mathbf{LR}, the relaxed linking number, and with standard deviation \sim \sqrt{RT/2K}. Proteins which regulate the supercoiling and maintain this distribution are known as topoisomerases. For high temperatures C[\mathcal{LK},T]\rightarrow 0.

Taking into account thermodynamic, mechanical and biophysical properties of fluctuating DNA molecules within a thermal bath is complex. However, the aim of this section is to focus on a tractable stochastic mathematical description of a fluctuating DNA ribbon, in terms of a ‘stochastic geometry and topology’ which arises from contact with a thermal bath, or more precisely, a thermal noise bath.
\item This introduces analogies with Brownian motion and turbulence, which are studied and described on length scales as stochastic processes and where the deeper underlying statistical mechanics can be ignored. Indeed, the stochastic fluctuations of the DNA double-helix structure when immersed in a thermal bath whether in-vitro or in vivo, is a form of Brownian motion, arising from bombardment of the double helix by a ‘warm sea’ of smaller surrounding molecules.

For proteins, which consist of mobile helices, strands and sheets, Brownian bombardment from water molecules will also play an important role in the folding of the protein. Simple Brownian motion of a particle in a thermal bath at temperature T is described by a linear SDE of the for

d\widehat{X}(t)=dX(t)dt + \rho d\widehat{\mathcal{W}}(t)

where X(t) is the position of the Brownian particle. The random field \widehat{\mathcal{W}}(t) is a thermal white noise which described the effect of the thermal Brownian bombardment of the particle by the smaller water molecules. Generally

\mathbb{E}\langle\widehat{\mathcal{W}}(t)\rangle=0 and by the

equipartition theory  \mathbb{E}\langle \widehat{\mathcal{W}}(t)\otimes\widehat{\mathcal{W}}(s)\rangle \sim kT \delta(t-s) , where \mathbb{E}\langle...\rangle is the statistical expectation, k is Boltzmann constant. The solution is the well-defined stochastic Ito integral

 X(t)=\int_{0}^{t}X(u)du + \rho\int_{0}^{t}d\widehat {\mathcal{W}}(u)= \int_{0}^{t}X(u)du + \rho\int_{0}^{t}\widehat{\psi}(u)du

Suppose by analogy, we consider the tangent vectors X_{i}(x(s)) and X_{j}(y(s)) along the edges of the DNA ribbon. If the DNA ribbon undergoes a Brownian bombardment from smaller molecules, these vectors will be randomly perturbed as the ribbon fluctuates randomly about its average geometry/topology.

d\widehat{X}_{i}(x(s))=dX_{i}(x(s))+d\widehat{\mathcal{F}}_{i}(x(s))

d\widehat{X}_{j}(y(s))=dX_{j}(y(s))+d\widehat{\mathcal{F}}_{j}(y(s)

where \widehat{\mathcal{F}}_{i}(x(s)) is a Gaussian random vector field defined at all points in a domain \mathbb{D} containing the ribbon. GRVFs are discussed in the Appendix. The linking number will not be affected by thermal fluctuations of the DNA ribbon but the linking number $\mathcal{LK}=\mathcal{TW}+\mathcal{WR}$ will be randomly partitioned into twist and writhe which will randomly change. The total remains constant however. This can be made more precise as follows:

Let R(h,\bar{h},\alpha,\beta,\epsilon,p:Lk,Sh\subset\mathbb{D} be a DNA ribbon with linking number Lk and superhelicity Sh and let
X_{i}(x) and X_{j}(y) be tangent vectors at x^{i}\in h and y^{j}\in \bar{h}. The Wilson lines along the helices or ribbon edges are

  W[h]=exp\bigg(i\int_{h}X_{i}(x)dx^{i}\bigg)=\exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2}Lk(h,\bar{h}))

 W[\bar{h}]=exp\bigg(i\int_{\bar{h}}X_{i}(x)dx^{i}\bigg)=exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h}))

If the domain \mathbb{D} is a thermal bath at some temperature T then the tangent vectors can be interpreted as random Gaussian vector fields \widehat{X}_{i}(x) and \widehat{X}_{j}(y) such that

  \widehat{X}_{i}(x)= X_{i}(x) + \widehat{\mathcal{F}}_{i}(x) \widehat{X}_{j}(y)= X_{j}(y) + \widehat{\mathcal{F}}_{j}(y)

where the GRVF \widehat{\mathcal{F}}_{i}(x) obeys the statistics
 \mathbf{E}(\widehat{\mathcal{F}}_{i}(x))=0,~~~~\mathbf{E}(\widehat{\mathcal{F}}_{i}(x)\otimes \widehat{\mathcal{F}}_{j}(y))=C_{ij}(x,y;\zeta)  and where C_{ij}(x,y;\zeta) is a 2-point function or covariance regulated at x=y and \zeta is a correlation length such that C_{ij}(x-y)\rightarrow 0 for |x-y|\gg \zeta.
(Appendix A.)

The stochastic Wilson lines along the helices are then
 W[h]=exp\bigg(i\int_{h}\widehat{X}_{i}(x)dx^{i}\bigg)=exp\bigg(i\int_{h}(\widehat{X}_{i}(x) +\widehat{\mathcal{F}}_{i}(x))dx^{i}\bigg)=exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h}))exp\bigg(i\int_{h}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)

W[\bar{h}])=exp\bigg(i\int_{h}\widehat{X}_{i}(x)dx^{i}\bigg)=exp\bigg(i\int_{h}(\widehat{X}_{i}(x) +\widehat{\mathcal{F}}_{i}(x))dx^{i}\bigg)=exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h}))exp\bigg(i\int_{\bar{H}}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)

The stochastic expectations of the Wilson lines are
 \mathbf{W}[h]=\mathbf{E}(W[h])=\mathbf{E}\bigg(exp\bigg(i\int_{h}\widehat{X}_{i}(x)dx^{i}\bigg)\bigg)=exp\bigg(i\int_{h}(\widehat{X}_{i}(x) +\widehat{\mathcal{F}}_{i}(x))dx^{i}\bigg) =  exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h}))\mathbf{E}\bigg(exp\bigg(i\int_{h}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)\bigg)

\ \mathbf{W}[\bar{h}]=\mathbf{E}(W[H])=\mathbf{E}\bigg(exp\bigg(i\int_{h}\widehat{X}_{i}(x)dx^{i}\bigg)\bigg)= exp\bigg(i\int_{h}(\widehat{X}_{i}(x) +\widehat{\mathcal{F}}_{i}(x))dx^{i}\bigg)=exp(\tfrac{1}{2}i(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h}))\mathbf{E}\bigg(exp\bigg(i\int_{\bar{h}}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)\bigg)

 \mathbf{W}[h]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h})\mathbf{\mathcal{W}}(h)

 \mathbf{W}[\bar{h}]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h})\mathbf{\mathcal{W}}(\bar{h}

The stochastic integral

 \mathbf{E}(exp(i\int_{h}\widehat{\mathcal{F}}_{i}(x)dx^{i})

can be evaluated by a cluster integral decomposition, exploiting the Gaussian property of the random field \widehat{\mathcal{F}}_{i}(x). This is briefly stated and proved in the following theorem. Let \widehat{\mathcal{F}}_{i}(x) be a random Gaussian vector field defined for all  x^{i}\in\mathbb{D} with respect to a probability space (\Omega, F,\mathbf{P}) and having a 2-point function C_{ij}(x,y;\zeta). Then the expectations of the stochastic Wilson loops along h and \bar{h} are then

\mathbf{E}\bigg(exp\bigg(i\int_{h}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)\bigg)=exp\bigg(-\int_{h}\int_{\bar{h}}C_{ij}(x,y;\zeta)dx^{i}dy^{j}\bigg)  latex  \mathbf{\mathcal{W}}[\bar{h}]=\mathbf{E}\bigg(exp\bigg(i\int_{\bar{h}}\widehat{\mathcal{F}}_{i}(x)dx^{i}\bigg)\bigg)= exp\bigg(-\int_{h}\int_{\bar{h}}C_{ij}(x,y;\zeta)dx^{i}dy^{j}\bigg)$

\mathbf{W}[h]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h})exp\bigg(-\int_{h}\int_{h}C_{ij}(x,y;\zeta)dx^{i}dy^{j}\bigg)

\mathbf{W}[\bar{h}]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h})\exp\bigg(-\int_{h}\int_{\bar{h}}C_{ij}(x,y;\zeta)dx^{i}dy^{j}\bigg)

 \mathbf{W}[h]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2}Lk(h,\bar{h})\exp (-Q_{ij})

 \mathbf{W}[\bar{h}]=\mathbf{E}(\widehat{W}[h])=exp(\tfrac{1}{2}(\alpha^{2}+p^{2})^{1/2} Lk(h,\bar{h})\exp(-Q_{ij})

where

Q_{ij}(x,y;\zeta)=\int_{h}\int_{\bar{h}} C_{ij}(x,y;\zeta)dx^{i}dy^{j}. If Q_{ij}(x,y;\zeta)\rightarrow 0 for |x-y|\gg \zeta then exp(-Q_{ij})\rightarrow 1 as

Deep within your molecular biology are the paradoxical molecules called “Topoisomerase II”, which are potentially very dangerous in that they continuously and literally splice up your entire genome––and yet no life on Earth could ever exist or function without them. They temporarily induce immense damage to all your DNA…and then repair it perfectly. How this molecule functions on a near-atomic scale remains mostly a mystery but ALL life on Earth requires this molecule to exist and function, as do cancer cells. And it once also seemed that the key to finally curing cancer was to shut down this molecule so it cannot perform its function on DNA within cancerous cells.

In order for you to exist from microsecond to microsecond your entire genome–the sum total of your DNA–must be continuously maintained and regulated and repaired. This is done through an arsernal of enzymes, most crucially Topoisomerase II. The function of this enzyme/protein is to perform a crucial mathematical function: to remove knots and tangles and relieve overwinding/torsional stresses (positive supercoiling)  within DNA. If DNA has knots or tangles, then tracking enzymes like RNA polymerase (responsible for protein transcription and renewing all your proteins) and DNA polymerase (which is initially responsible for DNA replication and thus cell division) could never function. They move along DNA like beads on a string and cannot proceed if they enounter knots or tangles.

However Topoisomerase II removes the DNA knots by performing a very dangerous “double-strand break” on your DNA.: Much like Alexander the Great’s solution to the Gordian Knot, they cut it open like a sword. The cut knot is ‘unknotted’ or removed and then the strand is resealed. However, it does more than that: it ensures DNA remains “underwound” or “negatively supercoiled” with a specific ‘Gaussian linking number’, (A formula first written down by Gauss in 1820.)). The DNA of all life on Earth shares this property.This molecule seems to know all the mathematical tricks from knot theory textbooks. However, at any instant, these molecules are splicing up your entire genome/DNA like scissors cutting string into billions of pieces…yet at the same time they reseal it all back together with all the knots removed…and so you continue to function and exist.

The class of cancer drugs called “Doxyrubicins” were designed to trap or ‘freeze’ Topoisomerase II so that the cancer DNA would remain knotted and tangled. Hence, it would not be able to replicate itself and the cancer cells would die. However, the cancer cells which survived the first rounds of treatment were able to unknot themselves with fewer Topoisomerase enzymes. Also, new genes can be turned on which enables the cells to expell the drug. Neverthleless, while not a cure, they have proven themselves to be effective drugs.”

……………………………………………………………………………………………………..(Follows closely: 738–748 Nucleic Acids Research, 2009, Vol. 37, No. 3) Topoisomerase II is an essential enzyme that is required for virtually every process that requires movement of DNA within the nucleus or the openingof the double helix. This enzyme helps to regulate DNA under- and overwinding and removes knots and tangles from the genetic material. In order to carry out its critical physiological functions, topoisomerase II generates transient double-stranded breaks in DNA. Consequently, while necessary for cell survival, the enzyme also has the capacity to fragment the genome. The DNA cleavage/ligation reaction of topoisomerase II is the target for some of the most successful anticancer drugs currently in clinical use. However, this same reaction also is believed to trigger chromosomal translocations that are associated with specific types of leukemia. This article will familiarize the reader with the DNA cleavage/ligation reaction of topoisomerase II and other aspects of its catalytic cycle. In addition, it will discuss the interaction of the enzyme with anticancer drugs and the mechanisms by which these agents increase levels of topoisomerase II-generated DNA strand breaks.

A number of enzymes that catalyze essential physiological processes also have the capacity to damage the genome during the course of their normal activities. For example, while the cell requires DNA polymerases to copy the genetic material, these enzymes insert an incorrect base approximately every 107 nt (1). Consequently, in the absence of mismatch repair pathways, human DNA
polymerases would generate several hundred mutations every round of cell division. Furthermore, while DNA glycosylases initiate base-excision repair pathways, these enzymes can convert innocuous lesions to abasic sites
with far greater mutagenic potential (2). Finally, while cytochrome P450 enzymes play critical roles in detoxification pathways, they sometimes convert inert xenobiotic chemicals to compounds with mutagenic properties (3). Of all the enzymes required to sustain cellular growth, topoisomerase II is one of the most dangerous (4–8). As discussed below, this enzyme unwinds, unknots and
untangles the genetic material by generating transient double-stranded breaks in DNA (8–12). Although the cell cannot survive without topoisomerase II, the strand breaks that the enzyme generates have the potential to trigger cell death pathways or chromosomal translocations.

This article focuses on the mechanism by which topoisomerase
II cleaves the genetic material, the ability to exploit this reaction for the chemotherapeutic treatment of human cancers and the role of this reaction in triggering specific types of leukemia.
DNA TOPOLOGY
The existence of topoisomerases is necessitated by the structure of the double helix. Each human cell contains 2m of DNA that are compacted into a nucleus that is 10 mm in diameter (14,15). Because the genetic material is anchored to the chromosome scaffold and the two strands of the double helix are plectonemically coiled, accessing the genome is a complex topological challenge (11,12,16–18). Topological properties of DNA are those that can only be changed when the double helix is broken (12). Two aspects of DNA topology significantly affect nuclear processes. The first deals with topological relationships between the two strands of the double helix.

In all living systems, from bacteria to humans, DNA is globally underwound
(i.e. negatively supercoiled) by 6% (12,19–21). This is important because duplex DNA is merely the storage form for the genetic information. In order to
replicate or express, DNA must be separated. Since global underwinding of
the genome imparts increased single-stranded character to the double helix, negative supercoiling greatly facilitates strand separation (12,16–18). While negative supercoiling promotes many nucleic acid processes, DNA overwinding (i.e. positive supercoiling) inhibits them. The linear movement of tracking enzymes, such as helicases and polymerases, compresses the turns of the double helix into a shorter region (Figure 1) (12,19–21). Consequently, the double helix becomes increasingly overwound ahead of tracking systems. The positive supercoiling that results makes it more difficult to open the two strands of the double helix and ultimately blocks essential nucleic acid processes (10,12,16–18). The second aspect of DNA topology deals with relationships between separate DNA segments. Intramolecular knots (formed within the same DNA molecule) are
generated during recombination, and intermolecular tangles (formed between daughter DNA molecules) are produced during replication (Figure 1) (8,10,12,17).

DNA knots block essential nucleic acid processes because they make it impossible to separate the two strands of the double helix. Moreover, tangled DNA molecules cannot be segregated during mitosis or meiosis (8,10,12,17). Consequently, DNA knots and tangles can be lethal to cells if they are not resolved.
DNA TOPOISOMERASES
The topological state of the genetic material is regulated by enzymes known as topoisomerases (8,10,11,22,23). Topoisomerases are required for the survival of all organisms and alter DNA topology by generating transient breaks in the double helix (8,10,11,22,23). There are two major classes of topoisomerases, type I and type II, that are distinguished by the number of DNA strands that they cleave and the mechanism by which they alter the topological properties of the genetic material (8,10,11,22,23). Eukaryotic type I topoisomerases are monomeric
enzymes that require no high-energy cofactor (11,22,24). Type I enzymes are organized into two subclasses: type IA and type IB. These enzymes alter topology by creating transient single-stranded breaks in the DNA, followed by passage of the opposite intact strand through the break (type IA) or by controlled rotation of the helix around the break (type IB) (11,22,24). Type IA topoisomerases need divalent metal ions for DNA scission and attach covalently to the 50-terminal phosphate of the DNA (11,22,24). In contrast, type IB enzymes do not require divalent metal ions and covalently link to the 30-terminal phosphate (11,22,24).

As a result of their reaction mechanism, type I topoisomerases can modulate
DNA under- and overwinding, but cannot remove knots or tangles from duplex DNA. A number of excellent review articles on type I topoisomerases have
appeared recently (22,24,25). Consequently, these enzymes will not be discussed further in this article. Eukaryotic type II topoisomerases function as homodimers and require divalent metal ions and ATP for complete catalytic activity (5,8,26–28). These enzymes interconvert different topological forms of DNA by a ‘double-stranded DNA passage reaction’ that can be separated into a number of discrete steps (5,8,26–28). Briefly, type II topoisomerases (i) bind two separate segments of DNA, (ii) create a double-stranded break in one of the segments, (iii) translocate the second DNA segment through the cleaved nucleic acid ‘gate’, (iv) rejoin (i.e. ligate) the cleaved DNA, (v) release the translocated segment through a gate in the protein and (vi) close the protein gate and regain the ability to start a new round of catalysis (5,26–34). Because of their double-stranded DNA passage mechanism, type II topoisomerases can modulate DNA supercoiling and also can remove DNAknots and tangles.

TOPOISOMERASE II
Lower eukaryotes and invertebrates encode only a single type II topoisomerase, topoisomerase II (35–38). In contrast, vertebrate species encode two closely related Nuclear processes induce changes in DNA topology. DNA replication is used as an example. Although chromosomal DNA is globally underwound in all cells, the movement of DNA tracking systems generates positive supercoils. As shown in (A) chromosomal DNA ends are tethered to membranes or the chromosome scaffold (represented by the red spheres) and are unable to rotate. Therefore, the linear movement of tracking systems (such as the replication machinery represented by the yellow bars) through the immobilized double helix
compresses the turns into a shorter segment of the genetic material and
induces acute overwinding (i.e. positive supercoiling) ahead of the fork
(B). In addition, the compensatory underwinding (i.e. negative supercoiling)
behind the replication machinery allows some of the torsional stress that accumulates in the prereplicated DNA to be translated to the newly replicated daughter molecules in the form of precatenanes (C). If these precatenanes are not resolved, they ultimately lead to the formation of intertwined (i.e. tangled) duplex daughter chromosomes. Adapted from ref. 10. Nucleic Acids Research, 2009, Vol. 37, No. 3 739 isoforms of the enzyme, topoisomerase IIa and topoisomerase IIb. These isoforms differ in their protomer molecular masses (170 versus 180 kDa, respectively) and are encoded by separate genes (8,10,22,28,39–46).

Topoisomerase IIa and topoisomerase IIb display a high degree (70%) of amino acid sequence identity and similar enzymological characteristics. One notable difference between the two isoforms is that topoisomerase IIa relaxes (i.e. removes) positive superhelical twists 10 times faster than it does negative in vitro, while the b isoform is unable to distinguish the geometry of DNA supercoils during DNA relaxation (47). Topoisomerase IIa and topoisomerase IIb have distinct patterns of expression and separate cellular functions. Topoisomerase IIa is essential for the survival of proliferating cells, and protein levels rise dramatically during periods of cell growth (48–51). The enzyme is further
regulated over the cell cycle, with protein concentrations peaking in G2/M (50,52,53). Topoisomerase IIa is associated with replication forks and remains tightly bound to chromosomes during mitosis (9,51,54–56). Thus, it is believed to be the isoform that functions in growth-related processes, such as DNA replication and chromosome segregation (10,51). Topoisomerase IIb is dispensable at the cellular level but appears to be required for proper neural development (57–59).

Expression of topoisomerase IIb is independent of proliferative status and cell cycle, and the enzyme dissociates from chromosomes during mitosis (54,60,61). Topoisomerase IIb cannot compensate for the loss of topoisomerase IIa in mammalian cells, suggesting that these two isoforms do not play redundant roles in replicative processes (51,60,62,63). Although the physiological functions of topoisomerase IIb have yet to be defined, recent evidence indicates involvement in the transcription of hormonally or developmentally regulated genes (63,64). Much of what we understand regarding the mechanism of action of type II enzymes comes from experiments with topoisomerase II from species that express only a single form of the protein. Consequently, eukaryotic type II
topoisomerases will be referred to collectively as topoisomerase
II, unless the properties being discussed are specific to either the a or b isoform.
TOPOISOMERASE II-MEDIATED DNA CLEAVAGE
AND LIGATION
The ability of topoisomerase II to cleave and ligate DNA is central to all of its catalytic functions (5,8,11,27). All topoisomerases utilize active site tyrosyl residues to mediate DNA cleavage and ligation. Since type II enzymes cut both strands of the double helix, each protomer subunit contains one of these residues (Tyr805 and Tyr821 in human topoisomerase IIa and topoisomerase IIb,
respectively). Topoisomerase II initiates DNA cleavage by the nucleophilic
attack of the active site tyrosine on the phosphate of the nucleic acid backbone (Figure 2) (11,23,26,27). The resulting transesterification reaction results in the
formation of a covalent phosphotyrosyl bond that links the protein to the newly generated 50-terminus of the DNA chain.

It also generates a 30-hydroxyl moiety on the opposite terminus of the cleaved strand. The scissile bonds on the two strands of the double helix are staggered and are located across the major groove from one another. Thus, topoisomerase II generates cleaved DNA molecules with four-base 50-single-stranded cohesive ends, each of which is covalently linked to a separate protomer subunit of the enzyme (65–67). The covalent enzyme–DNA linkage plays two important roles in the topoisomerase II reaction mechanism. First, it conserves the bond energy of the sugar-phosphate DNA backbone. Second, because it does not allow the cleaved DNA chain to dissociate from the enzyme, the protein– DNA linkage maintains the integrity of the genetic material during the cleavage event (11,23,26,27). The covalent topoisomerase II-cleaved DNA reaction intermediate is referred to as the ‘cleavage complex’ and is critical for the pharmacological activities of the enzyme, which are discussed later in this article.

Although topoisomerase II acts globally, it cleaves DNA at preferred sites (68). The consensus sequence for cleavage is weak, and many sites of action do not
conform to it (68). Ultimately, the mechanism by which topoisomerase II selects DNA sites is not apparent, and it is nearly impossible to predict de novo whether a given DNA sequence will support scission. Most likely, the specificity of topoisomerase II-mediated cleavage is determined by the local structure, flexibility, or malleability of the DNA that accompanies the sequence, as opposed
to a direct recognition of the bases that comprise that sequence (69). Beyond the nucleophilic attack of the active site tyrosine on the DNA backbone, the details of topoisomerase II mediated DNA cleavage are not well defined.

Pycnonuclear reactions occur in high density matter when nuclei are frozen in lattice structures. Such matter characterizes the cores of white dwarfs, and the crusts of accreting neutron stars. While pycnonuclear reactions in white dwarf matter are dominated by 12C and 16O induced fusion processes, pycnonuclear reactions in neutron star crust matter are dominated by fusion between very neutron rich carbon, oxygen, and neon isotopes at the limits of stability. In a post X-ray burst neutron star binary system ashes from the burst are forced deep into the crust by the weight of freshly accreted matter. With increasing density the remnant abundance distribution is changed by a host of electron capture reactions and neutron emissions, leading to the effective processing of the ashes into very neutron rich nuclei in the carbon to magnesium range. Near nuclear matter densities the nuclei are forced into a lattice, surrounded by a neutralizing degenerate electron gas. Though held in a solid lattice, reactions are still possible due to the reduction of Coulomb potential by screening and the overlap of nuclear wave functions between lattice sites. This type of reaction is highly density dependant, and unlike the more familiar thermonuclear counterpart, is not strictly temperature dependent

A famous result of Witten is a derivation of the Gaussian linking number from an Abelian Chern-Simons theory and the Jones polynomials from the non-abelian theory, thus providing a bridge between quantum theory and knot theory. The Chern-Simons theory is a topological quantum field theory:  the Hamiltonian  H has only zero eigenstates and the Hilbert space is finite. For a 3-manifold  \mathbb{M}=\mathbb{R}^{3} the VEV of the Wilson loop is

\langle W[K]\rangle = \int DA\exp(iS_{cs})W[K]   = {\displaystyle \int} \mathcal{D}A\exp \bigg(\frac{ik}{4\pi}{\displaystyle \int_{\mathbb{M}}} tr({A} \wedge  d{A}+\frac{2}{3}{A} \wedge  A \wedge A)\bigg) exp\bigg(i\int A\bigg)

where

S_{CS}=(ik/4\pi g){\displaystyle \int_{\mathbb{M}}}Tr({A}\wedge d{A}+\frac{2}{3} {A} \wedge A \wedge A\wedge A)

for k \in \mathbb{Z} and g=1, and  A is an  SU(2) connection or gauge field. The partition function Z(\mathbb{M}) is itself an invariant (the Witten invariant) of the 3-manifold and can be calculated analytically. Witten (refs) demonstrated that this connects with the Jones polynomial of knot theory and the simple non-Abelian CS theory reproduces the Gaussian linking number

  \langle W[K] \rangle  =  \int \mathcal{D}A exp \bigg (\frac{ik}{4\pi}{\displaystyle \int_{\mathbb{M}} }tr({A}\wedge d{A})\bigg) exp\bigg(i {\displaystyle \int} A\bigg)=exp(\mathbf{\mathcal{L}k}[K, K'])

The Gaussian linking number \mathcal{GL}(K,\ \bar{K}) of the curves (K,\bar{K}) describes the number of times one curve intersects the surface bounded by the other.

  \mathcal{GL}(K,\bar{K})= \mathcal{LK} =\frac{1}{4\pi}{\displaystyle \int_{K}}dx^{i} {\displaystyle \int_{\bar{K}}}dy^{j}\epsilon_{ijk}\frac{(x-y)^{k}}{|x-y|^{3}}

and for closed curves or knots  (C,\bar{C}) forming a link

 \mathcal{GL} (K,\bar{K})=\mathcal{LK}\frac{1}{4\pi}{\displaystyle \oint_{\mathfrak{S}}}dx^{i}{\displaystyle \oint_{\bar{C}}}dy^{j}\epsilon_{ijk}\frac{(x-y)^{k}}{|x-y|^{3}}

Equivalently in terms of the paremetrization s \in[0,1],

  \mathcal{GL}(K,\bar{K})=\mathcal{LK}= \mathcal{TW}_{o} =  \frac{1}{4\pi}{\displaystyle \int_{K}} {\displaystyle \int_{\bar{K}}}ds  ds'   [\frac{d\vec{y}(s')}{ds'} \wedge \frac{d\vec{y}(s')}{ds'}]   \frac{(\vec{x}(s)-\vec{y}(s'))^{k}} {|\vec{x}(s)-\vec{y}(s')|^{3}}

Although this is a topological invariant it ceases to be so for the self-linking \mathcal{GL}(K,\bar{K}) for when the curves coincide. However,  \mathcal{GL}(K, \bar{K}) can still always describe the properties of a ribbon with edges  K and \bar{K} that self-entangles or knots about its main helical axis K_{h}

There is a natural connection between basic knot theory and magnetostatics. Historically, the Gauss Linking Number arose from the Biot-Savart law of magnetostatics, found in 1820. If two currents \mathcal{J} and \bar{\mathcal{J}} flow in two loops K and \bar{K'} then the ‘action-at-a -distance’ formula for the mutual energy is

 E\sim -\mathcal{J}\bar{\mathcal{J}}\int_{K}\int_{\bar{K}}d\mathbf{r}.d\mathbf{r'} \frac{1}{|\mathbf{r}-\mathbf{r'}|}

which is essentially a knot energy functional. The power of Maxwell’s theory subsequently explained the result in terms of a vector potential \mathbf{A}(\mathbf{r}), with the magnetic field \mathbf{B}(t) essentially the ‘helicity’ \mathbf{B}(\mathbf{r})=\nabla\wedge \mathbf{A}(\mathbf{r}). Suppose  \mathbb{D}\subset\mathbb{R}^{3} and let (K,K')\subset\mathbb{D} be curves or knots within domain \mathbb{D} bounding 2-surfaces \Sigma and \Sigma' such that K=\partial\Sigma and k'=\partial \Sigma'. Let  \mathbf{B}=\mathbf{B}(x) be a magnetic field defined for all x\in\mathbb{D}. The field obeys the Maxwell equations so that

\nabla \wedge \mathbf{B}=\mu_{o}\mathcal{J},~~~\nabla . \mathbf{A}=0  where \mathcal{J} is the current density within the wire. Ampere’s Law follows from Stoke’s Theorem as

  \int_{K}\mathbf{B}.  dl  = {\displaystyle \int_{K}} (\nabla \wedge \mathbf{B}).ds = \mu_{o} {\displaystyle \int_{\Sigma}}   \mathcal{J}. ds=\mu_{o} I 

If the wire pierces the surface \Sigma m times then {\displaystyle \int_{K}} \mathbf{B}. d\mathbf{l}=\mu_{o} m I. If we know introduce the vector potential  \mathbf{A} then  B=\nabla\wedge \mathbf{A} and the Coulomb gauge is  \nabla \mathbf{A}=0 then  B=\nabla \wedge \nabla\wedge A = \nabla^{2}A=\mu_{o}\mathcal{J}. The solution for the vector potential is

  \mathbf{A}(\vec{r})=\frac{\mu_{o}}{4\pi}\int_{\mathbb{D}} d^{3}\mathbf{r'}\frac{J(\vec{r}')}{|\mathbf{r}-\vec{r}'|} = \frac{\mu_{o}}{4\pi}{\displaystyle \int_{K'}} \frac{I d\mathbf{l}'}{|\vec{r}-\vec{r}'|})

since  \mathcal{J}(\mathbf{r}')d^{\mathbf{r}'}=\mathcal{I} d\mathbf{l}' as the current only has support within the wire. The magnetic field is then

 \mathbf{B}(\mathbf{r})=\nabla\wedge\mathbf{A}(\mathbf{r})=\frac{I\mu_{o}}{4\pi}\nabla \wedge {\displaystyle \int_{K'} \frac {d\mathbf{l}}{|\mathbf{r}-\mathbf{r}'|}}  =-\frac{I\mu_{o}}{4\pi}{\displaystyle  \int_{K'}}\frac{(\mathbf{r}-\mathbf{r}^{\prime})\wedge d\mathbf{l}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}} 

Introducing a second curve  K and integrating  \mathbf{B} around  K gives

 \int_{K}\mathbf{B}(\mathbf{r}).d\mathbf{l}={\displaystyle \int_{K}}(\nabla\wedge\mathbf{A})(\mathbf{r}).d\mathbf{l}=-\frac{I\mu_{o}}{4\pi} {\displaystyle \int_{K}\int_{K'}}\frac{(\mathbf{r}-\mathbf{r}^{\prime})\wedge d\mathbf{l}^{\prime}.d\mathbf{l}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}=\mu_{o}\pi m 

It can then be seen that

 \mathbf{\mathcal{GL}}(K,K')=m=-\frac{1}{4\pi}{\displaystyle \int_{K}}{\displaystyle \int_{K'}}\frac{(\mathbf{r}-\mathbf{r}^{\prime})\wedge d\mathbf{l}^{\prime}.d\mathbf{l}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}=\frac{1}{4\pi} {\displaystyle \int_{K}\int_{K^{\prime}}} \epsilon_{ijk}dx^{i}dy^{j}\frac{(x-y)^{3}}{|x-y|^{3}} 

The rhs of (-) is the modern expression of the linking number but if \mathbf{r}=(x,y,z) and \mathbf{r}'=(x^{\prime},y^{\prime},z^{\prime}) then (-) it can be expressed as

 m=-\frac{1}{4\pi}{\displaystyle \int\int} \frac{(x-x')(dydz^{\prime}-dzdy^{\prime})+(y^{\prime}-y)(dzdx^{\prime}-dxdz^{\prime}+ (z-z^{\prime})(dxdy^{\prime}-dydx^{\prime})} {|(x^{\prime}-x)^{2}+(y^{\prime}-y)^{2}+(z^{\prime}-z)^{2}|^{3/2}} 

(Gauss C. F. 1. Zur mathematischen theorie der electrodynamischen Wirkungen (1833), Werke. K¨oniglichen Gesellschaft der Wissenschaften zu G¨ottingen, 5, 605, (1877).)
This is the original expression first presented by Gauss.

While this historical origin is very far removed from molecular biology, \mathcal{GL}(K,\bar{K}) also provides a natural description for a self-linking or self-entangling of a ribbon of DNA exhibiting tertiary structure or superhelicity. Here, the DNA ribbon self-intersects or knots about its own helical axis, leading to writhe and twist. A description of this is accomodated by the Calagareanu-White theorem which follows from the Gauss linking number, and first applied to a description of DNA topology by (refs.) Laboratory observed topological states of DNA are consistent with the C-W theorem

Let   R(K,\bar{K},K_{h},\beta,\pm \pm n)  be a (twisted) ribbon with edges (K,K') and width \beta, then if  \mathcal{LK} is identified with  \mathcal{GL}(K,K'), then the linking number is the sum of twist and writhe about the helical axis K_{h} such that in the limit as K\rightarrow K_{h} and \bar{K}\rightarrow K_{h} as \beta\rightarrow 0

 \mathbf{\mathcal{LK}}(K,K')=\mathbf{\mathcal{TW}}(K)+\mathbf{\mathcal{WR}}(K_{h},K_{h}) 

where the twist is defined as

\ \mathbf{\mathcal{TW}}(K_{h})=\frac{1}{2\pi}\int_{\mathfrak{S}} ds \frac{d\vec{x}(s)}{ds}.\left[\vec{N}(s)\wedge\frac {d\vec{N}(s)}{ds}\right]/[d\vec{x}(s)/ds]

and the writhe is

  \mathcal{WR}(K_{h})  = \frac{1}{4\pi}  {\displaystyle \int_{K_{h}}} ds     {\displaystyle \int_{K_{h}}}ds'    [\frac{d\vec{x}(s)}{ds}  \wedge    \frac{d\vec{y}(s')}{ds'}.\frac{\vec{x}(s)-\vec{y}(s')}{|\vec{x}(s)-\vec{y}(s')|}

where  \vec{N}(x(s))=\vec{N}(s) is the unit vector pointing to  y(s) and   y(s)=x(s)+\beta  \vec{N}(s). If   \vec{N}(s) is the normal vector then   \mathbf{T}_{w} is the total integrated torsion of the single curve  K_{h}

Let  \mathbb{M}^{2+1}=\mathbb{R}^{2}\times\mathbb{R} be a (2+1)-dimensional manifold and consider an Abelian \mathbb{G}=U(1) Chern-Simons gauge field with cpts. A_{\mu}(x)=(A_{o}(x),A_{1}(x),A_{2}(x)). The basic action is

   S_{CS}[A]=\frac{k}{2} \int_{\mathbb{M}} d^{3}x \epsilon^{\mu\nu\rho}A_{\mu}(x)\partial_{\nu} A_{\mu}(x) = \frac{k}{2}{\displaystyle \int_{\mathbb{M}}} A  \wedge  dA 

which is invariant under gauge transformations. This is also analogous to the helicity integrals of (-) except here the integral is over \mathbb{M}^{2+1} and not \mathbb{R}^{3}. One defines a product of Wilson loops

 W[\lbrace q_{i}\rbrace;A]]=\prod_{i=1}^{r} exp\left(iq_{i}\int_{\mathfrak{S}_{i}}A_{\mu}(x)dx^{\mu}\right)\equiv exp\left(iq_{i}\int_{\mathfrak{S}_{i}}A\right) 

for a set of curves \lbrace K_{i}\rbrace, the trajectories of particles with charges \lbrace Q_{i}\rbrace, with $i=1…n$. Gauge invariance requires that Q_{i} are integers (ref Yang.) The vacuum expectation values of the Wilson loops are given by a path integral

 \left\langle W[\lbrace q_{i}\rbrace;A]]\right\rangle =\frac{\left\langle\psi_{o}\|W[\lbrace q_{i}\rbrace;A]]\|\psi_{o}\right\rangle}{\left\langle\psi_{o}|\psi_{o}|\right\rangle} = \frac{{\displaystyle\int}\mathcal{D}A W[\lbrace q_{i}\rbrace ;A]]\exp(iS_{CS}[A])}{{\displaystyle\int} \mathcal{D}A exp(iS_{CS}[A])}

where |\psi_{o}\rangle is the vacuum state and \int\mathcal{D}A is the functional integral over the gauge field. The calculation can be done exactly (refs) and is

  \langle W[q_{i} ; A]  \rangle  =  exp\bigg(\frac{i}{2k}   \sum_{i,j} q_{i}  q_{j}  \mathcal{GL} (K_{i},K_{j}) \bigg)  

where   \mathcal{GL}(K_{i}, K_{j})  is the linking number. The non-abelian C-S theory can also be solved and connects with the Jones polynomial for knot theory.

To first order and using a framing such that    y_{\mu}(s)=x^{\mu}(s)+\beta n^{\mu}(s), then the self intersection is avoided since every curve   K_{i} is smeared out into thin ribbon of width   \beta with edges  K_{i} and  \bar{K}_{i} (ref Witten.) If  q_{i}=q_{j}=Q then

  \bigg( \langle W(K)  \bigg\rangle = exp \bigg(\frac{i q^{2}}{2K}\mathcal{GL}(K,\bar{K})\bigg) 

which describes a ribbon with edges  (K,\bar{K}). This is useful from the point of view of quantum computing where the ribbon edges can represent worldlines of anyons and braiding is achieved by twisting the ribbon. A twist in the ribbon gives a linking of the worldlines and is registered within the Wilson-loop expected values as a phase shift. For example, if a ribbon

  R(K, \bar{K},\beta)  is twisted by   +2\pi then   \mathcal{GL} (K, \bar{K}) \rightarrow \mathcal{GL}(K,\bar{K})+1 . Then

 \langle  W(K) \rangle \rightarrow   exp \bigg(\frac{i q^{2}}{2K}\bigg)\langle  W(K)  \rangle 

Again, while this is a powerful representation of the linking number in relation to fluctuating curves, ribbons and knot theory, it is a quantum theory. But for molecular biology–and indeed other applications such as knotted filament
lines/vortices in classical statistical turbulence and polymer fluctuations–it would be useful to have a similar mathematical representation but from a \emph{classical stochastic theory}. This can then represent classically fluctuating polymers and ribbon systems that have an inherent stochastic or random structure and behavior, but which fluctuate or writhe about an average geometry/topology. This is highly relevant to large molecules like DNA, RNA, proteins and other structures within cells or in-vitro solutions, which are essentially immersed in a ‘warm sea’ of smaller molecules which induce a Brownian-motion-type bombardment of the larger structure.
\subsection{Knot energy functionals}

The GLN arose originally from magnetostatics. Another form of knot functional can arise from electrostatics by considering charged knots or strings. Suppose we have a uniformly charged knot  K within a non-conducting fluid filling a domain  \mathbb{D}. Then a point   x\in K and a point  y can experience a Coulomb force. The configuration of the knot may the evolve to decrease its electrostatic energy to a minima. Charged polymers in a solution. Knot energy functionals for a single knot, in the mathematical sense, were considered in relation to an electrostatic anology with a modified Coulomb potential (refs.)

Let   K\subset\mathbb{D} \subset\mathbb{R}^{3} be a charged knot and let   \mathcal{D}_{K}(x,y) be the shorter arc length between  x\in K and y\in \bar{K}. The unregulated knot energy functional for K is then

  \overline{\mathbf{\mathcal{K}ef}}(K)=\delta_{ij}\int_{K}\int_{K}\frac{dx^{i}dy^{j}}{|x-y|^{2}}

which blows up at  x=y. To regulate   \overline{\mathcal{K}ef}(K,\bar{K}), the following are defined

  V^{j}(K:x,\epsilon)  = \int_{K,  y\in K, \mathcal{D}_{K}(x,y)\ge \epsilon}\frac{dy^{j}}{|x-y|^{2}}

  E(K,\epsilon)=\delta_{ij}\int_{K}V^{j}(K:x,\epsilon)dx^{i}=\delta_{ij} \int_{K} \int_{K,\mathcal{D}_{K}(x,y)\ge  \epsilon}  \frac{dx^{i}dy^{j}}{|x-y|^{2}} 

These have the limits (ref)

V^{j}(K:x,\epsilon)=lim_{\epsilon\rightarrow 0}\bigg(\int_{K, y\in K,\mathcal{D}_{K}(x,y)\ge \epsilon}\frac{dy^{j}}{|x-y|^{2}}-\frac{2}{\epsilon}\bigg) + 4

 E(K,\epsilon)=\delta_{ij} \int_{K}V^{j}(K:x,\epsilon)dx^{i}=\delta_{ij} \int_{K} \int_{K,\mathcal{D}_{K}(x,y)\ge \epsilon} \frac{dx^{i}dy^{j}}{|x-y|^{2}}-\frac{2}{\epsilon}\bigg) + 4

The term  V^{j}(K:x,\epsilon)=\lim_{\epsilon\rightarrow 0}\bigg(\int_{K, y\in K,\mathcal{D}_{K}(x,y)\ge \epsilon}\frac{dy^{j}}{|x-y|^{2}}-\frac{2}{\epsilon}\bigg) + 4 is the “voltage” at x when the subarc  \lbrace y\in K |D_{K}(x,y)\ge \epsilon \rbrace is charged. The regulated knot energy functional of O’Hara is then (ref)

\mathcal{K}ef(K)= E(K,\epsilon)=\delta_{ij}\bigg( \int_{K} \int_{K,\mathcal{D}_{K}(x,y)\ge \epsilon} \frac{dx^{i}dy^{j}}{|x-y|^{2}}-\frac{1}{(\mathcal{D}_{K}(x,y))^{2}}\bigg)dx^{i}dy^{j} 

For open knots, the  +4 term can be dropped and for  K a circle or straight line  \mathcal{K}ef(K)=0 .

Using the electrostatic analogy, knot energy type functionals can be defined for a pair of open or closed curves or knots K,\bar{K})\in\mathbb{D} Let (x,y)\in \mathbb{D} and let (q_{1},q_{2}) be a pair of charges at (x,y). The Coulomb law is then

\mathcal{F}_{q_{1}q_{2}} \sim \frac{1}{4\pi}\frac{q_{1}q_{2}}{|x-y|^{2}}

If the charges are smeared out into strings or curves (K,\bar{K}) then q_{1}=\int_{K}Q_{1}dx and q_{2}=\int_{\bar{K}}Q_{2}dy, where (\mathcal{Q}_{1},\mathcal{Q}_{2}) are charge densities. This gives

\mathcal{F}_{\mathcal{Q}_{1},\mathcal{Q}_{2}} \sim \frac{1}{4\pi}\int_{K}\int_{\bar{K}}\frac{\mathcal{Q}_{1}\mathcal{Q}_{2}}{|x-y|^{2}}dxdy

Setting \mathcal{Q}_{1}=\mathcal{Q}_{2}=1 we can define a non-regulated knot energy functional of the form

  \mathbf{\mathcal{K}ef}(K,\bar{K})=\frac{\delta_{ij}}{4\pi}\int_{K}\int_{\bar{K}}\frac{ dx^{i}dy^{j}}{|x-y|^{2}} 

Introducing a correlation length \zeta, the dimensionless form is

\mathbf{\mathcal{K}ef}(K,\bar{K})=\frac{\delta_{ij}\zeta^{2}}{4\pi}\int_{K}\int_{\bar{K}}\frac{ dx^{i}dy^{j}}{|x-y|^{2}}

We follow closely the formalism developed by Reuter and Meyer (ref), which gives a set of modified Einstein equations for a running Newton constant.  We will be concerned with how the ultraviolet behaviour involving black hole singularities can potentially be altered by a running Newton constant. To avoid confusion the following notation is utilised: \mathbf{G}_{\mu\nu} denotes the Einstein tensor, G is the usual Newton constant and the scalar field \mathcal{G}(x) is a position dependent Newton constant.

The action for pure Einstein gravity on a spacetime manifold \mathbb{M}^{3+1}=\mathbb{M}^{3}\times\mathbb{R} is

  S_{EH}[g,G]=\frac{1}{16\pi G}{\displaystyle \int_{\mathbb{M}^{3+1}}} d^{4}x\sqrt{-g}\mathbf{R}

For gravity interacting with fluids or matter fields, one can add another term to the action to account for a source term \mathbf{T}_{\mu\nu}, so that

 S_{tot}=\frac{1}{16\pi G} {\displaystyle \int_{\mathbb{M}^{3+1}}}d^{4}x\sqrt{-g}\mathbf{R}+S_{M}[g] 

and

\mathbf{T}_{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta S_{M}[g]}{\delta g_{\mu\nu}} .

Varying the total action \delta S_{tot}=0 gives the usual Einstein field equations

 \mathbf{G}_{\mu\nu}=\mathbf{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{R}= 8\pi G\mathbf{T}_{\mu\nu} 

Since G is constant, it remains part of a prefactor outside the integration over the Lagrangian and thus has a static rather than a dynamical role. If G is now elevated to the status of a scalar field \mathcal{G}(x) then \mathcal{G}(x) must now be brought under the action integral

 S_{EH}[g,G]=\frac{1}{16\pi}{\displaystyle \int_{\mathbb{M}}^{3+1}} d^{4}x \sqrt{-g}\frac{\mathbf{R}}{\mathcal{G}(x)}+S_{M}+S_{\Phi}

where S_{\Phi} is now a source for the scalar field \mathcal{G}(x). This will give rise to a scalar-tensor gravity theory that modifies the pure Einstein gravity theory.

The first such theory was the Brans-Dicke scalar-tensor theory. This theory originally sought to makes Mach’s principle compatible with a modified general relativity. The Brans-Dicke scalar field \phi(x) was introduced as the inverse of a position-dependent Newton constant so that \phi(x)=\mathcal{G}^{-1}(x). The action for the theory is

 S_{BD}=\frac{1}{16\pi}{\displaystyle \int} d^{4}x \sqrt{-g}(\phi\mathbf{R}-\omega\phi^{-1}\partial_{\mu}\phi \partial^{\mu}\phi)~+~S_{M} 

Varying the action gives the field equations

 \mathbf{G}_{\mu\nu}=8\pi(\phi^{-1}\mathbf{T}_{\mu\nu}+\mathbf{T}_{\mu\nu}^{BD})+\phi^{-1}(\nabla_{\mu}\nabla_{\nu}\phi-\square \phi)

and a scalar wave equation

 \square\phi(x) = \square G^{-1}(x)=\frac{8\pi}{3+2\omega}\mathbf{T}_{\mu}^{\mu}(x)

Here,  \mathbf{T}_{\mu\nu} is the source tensor arising from the matter while \mathbf{T}_{\mu\nu}^{BD} is the source tensor for the field \Phi (refs.)

 \mathbf{T}_{\mu\nu}^{BD}=\frac{\omega}{8\pi \phi}[\nabla_{\mu}\phi\nabla_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\nabla_{\alpha}\phi \nabla ^{\alpha}\phi]

The key idea in the BD theory is that the position-dependent Newton constant is governed by a simple equation of motion, the Klein-Gordon equation. In this more general theory of a position-dependent Newton constant(-),the BD theory can arise as a special limiting case.

Varying the modified action gives

  \frac{1}{16 \pi}{\displaystyle \int_{\mathbb{M}}} (-\mathbf{R}_{\mu\nu}+\frac{1}{2}   g_{\mu\nu}\mathbf{R}) \mathcal{G}(x))  +  (- \nabla_{\mu} \nabla_{\nu} + g_{\mu\nu}  \square) \mathcal{G}^{-1}(x)) \delta  g^{\mu \nu}  

The field equations that arise from the action (-) will have the general form

   \mathbf{R}_{\mu\nu} - \frac{1}{2}g_{\mu\nu} \mathbf{R})\mathcal{G}(x))= \mathcal{G}(x)\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square)   \mathcal{G}^{-1}(x))

or

 \mathbf{G}_{\mu\nu}=  \mathbf{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{R} = 8\pi \mathcal{G}(x)\mathbf{T}_{\mu\nu}+ \mathbf{\Psi}_{\mu\nu}+ \mathbf{\Phi}_{\mu\nu} 

or

 \mathbf{G}_{\mu\nu}=\mathbf{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{R}=8\pi  \mathcal{G}(x)(\mathbf{T}_{\mu\nu}+ \widehat{\mathbf{\Psi}}_{\mu\nu}+\widehat{\mathbf{\Phi}}_{\mu\nu} )

where \mathbf{\Phi}_{\mu\nu} is the source term for the scalar field \mathcal{G}(x) and \mathbf{\Psi}_{\mu\nu} is the tensor that arises from the position dependence of \mathcal{G}(x) within the variation of the action. The Bianchi identities and the conservation of energy and momentum demands that

  \nabla^{\mu}\mathbf{G}_{\mu\nu}=8\pi \nabla^{\mu} (\mathcal{G}(x)\mathbf{T}_{\mu\nu})+\nabla^{\mu}\widehat{\mathbf{\Psi}}_{\mu\nu}+\nabla^{\mu}\widehat{\mathbf{\Phi}}_{\mu\nu}=0

In the absence of matter

 \mathbf{G}{\mu\nu}=\Psi_{\mu\nu}+\mathbf{\Phi}_{\mu\nu}

  \nabla^{\mu}\mathbf{G}_{\mu\nu}= \nabla^{\mu}\widehat{\mathbf{\Psi}}_{\mu\nu}+\nabla^{\mu}\widehat{\mathbf{\Phi}}_{\mu\nu}=0

Explicitly, the tensor \mathbf{\Psi}_{\mu\nu}(x) is

 \widehat{\mathbf{\Psi}}_{\mu\nu}(x)=\mathcal{G}(x)(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square)G^{-1}(x) 

or

  \widehat{\mathbf{\Psi}}_{\mu\nu}=\mathcal{G}(x)^{-2}[2\nabla_{\mu}\mathcal{G}(x)\nabla_{\nu}\mathcal{G}(x) -\mathcal{G}(x)\nabla_{\mu}\nabla_{\nu}\mathcal{G}(x)-2g_{\mu\nu} \nabla_{\mu}\mathcal{G}\nabla^{\mu}\mathcal{G}(x)-\mathcal{G}(x)\square\mathcal{G}(x)]

The trace and covariant divergence are

 \widehat{\mathbf{\Psi}}_{\mu}^{\mu}=\frac{3}{\mathcal{G}^{3}(x)}[\mathcal{G}(x)\square \mathcal{G}(x)-2\nabla_{\mu}\nabla^{\mu}\mathcal{G}]

 \nabla^{\mu}\widehat{\mathbf{\Psi}}_{\mu\nu}=\frac{1}{\mathcal{G}(x)}\nabla^{\mu} \mathcal{G}(x)(\mathbf{\Psi}_{\mu\nu}-\mathbf{R}_{\mu\nu}] 

The consistency or constraint condition also gives

  \frac{1}{\mathcal{G}(x)}  \nabla^{\mu}\mathcal{G}(x)(\widehat{\Psi}_{\mu\nu}-R_{\mu\nu})+ \nabla^{\mu}\widehat{\Phi}_{\mu\nu}+8\pi (\nabla_{\mu}\mathcal{G}(x)\mathbf{T}^{\mu}_{\nu}=0

However, the following alternative expression for the consistency condition also provides the link to the original BD theory:

 \frac{3}{2}\frac{\nabla_{\nu}\mathcal{G}(x)}{\mathcal{G}(x)^{2}}[\mathcal{G}(x)\square \mathcal{G}(x)-2(\nabla_{\mu}\mathcal{G}\nabla^{\mu}\mathcal{G}(x)))]+\nabla^{\mu}  \widehat{\Phi}_{\mu\nu} -\frac{\nabla^{\mu}\mathcal{G}(x)}{\mathcal{G}(x)}\widehat{\Phi}_{\mu\nu}+4\pi \mathbf{T}\nabla_{\nu}\mathcal{G}(x)=0 

where \mathbf{\Phi}_{\mu\nu}=\mathbf{\Phi}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\Phi_{\alpha}^{\alpha} and \mathbf{T}=\mathbf{T}_{\mu}^{\mu}. When \mathbf{T}=0, then the energy-momentum tensor is traceless giving

\ \frac{3}{2}\frac{\nabla_{\nu}\mathcal{G}}{\mathcal{G}(x)^{2}}[\mathcal{G}(x)\square \mathcal{G}(x)-2(\nabla_{\mu}\mathcal{G}(x)\nabla^{\mu}\mathcal{G}(x)))]+\nabla^{\mu}\widehat{\Phi}_{\mu\nu} -\frac{\nabla^{\mu}\mathcal{G}9x)}{\mathcal{G}}\widehat{\Phi}_{\mu\nu}=0

The tensor  \widehat{\mathbf{\Phi}}_{\mu\nu} which solves this equation has the form

 \mathbf{\Phi}_{\mu\nu}^{BD}= \bigg(-\frac{3}{2}\bigg)\frac{1}{8\pi \mathcal{G}(x)^{3}}[\nabla_{\mu}\mathcal{G}(x)\nabla_{\nu}\mathcal{G}(x)-\frac{1}{2}g_{\mu\nu}\nabla_{\mu} \mathcal{G}(x)\nabla^{\mu}G(x)] 

For the choice of the Brans-Dicke parameter \omega=\frac{3}{2} and \mathbf{\Phi}=\mathcal{G}^{-1}, this is then equivalent to the original Brans-Dicke source tensor (-) The action S_{\Phi} this tensor is

 S _{\Phi}=\frac{3}{32\pi}{\displaystyle \int} d^{4}x \sqrt{g}\mathcal{G}^{-2}\nabla_{\mu}\mathcal{G}\nabla^{\nu}\mathcal{G}

The modified Einstein equations are now formulated on a spherically symmetric metric. This is the standard static isotropic form

ds^{2}=-X(r)dt^{2}+Y(r)dr^{2}+r^{2}(d\theta^{2}+ sin^{2}(\theta)d\varphi^{2} 

For the Killing vectors K^{\mu} we require that K^{\mu}\partial_{\mu}G(x)=0 so that G=G(r) and so G can vary with the radial coordinate r only. The radial derivatives will be denoted D_{r}=d/dr and D_{rr}=d^{2}/dr^{2}. It will be instructive to initially formulate the full generic modified field equations on this metric, and then solve for the specific case when there is no matter and no cosmological constant term. This will illucidate how the solution of the modified vacuum equations compares to the Standard Schwarzschild vacuum solution and the interior Schwarzschild solution. The full modified Einstein equations are

 \mathbf{G}_{\mu\nu}=\mathbf{R}_{\mu\nu}-\frac{1}{2} g_{\mu\nu}\mathbf{R} = -g_{\mu\nu}\Lambda    + 8 \pi \mathcal{G}(x) \mathbf{T}_{\mu\nu}(x) + \widehat{\Psi}_{\mu\nu}(x) + \widehat{\Phi}_{\mu\nu}(x) = \mathbf{T}tot_{\mu}^{\nu}  

where \widehat{\mathbf{\Psi}}_{\mu\nu}=\mathbf{\Psi}_{\mu\nu}/8\pi \mathcal{G}(x) and \widehat{\mathbf{\Phi}}_{\mu\nu}=\mathbf{\Phi}_{\mu\nu}/8\pi \mathcal{G}(x). The source tensor for fluid system in hydrostatic equilibrium has the usual form

 \mathbf{T}_{\mu\nu}=p g_{\mu\nu}+(p+\rho)\xi_{\mu}\xi^{nu}

with an equation of state p=p(\rho) so that  \mathbf{T}_{\mu\nu}=diag[-\rho(r),p(r),p(r),p(r)]. The conservation equation arising from the Bianchi identities is \nabla^{\mu}\mathbf{T}_{\mu\nu}=0 and leads to the condition

p^{\prime}(r)+\frac{X^{\prime}}{2X}(\rho+p)=0 

which is essentially the Tolman-Oppenheimer-Volkoff equation for relativistic stellar structure.

Next,the source tensor arising from the spatial dependence of Newton’s constant is expressed in terms of density and pressures, in the same form as for the matter fluid source so that

\\mathbf{\Psi}_{\mu\nu}=\frac{\widehat{\mathbf{\Psi}}_{\mu\nu}}{8\pi \mathcal{G}}=diag[-\mathcal{D}(r),\mathcal{P}(r),\mathcal{P}_{a}(r),\mathcal{P}_{a}(r)]

Here, \mathcal{D}=-\mathbf{\Psi}_{t}^{t} is an energy density induced by the position-dependence of Newton’s constant, and \mathcal{P}(r)=\mathbf{\Psi}_{r}^{r} is the associated radial pressure and
 \bar{\mathcal{P}}(r)=\mathbf{\Psi}_{\theta}^{\theta}=\mathbf{\Psi}_{\varphi}^{\varphi}, are the angular or tangential pressures.

For the metric (-), the density and pressures are induced within the vacuum due to the running of the Newton constant are(refs)

  \mathcal{D}(r)=-\mathbf{\Psi}_{t}^{t}=\frac{1}{8\pi \mathcal{G}(r)}\bigg( \frac{2}{rY} \mathcal{G}^{-1}(r)\mathcal{G}^{\prime}(r)

   + \frac{Y^{\prime}}{2Y^{2}}  \mathcal{G}^{-1}(r)  

    \mathcal{G}^{\prime}(r) +      \frac{2}{Y} \mathcal{G}(x)^{-2}(r) \mathcal{G}(x)^{\prime}(r)  \mathcal{G}^{\prime}(r) +\frac{1}{Y}\mathcal{G}^{-1}(r)\mathcal{G}^{\prime\prime}(r)\mathcal{G}(r)\bigg)

  \mathcal{P}(r)=\mathbf{\Psi}_{r}^{r}=\frac{1}{8\pi \mathcal{G}(r)}\bigg( \frac{X^{\prime}}{XY}\mathcal{G}^{-1}(r)\mathcal{G}^{\prime}(r)+ \frac{2}{Y}\mathcal{G}^{-1}(r)\mathcal{G}^{\prime}(r) \bigg)

  \bar{\mathcal{P}}(r)=\mathbf{\Psi}_{\theta}^{\theta}=\mathbf{\Psi}_{\varphi}^{\varphi} =\frac{1}{8\pi \mathcal{G}(r)}\bigg(\frac{1}{rY}\mathcal{G}^{-1}(r)\mathcal{G}^{\prime}(r) -\frac{Y^{\prime}}{2Y^{2}}G^{-1}(r)\mathcal{G}^{\prime}(r)\nonumber\\[0.1in]+ \frac{X^{\prime}}{2XY}\mathcal{G}^{-1}(r)G^{\prime}(r)-\frac{2}{Y}\mathcal{G}^{-2}(r)\mathcal{G}^{\prime}(r) \mathcal{G}^{\prime}(r) -\frac{1}{Y}\mathcal{G}^{-1}(r)\mathcal{G}^{\prime\prime}(r)\bigg)

The \Phi_{\mu\nu} source tensor will be left arbitrary for the moment, but it is sufficient that it has same diagonalised form

\mathbf{\Phi}_{\mu}^{\nu}=diag[-\mathscr{D}(r),\mathscr{P}(r),\bar{\mathcal{P}}(r),\bar{\mathcal{P}}(r)]

The final contribution can come from a standard cosmological constant term \Lambda or vacuum ‘dark energy’ with no spatial dependence such that \mathbf{\Theta}_{\mu}^{\nu}=diag[-\rho_{v},p_{v},p_{v},p_{v}] with \rho_{v}=\Lambda/8\pi \mathcal{G} and p_{v}=\Lambda/8 \pi \mathcal{G}. The total source tensor is then

 \mathbf{T}_{\mu}^{\nu}=diag[-D(r),P(r),\bar{P}(r),\bar{P}(r)] 

so that

D(r)=\rho(r)+\rho_{v}+\mathcal{D}(r)+\mathscr{D}(r)

P(r)=p(r)+p_{v}+\mathcal{P}(r)+\mathscr{P}(r)

 \bar{P}(r)=p(r)+p_{v}+\bar{\mathcal{P}}(r)+\bar{\mathscr{P}}(r)

The modified Einstein equations can then be written as

 \mathbf{G}_{\mu\nu}=8\pi G\mathbf{T}_{\mu\nu}^{(tot)}

For the metric (-) the nonvanishing components of the Einstein tensor are
\mathbf{G}_{t}^{t},\mathbf{G}_{r}^{r},\mathbf{G}_{\theta}^{\theta},\mathbf{G}_{\varphi}^{\varphi} so that the full field equations are

 \mathbf{G}_{t}^{t}=-8\pi \mathcal{G} D(r)=-8\pi \mathcal{G}[\rho(r)+\rho_{v}+\mathcal{D}(r)+\mathscr{D}(r)]

\mathbf{G}_{r}^{r}=~8\pi\mathcal{G} P(r)=~8\pi G[p(r)+p_{v}+\mathcal{P}(r)+\mathscr{P}(r)]

 \mathbf{G}_{\theta}^{\theta}=\mathbf{G}_{\varphi}^{\varphi}=8\pi P_{a}(r)=8\pi \mathcal{G}[p(r)+p_{v}+\mathcal{P}_{a}(r)+\bar{\mathcal{P}}(r)]

In terms of the metric functions X(r) and Y(r)

 \mathbf{G}_{t}^{t}=\frac{1}{r^{2}Y}-\frac{1}{r^{2}}-\frac{Y^{\prime}}{r Y^{2}}=-8\pi G(\rho(r)+\rho_{v}+\mathcal{D}(r)+\mathscr{D}(r)) 

 \mathbf{G}_{r}^{r}=\frac{1}{r^{2}Y}-\frac{1}{r^{2}}-\frac{A^{\prime}}{r Y^{2}}=~8\pi \mathcal{G}(r)(p(r)+p_{v}+\mathcal{P}(r)+\mathscr{P}(r)) 

 \mathbf{G}_{\theta}^{theta}=\mathbf{G}_{\varphi}^{\varphi}= -\frac{Y^{\prime}}{2rY^{2}}+\frac{X^{\prime}}{2r XY}-\frac{X^{\prime} Y^{\prime}}{4B^{2}X}-\frac{(X^{\prime})^{2}}{4YX^{2}}+ \frac{\partial_{rr}X^{ \prime\prime}}{2XY}=8\pi \mathcal{G} [p(r)+p_{v}+\mathcal{P}_{a}+\bar{\mathscr{P}}]

The conservation of energy and momentum requires that \nabla^{\mu}\mathbf{T}_{\mu\nu}=0 which gives

 \frac{dP(r)}{dr}+ \frac{X^{\prime}}{2X}(D+P))+\frac{2}{r}(\mathcal{P}-\bar{\mathcal{P}})=0

Given the formal properties of GRVFs on \mathbf{R}^{3} and the existence of stochastic integrals, Wilson lines can be constructed.

Let  K \in\mathbf{D}\subset\mathbf{R}^{3} be an open string/knot or smooth curve in the simply connected region \mathbf{D} with Euclidean coordinates x^{i}(s)\in K \subset \mathbf{D}, where s \in [0,1] and with ends x^{i}(0)=\alpha and x^{i}(1)=\beta. Now let X_{i}(x) with \mathcal:\mathbf{R}^{3}\rightarrow\mathbf{R}^{3} be a smooth non-random vector field defined for all x\in\mathbf{D}. Then the classical Wilson-line holonomy is

 W[K]=\mathcal{P} exp\bigg(\mu{\displaystyle \int_{0}^{s}}d\bar{s}{X}_{i}(x(\bar{s})) \frac{dx^{i}}{d\bar{s}}d\bar{s}\bigg)= \vec{\mathcal{P}} exp\bigg ({\displaystyle \int_{K}}{X}_{i}(x)dx^{i} \bigg)

and for closed strings/curves K

W[K]=\vec{\mathcal{P}} exp\bigg({\displaystyle \oint_{\gamma}}X_{i}(x)dx^{i}\bigg)

Proof:
Let W(s) with W:[0,1]\rightarrow \mathbb{R}^{+} be a smooth differentiable function on K. Given the line integral for L[K]={\displaystyle \int_{K}}{X}_{i}(x)dx^{i}, define a second integral

\Psi(K,W(s))={\displaystyle \int_{K}}{X}_{i}(x)W[s]dx^{i} ={\displaystyle \int_{0}^{s}}d\bar{s}{X}_{i}(x(\bar{s}))\frac{dx^{i}(\bar{s})}{d\bar{s}} {W}(\bar{s})

then this is equivalent to the differential equation

 {\displaystyle \frac{d\Psi(K,W[s])}{ds}=X_{i}(x(s))\frac{dx^{i}(s)}{ds}W[s]}

Choosing \Psi(K, W[s])=W[s] then

 {\displaystyle \frac{dW(s)]}{ds}=X_{i}(x(s))\frac{dx^{i}(s)}{ds}W(s)}

The integral solution is

  W(s)=W(0)+{\displaystyle \int_{0}^{s}} ds_{1} X_{i}(x(s_{1}))\frac{dx^{i}(s)}{ds_{1}}W[s_{1}]

and describes transport of the function W[s] along the curve K. Performing one iteration  W[s]=W[0]

+\int_{0}^{s}ds_{1}{X}_{i}(x(s_{1}))\frac{dx^{i}(s)}{ds}\bigg( W[0]+{\displaystyle \int_{0}^{s_{1}}}ds_{2}{X}_{i}(x(s_{2}))\frac{dx^{i}(s_{2})}{ds_{2}} W[s_{2}]\bigg)

Successive iterations can be continued indefinitely and the solution W(s) is a convergent infinite series; that is, an exponential so that the general solution is

 W(s)={P} exp\bigg (\int_{0}^{s}d\bar{s}{X}_{i}(x(\bar{s})){dx^{i}(\bar{s})} {d\bar{s}}\bigg)

 = \sum_{M=0}^{\infty}\frac{1}{m!}{\displaystyle \int_{0}^{s}}ds_{1}...{\displaystyle \int_{0}^{s_{M-1}}}ds_{M} \mathcal{P}[ X_{i_{1}}(x(s_{1}))\times...\times X_{i_{n}}(s_{n})\frac{dx^{i_{1}}(s_{1})}{ds_{1}}\times...\times \frac{dx^{i_{M}}(s_{M})}{ds_{M}}

where P is a ‘normal-ordering’ operator. This is then a classical Wilson-line operator in \mathbf{R}^{3}.

 W[K]=\vec{\mathcal{P}} exp\bigg(\mu{\displaystyle \int_{0}^{s}}d\bar{s}X_{i}(x(\bar{s}))\frac{dx^{i}(\bar{s})}{d\bar{s}}\bigg)= \vec{\mathcal{P}} exp\bigg({\displaystyle \int_{K}}{X}_{i}(x)dx^{i}\bigg) 

and for closed curves \mathfrak{K}

  W[K] =\vec{\mathcal{P}}  exp \bigg({\displaystyle \oint_{K}} X_{i}(x)dx^{i}\bigg)  

the stochastic Wilson line follows from coupling a random field to  X_{i}(x)

Let  x\in\mathbb{D}\subset\mathbf{R}^{3}. Let  K be a curve or knot in \mathbb{D} so that  x\in K\subset\mathbb{D}. The vector field X_{i}(x) exists for all x\in\mathbb{D}. Let \widehat{\mathcal{F}}_{i}(x;\omega) be a random vector field defied with respect to a probability space (\Omega,\mathfrak{F},\mathbf{P}) such that \exists map \mathfrak{M}:(x,\omega)\rightarrow\mathcal{X}_{i}(x,\omega) for all x\in\mathbb{D} and \omega\in\Omega. Then the total field at any x\in\mathbb{D} is

 \widehat{X}_{i}(x)=X_{i}(x)+\widehat{\mathcal{X}}_{i}(x;\omega)

The total line integral along K is

 \widehat{L}[K]=L[K]+\widehat{\mathcal{L}}[K]= {\displaystyle \int_{K}}X_{i}(x)dx^{i}+{\displaystyle\int_{K}}\widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}

The corresponding stochastic Wilson line is then

  \widehat{W}[K]=exp\bigg({\displaystyle \int_{K}}X_{i}(x)dx^{i}+{\displaystyle \int_{K}}\widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}\bigg)

 =exp\bigg({\displaystyle\int_{K}}X_{i}(x)dx^{i}\bigg) exp{\displaystyle \int_{K}}\widehat{\mathcal{X}}_{i}(x,\omega)dx^{i} \bigg)

The expectation is

\ \mathfrak{W}[K]=\mathbf{E}(\widehat{W}[K])=exp \bigg({\displaystyle \int_{K}}X_{i}(x)dx^{i}+{\displaystyle \int_{K}} \widehat{\mathcal{X}}_{i}(x, \omega)dx^{i} \bigg)

= exp({\displaystyle \int_{K}}X_{i}(x)dx^{i})\mathbf{E} exp{\displaystyle \int_{K}} \widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}

= exp{\displaystyle \int_{K}}V_{i}(x)dx^{i}){\displaystyle \int_{\Omega}\int_{\Omega}} exp{\displaystyle \int_{K}} \widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}) d\mathbf{P}(\omega)d\mathbf{P}(\xi)

= {\displaystyle\int_{\Omega}\int_{\Omega}} exp{\displaystyle \int_{K}}\widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}) d\mathbf{P}(\omega)d\mathbf{P}(\xi)

=exp{\displaystyle \int_{K}}\widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}) d\mathbf{P}(\omega)d\mathbf{P}(\xi))

If the background field is homogeneous and constant then X_{i}(x)=X_{i} for  x\in\mathbb{D} so that \int_{K}X_{i}dx^{i}=X_{i}l=C . Hence $latex W[K]=W=const.$

The expectation of the stochastic Wilson line is

 \mathbf{\mathfrak{W}}[K]=\mathbf{E}(\widehat{\mathfrak{W}}[K])={\displaystyle \int_{\Omega}} exp{\displaystyle \int_{K}} \widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}d\mathbf{P}(\omega) 

For a closed curve or knot K

 \mathbf{\mathfrak{W}}[K]=\mathbf{E}(\widehat{\mathfrak{W}}[K])={\displaystyle \int_{\Omega}} exp{\displaystyle \oint_{K}} \widehat{\mathcal{X}}_{i}(x,\omega)dx^{i}d\mathbf{P}(\omega)

We can now present the main cluster integral expansion technique that enables the stochastic expectation of the Wilson loop to be evaluated, provided the random vector field \mathcal{X}_{i}(x;\omega) is Gaussian.Such expansions have also been considered in (-

Random Gaussian fields applied to fluctuating polymers and biopolymers, eg protein chains and DNA. First, non-conducting knotted charged strings or polymers in a non-conducting viscous fluid leads to the concept of knot energy functionals. The electrosatic energy will minimise without intersection due to Coulomb forces. One can then define a knot-energy funcional. Let (K,\bar{K}) be two such knots in a domain \mathbb{D} such that (K,\bar{K})\subset\mathbb{D}\subset\mathbb{R}^{3} . Let x^{i}\subset K and  y^{j}\subset \bar{K}. A regulated knot energy functional \mathbf{\mathcal{K}ef}(K,\bar{K}) is defined as

 \mathbf{\mathcal{K}ef}(K,\bar{K}) = - \delta^{ij}lim_{\epsilon \rightarrow 0} \bigg({\displaystyle \int_{K}}{\displaystyle \int_{\bar{K}}} {\displaystyle \frac{dx^{i} dy^{j} } {|x-y|^{2}}}-\frac{2}{\epsilon}\bigg)  

Similary, two knots carrying unit currents leads to the Bio-Savart law

 \mathbf{\mathcal{K}f}(K,\bar{K}) = \delta_{ij}\bigg({\displaystyle \int_{K}}{\displaystyle \int_{\bar{K}}} {\displaystyle \frac{dx^{i} dy^{j} } {|x-y|} }\bigg)  

\bf{Theorem}

Let \widehat{\mathcal{Q}}_{i}(x,\omega) be a random Gaussian vector field defined for all x\in\mathbb{D}\subset\mathbb{R}^{3} and with respect to a triplet \Omega,\mathcal{F},\mathbf{P}). Then if:

\mathbf{E}(\widehat{\mathcal{Q}}_{i}(x,\omega)={\displaystyle\int}_{\Omega}d\mathbf{P}(\omega)\widehat{\mathcal{Q}}_{i}(x,\omega) =0

\mathbf{E}(\widehat{\mathcal{Q}}_{i}(x,\omega)\otimes \widehat{\mathcal{Q}}_{j}(y,\omega)

={\displaystyle\int}_{\Omega}d\mathbf{P}{\displaystyle\int}_{\Omega}d\mathbf{P}(\omega)\widehat{\mathcal{Q}}_{i}(x,\omega)(\omega)\widehat{\mathcal{Q}}_{j}(y,\eta))=  {\displaystyle \frac{\zeta^{2}}{|x-y|^{2}}}

The expectation of the stochastic Wilson line of \mathcal{Q}_{i}(x,\omega) gives

\mathbf{E}(\widehat{W}[K])={\displaystyle \int_{\Omega}}d\mathbf{P}(\omega)exp\bigg({\displaystyle\int_{K}}  \widehat{\mathcal{Q}}_{i}(x,\omega)dx^{i}\bigg)  exp \bigg({\displaystyle \int_{K}\int_{\bar{K}}\delta_{ij}\zeta^{2} \frac{dx^{i}dy^{j}}{|x-y|^{2}}}\bigg)

\mathbf{E}(\widehat{W}[K])={\displaystyle \int_{\Omega}}d\mathbf{P}(\omega)exp  \bigg(i{\displaystyle\int_{K}}  \widehat{\mathcal{Q}}_{i}(x,\omega)dx^{i}\bigg)  exp \bigg(-{\displaystyle \int_{K}\int_{\bar{K}}\delta_{ij}\zeta^{2} \frac{dx^{i}dy^{j}}{|x-y|^{2}}}\bigg)

if:

\mathbf{E}(\widehat{\mathcal{Q}}_{i}(x,\omega)={\displaystyle \int}_{\Omega}d\mathbf{P}(\omega)\widehat{\mathcal{Q}}_{i}(x,\omega) =0

\mathbf{E}(\widehat{\mathcal{Q}}_{i}(x,\omega)\otimes \widehat{\mathcal{Q}}_{j}(y,\omega)

={\displaystyle\int}_{\Omega}d\mathbf{P}{\displaystyle\int}_{\Omega}d\mathbf{P}(\omega)\widehat{\mathcal{Q}}_{i}(x,\omega)(\omega)\widehat{\mathcal{Q}}_{j}(y,\eta))=  {\displaystyle \frac{\zeta}{|x-y|}}

The expectation of the stochastic Wilson line of  \widehat{\mathcal{Q}}_{i}(x,\omega) gives

\mathbf{E}(\widehat{W}[K])={\displaystyle \int_{\Omega}}d\mathbf{P}(\omega)exp\bigg({\displaystyle\int_{K}}\widehat{\mathcal{Q}}_{i}(x,\omega)dx^{i}\bigg)  exp \bigg({\displaystyle \int_{K}\int_{\bar{K}}\delta_{ij}\zeta \frac{dx^{i}dy^{j}}{|x-y|}}\bigg)

\mathbf{E}(\widehat{W}[K])={\displaystyle \int_{\Omega}}d\mathbf{P}(\omega)exp\bigg(i{\displaystyle\int_{K}}\widehat{\mathcal{Q}}_{i}(x,\omega)dx^{i}\bigg)  exp \bigg(-{\displaystyle \int_{K}\int_{\bar{K}}\delta_{ij}\zeta \frac{dx^{i}dy^{j}}{|x-y|}}\bigg)

 

Having established existence of the derivative of a SRVF, stochastic integrals can be defined as a mean-square Riemann integration. Let \widehat{\mathcal{U}}_{i}(x) be a RVF on a domain \mathbb{D}\subset\mathbb{R}^{3} and for all x^{i},y^{j}\in\mathbb{D}, let \psi(x,y) be a deterministic continuous and bounded function such that \psi:\mathbb{R}^{3}\times\mathbb{R}^{3}\rightarrow\mathbb{R}. The function can also be complex such that  \psi:\mathbb{R}^{3}\otimes\mathbb{R}^{3}\rightarrow\mathbb{C}. The stochastic volume integral is the mean-square Riemann integral

 \widehat{\mathcal{R}}_{i}(y)={\displaystyle\int_{\mathbb{D}}}\psi(y,x)\otimes\widehat{\mathcal{U}}_{i}(x,\omega)d^{3}x 

The integral exists if the limit of the Riemann sum exists
 \widehat{\mathcal{R}}_{i}^{(m)}(x)=\sum_{\zeta=1}^{m}\psi(y,x_{\zeta})\otimes \widehat{\mathcal{U}}_{i_{\zeta}}(x_{\xi},\omega)\Delta_{vol}(x_{\zeta}) 

where \Delta(x_{\zeta}) is a volume element centered on x_{\xi} and \sum_{\zeta}\Delta(x_{xi})=vol(\mathbb{D}). The integral then exists if \widehat{\mathcal{R}}_{i}(y)=\lim_{m\rightarrow \infty} \mathcal{I}^{(m)}_{i}(x)
Since \mathbf{E}(\widehat{\mathcal{U}}(x))=0 then \mathbf{E}\langle\widehat{\mathcal{Z}}(y)\rangle=0. When \psi(x,y)=1 then

 \widehat{\mathcal{R}}_{i}(y)={\displaystyle \int_{\mathbb{D}}}\widehat{\mathcal{U}}_{i}(x,\omega)d^{3}x

This definition leads to the following corollary (See Gilrod ref…).
A SRVF \widehat{\mathcal{U}}_{i}(x) is mean-square Riemann integrable iff

 \mathbf{E}\langle\widehat{\mathcal{R}}_{i}(y)\otimes\widehat{\mathcal{R}}_{j}(y')\rangle =

 {\displaystyle  \int_{\mathbb{D}}} {\displaystyle \int_{\mathbb{D}}} \psi(x, y) \psi(x', y')\mathbf{C}ov_{ij}(x,x')d^{3} x d^{3} x' <   lim_{m\rightarrow\infty} \sum_{\xi=1}^{m} \sum_{\eta=1}^{m} \psi(y, x_{\xi})\psi(y',x_{\eta}')  \mathbf{E}  (\widehat{\mathcal{U}}_{i_{\xi}}(x_{\xi})\otimes \widehat{\mathcal{U}}_{j_{\eta}}(x'))\Delta_{vol}(x_{\xi}) \Delta_{vol}(x_{\eta}))

For a Gaussian RVF this is equivalent to

 \mathbf{E}(\widehat{\mathcal{R}}_{i}(y)\otimes\widehat{\mathcal{R}}_{j}(y'))={\displaystyle \int_{\mathbb{D}}}{\displaystyle \int_{\mathbb{D}}}\psi(x,y)\psi(x',y')\mathbf{E}(\widehat{\mathcal{U}}(x)\otimes\widehat{\mathcal{U}}(x')) d^{3}xd^{3}x'<\infty

The definition can be extended to include line integrals. Let \widehat{\mathcal{U}}_{i}(x) be a GRVF on a domain \mathbb{D}\subset\mathbb{R}^{3} and for all x^{i}\in\mathfrak{S}\subset\mathbb{D} and y^{j}\in\bar{\mathfrak{S}}\subset\mathbb{D}. Let \psi(x,y) be a deterministic continuous and bounded function such that  \psi:\mathbb{R}^{3}\times\mathbb{R}^{3}\rightarrow\mathbb{R}.

\widehat{\mathcal{R}}(y)={\displaystyle \int_{\mathfrak{S}}}\psi(y,x)\otimes\widehat{\mathcal{U}}_{i}(x)dx^{i} 

The integral exists if the limit of the Riemann sum exists
\widehat{\mathcal{I}}^{(m)}(x)=\sum_{\xi=1}^{m}\psi(y,x_{\xi})\otimes\widehat{\mathbb{U}}_{i_{\xi}}(x_{\xi})\Delta x_{\xi} ,

where \mathcal{R}x_{\xi} is a line element centered on x_{\xi} and \sum_{\xi}\Delta x_{xi} is the length of the line. The stochastic line integral then exists if

 \widehat{\mathcal{R}}(y)=\lim_{m\rightarrow \infty}\widehat{\mathcal{R}}^{(m)}(x). If \psi(x,y)=1 then

\widehat{\mathcal{L}}(y)={displaystyle\int_{\mathfrak{S}}}\widehat{\mathcal{U}}_{i}(x)dx^{i}

\widehat{\mathcal{L}}(y)={\displaystyle \oint_{\mathfrak{S}}}\widehat{\mathcal{U}}_{i}(x)dx^{i}

For a GRVF, the linking \gamma\bigcap\gamma' of two curves or knots is described by the correlations

\mathbf{E}\langle\widehat{\mathcal{L}}(x)\otimes\widehat{\mathcal{L}}(y)\rangle={\displaystyle{\int_{\mathfrak{S}}}{\displaystyle\int_{\bar{\mathfrak{S}}}}\mathbf{E}\langle\widehat{\mathcal{U}}_{i}(x) \otimes\widehat{\mathcal{U}}_{i}(y)\rangle dx^{i}dy^{j}

 \mathbf{E}\langle\widehat{\mathcal{L}}(x)\otimes\widehat{\mathcal{L}}(y)\rangle =\oint_{\mathfrak{C}} \oint_{\bar{\mathfrak{\gamma}}}\mathbf{E}\langle \widehat{\mathcal{U}}_{i}(x)\otimes \widehat{\mathcal{U}}_{i}(y)\rangle dx^{i}dy^{j}  

The n-point correlation for a link comprised of n curves or knots

 \gamma_{1}\bigcap\gamma_{2}\bigcap...\gamma_{n-1}\bigcap\gamma_{n}

is

  \langle\widehat{\mathcal{L}}(x_{1})\otimes...\otimes\widehat{\mathcal{L}}(x_{n})\rangle = \int_{\mathfrak{S}_{1}}...\int_{\mathfrak{S}_{1}}d^{3} x_{1}...d^{3} x_{n} \mathbf{E} \langle\widehat{\mathcal{U}}_{i_{1}}(x_{1})\otimes...\otimes \widehat{\mathcal{U}}_{i_{1}}(x_{1})\rangle  

which can be expressed in a path integral form as

  \langle\widehat{\mathcal{L}}(x_{1})\otimes...\otimes\widehat{\mathcal{L}}(x_{n})\rangle= \int D_{n}x\mathbf{E}\langle\widehat{\mathcal{U}}_{i_{1}}(x_{1})\otimes...\otimes\widehat{\mathcal{U}}_{i_{1}}(x_{1})\rangle 

We can now examine the spectral properties of Gaussian random vector field.
Given the stochastic integral \widehat{\mathcal{Z}}(y)=\int_{\mathbb{R}^{3}}d^{3}x f(x,y)\otimes\widehat{\mathcal{U}}_{i}(x) of (-), then setting y=k with k  \in\mathbb{R}^{3} and  f(x,k)=\frac{1}{2\pi}^{3/2}\exp[-ik_{i}x^{i}] gives a stochastic Fourier integral.

\widehat{\mathcal{U}}_{i}(k)=\frac{1}{2\pi^{3/2}}\int d^{3}k\widehat{\mathcal{U}}_{i}(x)\exp[ik_{i}x^{i}] 

If there exists such a random vector field \widehat{\mathcal{U}}_{i}(k) then the random vector field
\widehat{\mathcal{U}}_{i}(k) is harmonizable. The integral exists in the mss iff

 \mathbf{E}(\widehat{\mathcal{U}}_{i}(x)\otimes\widehat{\mathcal{U}}_{j}(y))

=\frac{1}{2\pi^{3}} \int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}d^{3}x d^{3} y    exp[-i (x_{i} k^{i})]   exp[-ik_{j}y^{j}] \mathcal{C}ov_{ij}(x,y)<\infty

and  \mathcal{C} ov_{ij}(x, y) = \mathbf{E} (\mathcal{U}_{i}(x) \otimes \mathcal{U}_{j}(y) )  for a Gaussian RVF. Equation (-) defines the cross spectral density function, provided the integral exists.

Given the transform pair

 \widetilde{\mathcal{U}}_{i}(x)=(2\pi)^{-3/2}\int_{\mathbb{R}^{3}}d^{3}x \widetilde{\mathcal{U}}_{i}({x}) exp[-i k  x ]

  \widetilde{\mathcal{U}}_{i}(k)= (2\pi)^{-3} \int_{\mathbb{R}^{3}} d^{3}k \widetilde{\mathcal{U}}_{i}({k}) exp[+i k x]

there is a technical difficulty in that the integrals do not actually converge when integrated over all space. The problem is that the random function \widetilde{\mathcal{B}}({x},t) extends over all space and does not decay to zero as  |{x}|\rightarrow \infty, as is required for the existence of a classical Fourier transform. The resolution of this difficulty is essentially to define the transform of  \widetilde{\mathcal{U}}_{i}({x}) over a finite volume or domain \mathbb{D} of 3-space so that \mathbb{D}\subset \mathbb{R}^{3} and consider the limit as the volume grows large and approaches infinity.

A finite-range transform \widetilde{\mathcal{U}}_{i}^{[l]}({k}) and its
complex conjugate \widehat{\mathcal{U}}_{i}^{[l]*}({k}) can be defined as

  \widehat{\mathcal{U}}_{i}^{[l]}({k}) = (2\pi)^{-3}\int _{\mathbb{D}_{l}}d^{3}x~\widetilde{\mathcal{U}}_{i}({x},t) exp[-i k x]

\widehat{\mathcal{U}}_{i}^{[l]*}({k})= (2\pi)^{-3}\int _{\mathfrak{D}_{l}} d^{3}y \widehat{\mathcal{U}}_{i}({y})exp[+i{k}.{y}] 

with (x,y)\in\mathbb{D} and where \mathbb{D}_{l}\subset \mathbb{R}^{3} is a cubic volume defined by  x_{1}\in [-l,l], x_{2}\in [-l,l], x_{3}\in [-l,l] and 2 l defines the sides of the cube. Taking the product of the transform and its conjugate then averaging gives the spectral moment

 \mathbf{E}(\widehat{\mathcal{U}}^{[l]}({k})\circ\widehat{\mathcal{U}}^{l*}({k})) =(2\pi)^{-6}\int_{\mathfrak{D}_{l}}d^{3}x\int_{\mathbb{D}_{l}}d^{3}y G_{ij}({x}-{y},t)\exp[-i{k}.({x}-{y})]

We then consider the behavior of (36) in the limit as \mathbb{D}_{l}\rightarrow\infty or l\rightarrow \infty. The correlation G_{ij}(x,y) is a function of  \arrowvert{x}-{y}\arrowvert owing to its homogeneity, and there is a decorrelation for \arrowvert{x}-{y}\arrowvert\rightarrow \infty. The length scale over which  G_{ij}(x-y) has significant non-zero values is \mid\!\! {x}-{y}\!\!\mid =O(\ell_{*}) with \ell_{*} being the correlation length. Let  r=|x-y| be a vector separation of two points. Replacing y as an integration variable in (36) gives

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}^{[l]}({k})\otimes\widetilde{\mathcal{U}}_{j}^{[l*]}({k}))=(2\pi)^{-6} \int_{\mathbb{D}_{l}} d^{3}x \int_{\mathbb{D}_{l}'} d^{3}r G_{ij}({r},t) exp[-i{k}.{r}]

where \mathbb{D}_{l}' depends on on x and is defined by r_{1}\in [x_{1}-l,x_{2}+l], r_{2}\in [x_{2}-l,x_{2}+l], r_{3}\in[x_{3}-l,x_{3}+l]. The cube \mathbb{D}_{l} has sides  2L and is centred on  r= x. The integral is only significantly different from zero when $| r |=O(l)$ and this region lies well within \mathbb{D}_{l}', provided that l is many correlation lengths from the boundary of \mathbb{D}_{l}.Then as l\rightarrow \infty for most  x\in \mathbb{D}_{l}, the integration over r can be carried to infinity.

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}({k})\times\widetilde{\mathcal{U}}_{j}({k}))) \sim (2\pi)^{-3} \int_{\mathbb{D}_{l}}d^{3}x \psi_{ij}({k},t)=(l/\pi)^{3}\psi_{ij}(k) 

The approximation of the integral over r by \psi_{ij} applies everywhere except when x is of the order of a correlation length from the boundary of \mathbb{D}_{l}.

We are most interested in the case when (38) has different values of k. The spectral correlation is where the integration variable has been changed from y to r=x-y . Again. provided that |x| is many correlation lengths from the boundary of \mathbb{D}_{l} the integral over r can be
extended to infinity so that

 \mathbf{E}((\widetilde{\mathcal{U}}_{i}^{[l]}({k},t)\circ\widetilde{\mathcal{U}}_{j}^{[l*]}(k)]] )\sim\psi_{ij}(k',t)\Xi_{\ell}(k-k')

where
\Xi_{l}(k-k')= \frac{sin[(k_{1}-k_{1}')l]sin[(k_{2}-k_{2}')l]sin[(k_{3}-k_{3}')l]}{\pi^{3} (k_{1}-k_{1}')(k_{2}-k_{2}')(k_{3}-k_{3}')}

The function  \Xi_{\ell}({k}-{k}') has a maximum peak at  k=k' and becomes increasingly sharply peaked about k=k' as  l\rightarrow\infty with the integral \int d^{3}k\Xi_{l}(k)=1 always satisfied. It follows that \Xi_{l}(k-k') has all the properties of a delta function so that \Xi_{l}(k-k')\rightarrow \delta^{3}(k-k') as l\rightarrow\infty. This then leads
to the important result

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}({k})\circ \widetilde{\mathcal{U}}_{j}({k}')) =\psi_{ij}({k},t)\delta^{3}(k-k') 

We have replaced  \psi_{ij}(\vec{k}') with \psi_{ij}({k},t) since the delta function is effectively zero when x\ne {k}'. Since \widetilde{\mathcal{S}}_{i}(x,t) is real the complex conjugate of (-) gives \widetilde{\mathcal{U}}_{i}^{*}(k')=\widetilde{\mathcal{U}}_{j}(-k') so that the spectral correlation (in the infinite l limit) becomes

 \mathbf{E}({\mathcal{U}}_{i}(k)\circ\widetilde{\mathcal{U}}_{j}(k)) =\psi_{ij}(k,t)\delta^{3}(k+k') 

Spectral analysis of turbulence for example, involves Fourier analysis of a stochastic Navier-Stokes flow on \mathbb{R}^{3} so that \widetilde{\mathcal{U}}_{i}(x)=U_{i}(x)+\widetilde{\mathcal{V}}_{i}(x) The essential idea is to decompose the turbulent fluctuations about the mean flow into sinusoidal components and study the distribution of turbulent energy amongst the different wave vectors representing the different scales of turbulence.

The spectral function \psi_{ij}(k-k',t) is given by the Fourier transform
\mathcal{F}_{T}:\mathbb{R}^{3}\rightarrow \mathbb{R}^{3} of the stress tensor S_{ij}(x,y;t).

 \psi_{ij}(k,t)=(2\pi)^{-3}\int_{\mathbb{R}^{3}}d^{3}x S_{ij}(x,y,t)\exp[-k.y-x]=(2\pi)^{-3} *

 \int d^{3}x\mathbf{E}(\arrowvert\widetilde{\mathcal{U}}_{i}(x,t)\circ\widetilde{\mathcal{U}}_{j}(y,t)\arrowvert exp[-ik.r]

and depends on the wavenumber vector k and t as the turbulence evolves. The inverse transform is  R_{ij}(\vec{r},t)=\int d^{3}k \psi_{ij}({k},t)exp(i {k}.{r}), which is an integral over all 3-dimensional wavenumber space. Setting \vec{r}=0 and i=j in (-) gives the stress tensor of (-)

 \frac{1}{2}\mathbf{E}(\widetilde{\mathcal{U}}_{i}(x)\circ\widetilde{\mathcal{U}}_{j}({x}))=\frac{1}{2}S_{ij}(0,t)=\frac{1}{2}\int_{\mathbb{R}^{3}} d^{3}k \psi_{ij}(k,t)

The ‘turbulent kinetic energy’ can therefore be expressed as an integral over k. In turbulence theory for example, where \widehat{\mathcal{U}}_{i}(x) is a random Navier-Stokes flow, the function \psi_{ij}({k},t)/2 can be interpreted as the distribution of turbulent energy over the various wavenumbers. The stresses and spectral functions have the symmetries S_{ij}(r)=S_{ij}(-r) and \psi_{ij}(k,t)=\psi_{ji}(-k,t). Also \psi_{ij}^{*}(k,t)=\psi_{ji}(-k,t) leading to \psi_{ij}^{*}(k,t)=\psi_{ji}(k,t), which implies that \psi_{ij} is a Hermitian matrix. Each of the diagonal terms is therefore real as should be the case for a quantity interpreted as the distribution of energy with wavenumber. The magnitude of the wavenumber vector k=\mid {k}\mid has dimensions of inverse length. One can interpret the length scale \ell=(k^{-1}) as the spatial scale represented by wavenumber k . For example if \ell is an integral scale of turbulence the wavevectors \ell^{-1} represent the scales of turbulence that contain most of the kinetic energy. The spectral function then has its maximum values for \mid {k}\mid= O(l^{-1}) and decays to zero as |k| \rightarrow \infty. The stress tensor can then be represented as the inverse Fourier transform

 S_{ij}(x,y,t)=\mathbf{E}(\widetilde{\mathcal{U}}_{i}(x)\times\widetilde{\mathcal{U}}_{j}(y)) =\int_{\mathbb{R}^{3}} d^{3}k\int_{\mathbb{R}^{3}} d^{3}k' \psi_{ij}(k)\delta^{3}(k-\beta')\exp[i k'.(x-y)]

where \beta=  - k and the integration range is taken over all \mathbb{R} space since we take the spectral function to be such that the Fourier integrals now converge. Equation (44) becomes

 S_{ij}(x,y,t)=\mathbf{E}(\widetilde{\mathcal{U}}_{i}(x,t) \times\widetilde{\mathcal{U}}_{j}(y,t)) = (2\pi)^{-3}\int_{\mathbb{R}} d^{3}k \psi_{ij}({k},t)e^{i{\beta}.(x-y)} 

Now exp[i{\beta}.(x-y)]=exp[-ik.(x-y)]=\exp[ik.(y-x)] so that

 S_{ij}(x,y,t)=\mathbf{E})\lbrace\widetilde{\mathcal{U}}_{i}(x,t) \otimes\widetilde{\mathcal{U}}_{j}({y},t) \rbrace=(2\pi)^{-3}\int_{\mathbb{R}^{3}} d^{3}k~\psi_{ij}(k,t)exp[i {k}.(y-x)] 

Finally, for homogeneity and isotropy it is required that

 \frac{d}{d\xi_{j}}\mathbf{E}(\widetilde{\mathcal{U}}_{i}(x,t)\otimes\widetilde{\mathcal{U}}_{j}({x}+{\xi},t)) = (2\pi)^{-3}\int_{\mathbb{R}^{3}} d^{3}k ~ik_{j}\psi_{ij}(k,t)exp[ik_{j}\xi_{j}]=0 

so that the spectral function must satisfy k_{j}\psi{ij}({k})=0. The most general form is then

  \psi_{ij}(k)=\left(\delta_{ij}- \frac{k_{i}k_{j}}{k^{2}}\right)f(k) 

The 2-point correlations then depend on the ansatzes for the spectral function.
Let \psi_{ij}(k)=(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}} be the spectral function corresponding to the 2-point function

\mathbf{E}(\widetilde{\mathcal{U}}_{i}(x,t)\times\widetilde{\mathcal{U}}_{j}(y,t))= \frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}d^{3}k \psi_{ij}(k)\exp[ik^{i}(x-y)_{i}]

The ansatz for the spectral function and relationship between the wave vectors k_{i} then determines the form of the correlation for the GRVFS \widehat{\mathcal{U}}_{i}(x) such that

[1] If k_{i}k_{j}=0 then \psi_{ij}(k)=\delta_{ij} and

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}(x)\times\widetilde{\mathcal{U}}_{j}(y))=\frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}d^{3}k \delta_{ij}(k)\exp[ik^{i}(x-y)_{i}]=\delta_{ij}\delta^{(3)}(x-y) 

This correlation then describes a white-in-space noise.

[2]  If k_{i}k_{j}=\epsilon_{ijk}k^{k} and \psi_{ij}(k)=\frac{\epsilon_{ijk}k^{k}}{k^{2}} then

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}({x})\circ\widetilde{\mathcal{U}}_{j}(y))= \frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}d^{3}k \epsilon_{ijk}\frac{k^{k}}{k^{2}} exp[ik^{i}(x-y)_{i}]=\frac{1}{4\pi}\epsilon_{ijk}|x-y|^{k}\|x-y\|^{-3}

which is of the same mathematical form as the propagator in a Chern-Simons theory(ref). The following lemma on the possible 2-point correlations is also stated.

Let \widehat{\mathcal{F}}_{i}(x) be a Gaussian random field defined for all x\in\mathbb{D} with \mathbf{E}[\widehat{\mathcal{F}}_{i}(x)]=0. Let x^{i}\in\mathfrak{S}\subset\mathbb{D} and y^{j}\in\bar{\mathfrak{S}}\subset\mathbb{D}. If the 2-point function is a ‘white-in-space noise’ of the form

 \mathbf{E}[\widehat{\mathcal{F}}_{i}(x)\times\widehat{\mathcal{F}}_{i}(x)]= \delta_{ij}\delta^{3}(x-y)

then
[1] There exists a random Gaussian field \widehat{\mathcal{B}}_{i}(x) for all x\in\mathbb{D} with the 2-point correlation
 \mathbf{E}[\widehat{\mathcal{B}}_{i}(x)\times\widehat{\mathcal{B}}_{i}(y)]= \delta_{ij}|x-y|^{-1}
[2] There exists a random Gaussian field \widehat{\mathcal{G}}_{i}(x) for all x\in\mathbb{D} with the 2-point correlation

 \mathbf{E}[\widehat{\mathcal{G}}_{i}(x)\times\widehat{\mathcal{G}}_{i}(y)]= \frac{1}{4\pi}\epsilon_{ijk}(x-y)^{k}|x-y|^{-3}

Let  \mathbf{E}[\widehat{\mathcal{R}}_{i}(x)\otimes\widehat{\mathcal{R}}_{j}(y)]=K_{ij}(x-y) then

 K_{ij}(x-y)=-\frac{\delta_{ij}}{4\pi}\nabla_{(x)}^{2}(|x-y|^{-1})= \delta_{ij}\delta^{(3)}(x-y)

 K_{ij}(x-y) = -\frac{\delta_{ij}}{4\pi}\nabla_{(x)}^{2}  G(x-y)=\delta_{ij}\delta^{(3)}(x-y) 

which follows from the Poisson equation

  -\frac{1}{4\pi} \nabla_{(x)}^{2}|x-y|^{-1}=\delta^{(3)}(x-y) 

Hence

\mathbf{E}[\widehat{\mathcal{B}}_{i}(x)\otimes\widehat{\mathcal{B}}_{j}(x)]\equiv G_{ij}(x-y)=\delta_{ij}|x-y|^{-1}

or G(x-y)=|x-y|^{-1}. Next, given G(x-y) there is a Greens function C_{ij}(x-y) such that

 C_{ij}(x-y)=\frac{1}{4\pi}\epsilon_{ijk}\nabla^{k}G(x-y)= \frac{1}{4\pi}\epsilon_{ijk}\nabla^{k}|x-y|^{-1}

so that

 C(x-y)=\mathbf{E}[\widehat{\mathcal{G}}_{i}(x)\otimes\widehat{\mathcal{G}}_{j}(x)] =\frac{1}{4\pi}\epsilon_{ijk}(x-y)^{k}|x-y|^{-3}

Let (x,y)\in\mathbb{D}\subset\mathbb{R}^{3} and \widehat{\mathcal{U}}_{i}(x) and \widehat{\mathcal{U}}_{i}(x) are GRVFs on \mathbb{D}. As in (-), the spectral function choice \psi_{ij}(k)=\delta_{ij} gives a 2-point correlation for white noise.

\mathbf{E}(\widetilde{\mathcal{U}}_{i}({x})\times\widetilde{\mathcal{U}}_{j}({y})) =\frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}d^{3}k \delta_{ij}(k)\exp[ik^{i}(x-y)_{i}]=\delta_{ij}\delta^{(3)}(x-y) 

Let \mathbb{B}\subset\mathbb{D} be a ball of radius \xi and vol(\mathbb{B})=\tfrac{4}{3}\pi\xi^{3}. If (x,y)\in\mathbb{B} the delta-function is ‘smeared-out’ over \mathbb{B} into a very narrow and highly peaked rectangular or top-hat function of width \xi. Then

 \delta_{ij}\delta(x-y)\rightarrow \xi^{-2}\Pi[\xi^{-1}(|x-y|-\tfrac{1}{2}\xi)] 

where

 \xi^{-2}\Pi(\xi^{-1}|x-y|-\tfrac{1}{2}\xi)=1~for |x-y|\le \xi

 \xi^{-2}\Pi(\xi^{-1}|x-y|-\tfrac{1}{2}\xi)0 for x>\xi

so that the correlation vanishes outside \mathbb{B}. This is equivalent to choosing

\psi_{ij}(k)=\tfrac{1}{2\pi}^{1/2}\delta_{ij}sinc[k/2] = \tfrac{1}{2\pi}^{1/2}k^{-1}\sin[k/2]

for the spectral function so that the Fourier transform is

 \mathscr{F}[sinc[k/2]]=\mathbf{E}(\widetilde{\mathcal{U}}_{i}({x})\times\widetilde{\mathcal{U}}_{j}({y}))

 =\frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}d^{3}k\delta_{ij}~sinc[k/2] \exp[ik^{i}(x-y)_{i}]=\delta_{ij}\Pi[x-y]

and then making the rescaling |x-y|\rightarrow [|x-y|-\tfrac{1}{2}\xi]\xi^{-1}.

 \mathbf{E}(\widetilde{\mathcal{U}}_{i}({x})\times\widetilde{\mathcal{U}}_{j}({y}))= \frac{\alpha}{\xi^{2}}\delta_{ij}\Pi[\xi^{-1}(|x-y|-\tfrac{1}{2}\xi)]

Previous theorems established that there is zero probability of a  diffusion blowup such that \mathbf{P}[\psi(t)=\infty]=0 for all t\in[t_{\epsilon},\infty)=\mathbf{Y}\bigcup\mathbf{Z}\subset\mathbf{R}^{+} where

 \mathbf{Y}\bigcup\mathbf{Z} = [t_{\epsilon},t_{*})\bigcup (t_{*},\infty)

and t=t_{*} is the blowup time for the unperturbed system. From the perspective of consistency, the same conclusion should also follow from a Fokker-Planck equation, or Kolmogorov’s second equation, which is essentially an equation for conservation of probability. For a pure diffusion with non drift, the FP equation will reduce to a diffusion equation for the underlying probability density function.

Let \psi(t) be a generic diffusion satisfying athe nonlinear SDE

d\widehat{\psi}(T)=U[\psi(t)]dt +V[\psi(t)]\circ d\widehat{\mathcal{B}}(t) 

with  U:\mathbf{R}\times \mathbf{R}\rightarrow \mathbf{R} and  V:\mathbf{R}\times \mathbf{R}\rightarrow \mathbf{R}

some initial data  \psi(0),\widehat{\mathcal{B}}(0)=0. with \psi(t)\in [\psi(0),\infty). A pure drift-free diffusion is then

 d\widehat{\psi}(T)=V[\psi(t)]\circ d\widehat{\mathcal{B}}(t)=\lim_{U\rightarrow 0}U[\psi(t)]dt +V[\psi(t)]\circ d\widehat{\mathcal{B}}(t)

The corresponding nonlinear Fokker-Planck equation for a probability density   \mathscr{P}[\psi(T),t] is (refs)

 \partial_{t}\mathscr{P} [\psi(T),t]=-\partial_{\psi} V[\psi(t)]+\frac{1}{2}V[\psi(t)] \partial_{\psi} V[\psi(t)]]\mathscr{P}[\psi(t),t]+\frac{1}{2}\partial_{\psi}\partial_{\psi}[V[\psi(t)]^{2} \mathscr{P}[\psi(t),t]  

Taking the Ito interpretation

 \widehat{\psi}(t+\delta t)-\widehat{\psi}(T)=U[\psi(t)]\delta t +V[\psi(t)]{\displaystyle \int_{t}^{t+\delta t}}d\mathcal{B}(t')

The Fokker-Planck equation now becomes (ref)

\partial_{t}\mathscr{P}[\psi(t),t]=\partial_{\psi}U[\psi(t),t]\mathscr{P}[\psi(t),t]+\frac{1}{2}\partial_{\psi}\partial_{\psi}[V[\psi(t)]^{2}\mathscr{P}[\psi(t),t] 

If U[\psi(t)=0 then this is a nonlinear diffusion equation for \mathscr{P}[\psi(t),t]

 \partial_{t}\mathscr{P}[\psi(t),t]=+\frac{1}{2} \partial_{\psi}\partial_{\psi}[V[\psi(t)]^{2}\mathscr{P}[\psi(t),t] 

The general solution of a nonlinear Fokker-Planck equation cannot usually be found but a stationary solution is possible in the infinite-time relaxation limit such that

 \partial_{t}\mathscr{P}[\psi(t),t]=0 so that

 G\mathscr{P}[\psi(t)]=\partial_{\psi} V[\psi(t),t]\mathscr{P}[\psi(t]+\frac{1}{2} \partial_{\psi}\partial_{\psi}[V[\psi(t)]^{-2}\mathscr{P}[D(t]= 

where G is the generator of the diffusion as in (-). The stationary solution at t=\infty for the probability density is

 \mathscr{P}[\psi(t),t=\infty]= \bigg|\frac{C}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}}^{\psi(t)}\frac{U[\bar{\psi}(t)]}{(\Psi[\bar{\psi}(t)])^{2}}d\bar{\psi}(t)\bigg)  

where $ latex C$ is a constant. The normalised solution is

 {\displaystyle \int_{\psi(0)}^{\infty}}\mathscr{P}[\psi(t)]d\psi(t)={\displaystyle \int_{\psi_{0}}^{\infty}} \bigg|\frac{C}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]}{(V[\bar{\psi}(t)])^{2}}d\bar{\psi}(t)\bigg)=1

so that

  \mathscr{P}[\psi(t),t=\infty]= \frac{\bigg|\frac{1}{(V[\psi(t),t=\infty]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{iU[\bar{\psi}(t)]} {(V[\bar{\psi}(t)])^{2}}d \bar{\psi}(t)\bigg)}{{\displaystyle \int_{\psi_{0}}^{\infty}} \bigg|\frac{1}{(V[\psi(t)]^{2}} exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]}{(V[\bar{\psi}(t)])^{2}} d \bar{\psi}(t)\bigg)}  

Integrating over the probability density gives

\mathbf{P}[\widehat{\psi}(t),t=\infty]={\displaystyle \int_{\psi(0)}^{\psi(t)}}\mathscr{P}[\psi(t),t=\infty]dD(t)= {\displaystyle \int_{\psi(0)}^{\psi(t)}}\frac{\bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]} {(V[\bar{\psi}(t)])^{2}}d\bar{\psi}(t)\bigg)}{{\displaystyle \int_{\psi_{0}}^{\infty}} \bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]}{(V[\bar{\psi}(t)])^{2}} d\bar{\psi}(t)\bigg)} d\psi(t)

The probability that \psi(t)<\infty at t=\infty is then

\mathbf{P}[\widehat{\psi}(t)<\infty,t=\infty]={\displaystyle \int_{\psi(0)}^{\infty}}\mathscr{P}[\psi(t),t=\infty]d\psi(t)= {\displaystyle \int_{\psi(0)}^{\infty}}\frac{\bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]} {(V[\bar{\psi}(t)])^{2}}d\bar{\psi}(t)\bigg)}{{\displaystyle \int_{\psi_{0}}^{\infty}} \bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg| exp\bigg(2{\displaystyle \int_{\psi_{0}}^{\psi(t)}}\frac{U[\bar{\psi}(t)]}{(V[\bar{\psi}(t)])^{2}} d\bar{\psi}(t)\bigg)} d\psi(t)=1

For a pure diffusion with V[\psi(t)]=0 the solution reduces to

\mathbf{P}[\widehat{\psi}(t)<\infty]=\lim_{\psi(t)\rightarrow\infty}{\displaystyle \int_{\psi(0)}^{\psi(t)}}\mathscr{P}[\psi(t),t=\infty]=  \lim_{\psi(t)\rightarrow\infty}\frac{\displaystyle{\int_{\psi(0)}^{\psi(t)}}\bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg|}{{\displaystyle \int_{\psi_{0}}^{\infty}} \bigg|\frac{1}{(V[\psi(t)]^{2}}\bigg|}=1 T

Technical details of the FP equation are found in (refs.)

Given the SDE d\widehat{\psi}(t)=V[\psi(t)]\otimes d\widehat{\mathcal{B}}(t) defined for all t\in\mathbf{Y}\bigcup\mathbf{Z} with initial matter density  \psi_{\epsilon} the Fokker-Planck diffusion equation now becomes (ref)

 \partial_{t}\mathscr{P}[\psi(t),t]=+\frac{1}{2} \partial_{\partial}\partial_{\psi}((\psi(t))^{4}(\psi(t)-1))\mathscr{P}[\psi(t),t] 

For a stationary solution in the limit as t\rightarrow\infty

 G\mathbf{\mathscr{P}}[\psi(t),t=\infty] = \frac{1}{2}\partial_{\psi}\partial_{\psi}((\psi(t))^{4}(\psi(t)-1))\mathbf{\mathscr{P}}[\psi(t),t=\infty]=0

The probability that the diffusion \psi(t) is finite at t=\infty is then

  \mathbf{P}[\psi(T)<\infty,t=\infty]=\lim_{\psi(t)\rightarrow\infty} {\displaystyle \int_{\psi_{\epsilon}}^{\psi(t)}}\mathscr{P}[\psi(t),t=\infty]d\psi(t)

= \lim_{\psi(t)\rightarrow\infty}\frac{{\displaystyle \int_{\psi_{\epsilon}}^{\psi(t)}}\bigg|\frac{1}{(\psi(t))^{4}(\psi(t)-1)}\bigg|} {\int_{\psi_{\epsilon }}^{\infty} \bigg|\frac{1}{(\psi(t))^{4}(\psi(t)-1)}\bigg|}=1

This result can also be computed explicitly. Given the SDE d\widehat{\psi}(t)=V[\psi(t)]\circ d\widehat{\mathcal{B}}(t) defined for all
t\in\mathbf{Y}\bigcup\mathbf{Z} with initial matter density \psi_{\epsilon} the Fokker-Planck diffusion equation is as before(ref)

 \partial_{t}\mathscr{P}[[\psi(t),t] =+\frac{1}{2} \partial_{\psi}\partial_{\psi}\psi^{4}(t)(\psi(t)-1))\mathbf{\mathscr{P}}[\psi(t),t]

For a stationary solution in the limit as t\rightarrow\infty

 G\mathscr{P}[\psi(t),t=\infty]=\frac{1}{2}\partial_{\psi}\partial_{\psi}\psi^{4}(t))(\psi(t)-1))\mathscr{P}[\psi(t),t=\infty]=0

The stationary solution for the Fokker-Planck probability density is

  \mathscr{P}[\psi(t),t=\infty]=  \bigg(\frac{6(\psi_{\epsilon})^{3}} {6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}-2\psi_{\epsilon}-3\psi_{\epsilon}^{4}}\bigg) \frac{1}{(\psi(t))^{4}(\psi(t)-1)}

The probability that the diffusion \widehat{\psi}(t) is finite at \infty is now given by (-) as

 \mathbf{P}[\psi(t)<\infty,t=\infty]={\displaystyle\int_{\psi_{\epsilon}}^{\infty}}\mathscr{P}[\psi(t),t=\infty]   \bigg(\frac{6\psi_{\epsilon}^{3}}{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- \psi_{\epsilon}-3\psi_{\epsilon}^{4}}\bigg) {\displaystyle \int_{\psi_{\epsilon}}^{\infty}}\frac{d\psi(t)}{(\psi(t))^{4}(\psi(t)-1)}=1 

Performing the integral i (-)

  \mathbf{P}[\psi(t)<\infty,t=\infty]={\displaystyle \int_{\psi_{\epsilon}}^{\infty}}\mathscr{P}[\psi(t),t=\infty] =  \bigg(\frac{6\psi_{\epsilon}^{3}}{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- 2\psi_{\epsilon}-3\psi_{\epsilon}^{4}} \bigg)  

 \times\bigg(\frac{1}{3(\psi(t))^{3}}+\frac{1}{2(\psi(t))^{2}}+\frac{1}{\psi(t)}+ ln (1-(\psi(t))^{-1})\bigg)_{\psi(t)=\infty}  -\bigg(\frac{6\psi_{\epsilon}^{3}}{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- 2\psi_{\epsilon}-3\psi_{\epsilon}^{4}}\bigg) \bigg(\frac{1}{3(D_{\epsilon})^{3}} -\frac{1}{2(\psi_{\epsilon})^{2}}-\frac{1}{\psi_{\epsilon}}-ln (1-\psi_{\epsilon}^{-1})\bigg)  

The first term vanishes at \psi(t)=\infty since ln(1)=0 so that (-) reduces to

  \mathbf{P}[\psi(t)<\infty,t=\infty]={\displaystyle \int_{\psi_{\epsilon}}^{\infty}}\mathscr{P}[\psi(t),t=\infty]

  = -\bigg(\frac{6\psi_{\epsilon}^{3}}{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- 2\psi_{\epsilon}-3\psi_{\epsilon}^{4}}\bigg)

\times\bigg(\frac{1}{3(\psi_{\epsilon})^{3}} -\frac{1}{2(\psi_{\epsilon})^{2}}-\frac{1}{\psi_{\epsilon}}-ln (1-\psi_{\epsilon}^{-1})\bigg) =-\bigg(\frac{6\psi_{\epsilon}^{3}}{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- 2\psi_{\epsilon}-3\psi_{\epsilon}^{4}}\bigg)\times-\bigg(\frac{6\psi_{\epsilon}^{3} ln(\psi_{\epsilon}^{-1}-1)-6\psi_{\epsilon}^{3}- 2\psi_{\epsilon}-3\psi_{\epsilon}^{4}}{6\psi_{\epsilon}^{3}}\bigg)=1

Hence, the diffusion remains finite in the infinite-time relaxation limit.

General relativity breaks down at the space-time singularities associated with the Big Bang at the birth of the universe, and at the black-hole singularities which arise from the total gravitational collapse of massive stars which have exhausted their thermonuclear fuel. Global and topological techniques developed by Hawking and Penrose, demonstrated that singularities are totally unavoidable within pure general relativity provided specific conditions hold: namely that matter satisfies the strong energy condition (SEC) such that  \mathbf{T}_{\mu\nu}\xi^{\mu}\xi^{\nu}\ge 0 or \rho>0 and \sum_{i}p_{i}\ge 0, for density and pressures, and that the space-time manifold \mathbf{\mathcal{M}} must satisfy appropriate casuality conditions such as global hyperbolicity.

Gravitational collapse and singularity formation are unavoidable once a trapped surface \mathbf{\mathcal{T}} is formed around a sufficiently heavy collapsing star which is well above the Oppenheimer-Volkoff limit, and which has exhausted its fuel. The star can then no longer provide the thermal and radiative transport/pressure required to maintain the equilibrium condition required by the Tolman-Oppenheimer-Volkoff equation –or its Newtonian limit–and which kept it stable for billions of years. Most of the stars in the night sky are within the Oppenheimer-Volkoff limit and will ultimately collapse to super-dense white dwarfs or neutron stars following a red giant phase of nuclear helium burning. Since the (thermo) nuclear reaction rates R  for fusion of hydrogen (and heavier elements) in its core varies  as R \sim M^{3}, heavier stars burn through their fuel much more rapidly–a star 10 time heavier than the sun will burn through its fuel about a 1000 times faster. In reality, the heaviest stars will supernova but a superheavy central core can still completely collapse to a point.

It is remarkable that the Tolman-Oppenhimer-Volkoff equation follows directly from the energy conservation constraint imposed on the Einstein equations applied to a self-gravitating and isotropic spherically symmetric fluid or gas. The proper density, pressure and gravitational fields within a spherically symmetric star are described by the Einstein equations coupled to a perfect fluid-source tensor:

\mathbf{G}_{\mu\nu}=\mathbf{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=-8\pi G \mathbf{T}_{\mu\nu}

or

 \mathbf{R}_{\mu\nu}=-8\pi G \mathbf{T}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathbf{T}_{\gamma}{\gamma}

where \mathbf{T}_{\mu\nu}=pg_{\mu\nu}+(p+\rho)U_{\mu}U_{\nu} and U_{\mu} is the velocity 4-vector defined so that g^{\mu\nu}U_{\mu}U_{\nu}. Since the fluid is at rest U_{r}=U_{\theta}=U_{\phi}=0 and U_{t}=-(h(r))^{1/2}. It is a standard (albeit tedious) calculation to derive the Riemann tensor components and Einstein field equations for a spherical perfect-fluid star (refs) with an exterior metric

 ds^{2}=-f(r)dr^{2}+h(r)dr^{2}+r^{2}(sin^{2}\theta d\theta^{2}+d\varphi^{2})

so that

 \mathbf{G}_{tt}=-8\pi G \mathbf{T}_{tt}

 \mathbf {G}_{rr}=8\pi G \mathbf{T}_{rr}

 \mathbf{G}_{\theta\theta}=\mathbf{G}_{\varphi\varphi}= 8\pi  G \mathbf{T}_{\theta\theta}

The solution is the interior Schwarzchild metric:

ds^{2}=-\bigg(1-\frac{2GM(r)}{r}\bigg)dt^{2}+\bigg(1-\frac{2GM(r)}{r}\bigg)^{-1}dr^{2}+r^{2}d\Omega^{2}

where the interior mass up to radius r within the star is

M(r)=4\pi {\displaystyle \int_{0}^{r}}\rho(r')r'^{2}dr' .

For the exterior region r>R for a star of radius R, we have p=\rho=0 so the metric (-) smoothly joins into the vacuum Schwarzchild metric with total mass M=M(R)=4\pi \int_{0}^{R}\rho(r)r^{2}dr . Computing the total mass M of the star of radius R by integrating the density over its volume gives

M={\displaystyle \int_{o}^{R}}4\pi \rho(r)r^{2} (1-2GM(r)/r)^{-1/2}

The difference between this mass and the non-relativistic counterpart is

\delta M={\displaystyle \int_{o}^{R}}4\pi \rho(r)r^{2}\bigg( (1-2GM(r)/r)^{-1/2}-1\bigg)

which is due to gravitational binding energy. The interior mass term M(r), giving the mass up to given radius within the star, now “runs” within the interior region 0\le r\le R then M(r)\rightarrow 0 fast enough to remove the central singularity as r\rightarrow 0 and M(r)/r remains finite as r\rightarrow 0.

An important additional equation is the equation for hydrostatic equilibrium
which is a direct consequence of the energy conservation condition imposed on the matter source tensor \nabla^{\mu}\mathbf{T}_{\mu\nu}=0 v ia theEinstein equations. A standard computation yields the Tolman-Oppenheimer-Volkov equilibrium equation:

-r^{2}\frac{dp(r)}{dr}=GM(r)\rho(r)\left(1+\frac{p(r)}{\rho{r}}\right) \left( 1+\frac{4\pi r^{3}p(r)}{M(r)}\right)\left(1-\frac{2GM(r)}{r}\right)^{-1} 

For lower mass (Newtonian) stars–most of the stars in the night sky–the relativistic terms can be dropped giving the basic Newtonian equilibrium equation.

-r^{2}\frac{dp(r)}{dr}=GM(r)\rho(r)

This is one of the fundamental equations of stellar structure theory with general relativist corrections provided by the last three factors in the TOV  equation(ref). It is truly remarkable that this equation can be seen to arise from the purely geometrical Bianchi identify condition imposed on the Einstein tensor \nabla^{\nu}\mathbf{G}_{\mu\nu}=\nabla\mathbf{T}_{\mu\nu}=0. The relativistic factors in the Tolman-Oppenheimer-Volkoff equations are all positive definite so that for any value of r, Newtonian gravity seems to become amplified–relativistic corrections essentially strengthen the pull of gravitation.

When supplemented with an equation of state relating pressure to density,  this correctly describes a self=gravitiating spherically symmetric mass of isotropic fluid/gas in ‘gravi-hydrostatic’ equilibrium. They are integrated with respect to the boundary conditions M(0)=0, P(R)=\rho(R)=0. When its fuel is exhausted, the star can longer provide the pressure gradient required in the TOV equation so that

-r^{2}\frac{dp(r)}{dr}  \ne  GM(r)\rho(r)\left(1+\frac{p(r)}{\rho{r}}\right) \left( 1+\frac{4\pi r^{3}p(r)}{M(r)}\right)\left(1-\frac{2GM(r)}{r}\right)^{-1} 

and if above the Oppenheimer-Volkoff limit, the star collapses under its own weight to a point.  In more precise terms, a trapped surface forms around the star and null geodisic flows from the surface  are incomplete.

It is not necessary to solve the Einstein equations to prove the singularity theorems: they enter as a constraint condition imposed on matter via the strong energy conditions. Central to the proof of singularity theorems are the concepts of geodesic incompleteness, conjugate points, trapped surfaces, and the nonlinear Raychaudhuri equation for the expansion parameter of a congruence of geodesic flows. Additional singularity theorems remove some unwanted assumptions such as global hyperbolicity but the essential feature of each theorem is the concept of geodesic incompleteness in terms of the volume expansion \theta on spacetimes where the Einstein equations hold along with the strong or weak energy conditions for matter. The accepted standard definition of a singularity free spaceime is as follows:

A spacetime (\mathcal{M}, g) is singularity free if it is geodesically complete. A trapped surface \mathcal{T} is defined as follows: Let (\mathbf{\mathcal{M}},g) be a globally hyperbolic spacetime with Cauchy foliation \lbrace\mathbf{\Sigma}\rbrace . A compact 2-dimensional smooth space-like submanifold \mathbf{\mathscr{T}}\subset \Sigma , bounding a compact domain \mathbf{\mathcal{K}}\subset\mathbf{\Sigma} is a trapped surface if\theta\equiv tr(\chi)~<~0  on \mathbf{\mathcal{T}}. And \mathbf{\mathcal{T}}=\partial\mathbf{\mathcal{K}} with expansion \theta=Tr(\chi) and \chi_{\alpha}^{\beta}=\nabla_{\alpha}u^{\beta} for tangent vectors u^{\alpha}

At a trapped surface light cones will essentially “tip inward”. All trapped surfaces are entirely contained within black-hole regions of spacetime so that  \mathbf{\mathcal{T}}\in\mathbf{\mathbf{B}}, where the black-hole region \mathbf{B} is the boundary of the casual past of future-null infinity; that is \mathbf{B}=\mathbf{\mathcal{M}}-J^{-}(\mathscr{J}^{+}). The first singularity theorem was given by Penrose in 1965 in relation to trapped surfaces and future-directed null geodesics from the surface
 \bf{Theorem}
Let (\mathcal{M},\mathbf{g}) be a connected globally hyperbolic spacetime with Cauchy hypersurfaces \mathbf{\Sigma} and \mathbf{\mathcal{M}},\mathbf{g}) is a development of initial Cauchy data. let \mathbf{T}_{\mu\nu} describe matter on (\mathcal{M},\mathbf{g}) such that the Einstein equations hold. so that \mathbf{R}_{\mu\nu}=8\pi G \mathbf{T}_{\mu\nu}. For all null vectors n_{\mu}, the strong energy condition (SEC) q=\mathbf{T}_{\alpha\beta}u^{\alpha}u^{\beta}\ge 0 implies \mathbf{R}_{\alpha\beta}n^{\alpha}n^{\beta}\ge 0. Let \mathbf{\mathcal{T}} be a trapped surface, the boundary of a closed compact domain  \mathbf{\mathcal{K}}\subset\mathbf{\Sigma} with \chi_{o}=\theta_{o}<0. Then at least one inextendible future-directed orthogonal geodesic from \mathbf{\mathscr{T}} has affine length no greater than |\theta_{o}^{-1}\alpha| where \alpha=\frac{1}{2}.

There is also an initial-value Cauchy statement of the theorem in that the Cauchy development of initial data for an Einstein-matter system will always be incomplete. An important development of the theorems was the Cosmic Censure Conjecture which states thats Cauchy developments of initial data for an Einstein-matter system always leads to a singularity hidden behind a horizon or outermost trapped surface.

Let (\mathcal{M},\mathbf{g}) be a spacetime and let \mathcal{O}\subset\mathcal{M} be an open subset of \mathcal{M}. A congruence \mathcal{C} in \mathcal{O} is a set of curves or geodesics \lbrace \mathscr{G}\rbrace such that for each p\in\mathcal{O}, one and only one geodesic \mathscr{G}\subset\mathcal{C} passes through p. The associated tangent vectors for the congruence are u^{\alpha} so that u_{\alpha}u^{\alpha}=-1 for timelike geodesics and u_{\alpha}u^{\alpha}=0 for null geodesics. Defining \chi_{\alpha}^{\beta}=\nabla_{\alpha}u^{\beta} then \chi=\chi_{\alpha}^{\beta}=\theta .

In ordinary hydrodynamics, coordinates moving the flow or streamline can be chosen, then any observer moving with the stream observes rates of shear, separation, vorticity and volume variations relative to a chosen nearby flowline. By analogy an observer moving with a geodesic congruence uses a proper time s\in\mathbb{R}^{+} and has 4-velocity u^{\alpha} can also observe rates of shear, vorticity and volume variations with respect to neighboring geodesic flows.
Let \mathcal{V}\subset\Sigma be volume elements on spacelike Cauchy surfaces/slices \Sigma which are transverse to timelike geodesic congruences. Let \mathcal{A}\subset\mathcal{N} be area elements on null slices transverse to a null geodesic congruence. (Latin indices can be on these hypersurfaces):
[1] The expansion is  \theta_{ab}=h_{a}^{c}  h_{b}^{d} u_{cd}
[2] The volume expansion is \chi\equiv\theta=\theta_{ab}h^{ab}=\nabla_{a}u^{a}
[3] The shear is  \sigma_{ab}=\theta_{ab}-\frac{1}{3}h_{ab}\theta  
[4] The vorticity is  \omega_{ab}=h_{a}^{c} h_{bd}u_{cd}  

where h_{ab} is the spacelike positive-definite3-metric or 2-metric on $\latex \Sigma$ or \mathcal{N}. The trace of the extrinsic curvature of \Sigma is \chi=\chi_{\alpha}^{\alpha}=\nabla_{\alpha}u^{\alpha}\equiv\theta(s). The geodesic deviation equation or Jacobi equation gives the Raychaudhuri equation so that in terms of \chi

 \frac{d\chi(s)}{ds}=-\frac{1}{3}(\chi(s))^{2}-\mathbf{R}_{\alpha\beta}u^{\alpha}u^{\beta}-\sigma_{\alpha\beta}\sigma^{\alpha\beta}+\omega_{\alpha\beta}\omega^{\alpha\beta}

For null geodesics, with u_{\alpha}=\ell_{\alpha}, the Raychauhuri equation essentially describes geometric optics in curved space and has a similar form so that

\frac{d\chi(s)}{ds}=-\frac{1}{2}(\chi(s))^{2}-\mathbf{R}_{\alpha\beta}\ell^{\alpha}\ell^{\beta}-\sigma_{\alpha\beta}\sigma^{\alpha\beta}+\omega_{\alpha\beta}\omega^{\alpha\beta}

The energy conditions on matter require that via the Einstein equations

\mathbf{R}_{\alpha\beta}u^{\alpha}u^{\beta}=8\pi\mathbf{T}_{\alpha\beta}u^{\alpha}u^{\beta}\ge 0

The Raychuadhuri equation has very useful properties because it is a relatively simple nonlinear differential equation. It plays a central role in establishing the singularity theorems, but the Einstein-matter system itself is not solved explicitly–its only function is to impose a constraint or inequality via the energy conditions. The vorticity terms and the shear terms are usually dropped for hypersurface orthogonality so that the general basic form of the Raychaudhuri equation is the Riccati-type nonlinear differential equation

 \frac{d\chi(s)}{ds}=-\alpha(\chi(s))^{2}=q\equiv \alpha[(\chi(s))^{2}+q/\alpha]\equiv H(\chi(s))

where \alpha\in\lbrace\frac{1}{2},\frac{1}{3}\rbrace. Since q\ge 0 then

 \frac{d\chi(s)}{ds}\le -\alpha(\chi(s))^{2} 

so that for initial data \chi(0)=\chi_{o}<0, for the trace \chi_{o} of the extrinsic curvature of a trapped surface \mathcal{T} the solution blows up at s=\chi_{o}^{-1}/\alpha| or |\chi(\chi_{o}^{-1}/\alpha)|=\infty. It can be proved that this is more than a blowup in the caustic but that the spacetime is itself geodesically incomplete (refs.). For any points (p,q) a conjugate point is inflicted at q relative to p if \chi(q)=\infty.

Let \mathscr{G}\in\mathcal{C} be a geodesic of the congruence \mathcal{C}. If (p,q) \in\mathscr{G} with q\gg p then q is a conjugate point relative to p if \chi(q)=-\infty