TITLE: Proof that vector area is a boundary integral?
QUESTION [6 upvotes]: Let $M \subset \mathbb{R}^2$ be a closed topological disk and let $f: M \rightarrow \mathbb{R}^3$ be a smooth embedding; let $N$ be the corresponding unit normal field on $M$.  The vector area is defined as
$$A := \int_M N\ dA.$$
It is well-known that $A$ can also be expressed via the boundary integral
$$\frac{1}{2}\int_{\partial M} f \times df.$$
But in all the references I've seen the author merely invokes Stokes' theorem without providing further details.  (See, for instance, Sullivan, "Curvatures of Smooth and Discrete Surfaces," section 3.)
I'm not seeing it.  My question is: in excruciating detail, how do you get from point A to point B?
Thanks!

REPLY [5 votes]: Ok, finally figured it out -- the confusion basically arises from slightly nonstandard notation.  First, however, let me review some important facts about differential forms (mostly following the discussion on Wikipedia).  The application of a 2-form to a pair of vectors is often written as
$$ \alpha \wedge \beta(u,v) = \alpha(u)\beta(v) - \alpha(v)\beta(u). $$
Differential forms are often introduced in the scalar setting (e.g., $\alpha$ and $\beta$ are $\mathbb{R}$-valued) and so this expression makes sense: $\alpha(u)$ and $\beta(v)$ are both real-valued, so we can simply multiply them to get another real number.
More generally, however, suppose that $\alpha$ and $\beta$ take values in some vector space $E$.  Then there may no longer be a natural or obvious way to "multiply" vectors, and so in general the best we can do is write something like
$$ \alpha \wedge \beta(u,v) = \alpha(u) \otimes \beta(v) - \alpha(v) \otimes \beta(u), $$
where $\otimes$ denotes the outer or tensor product (which is defined for any vector space $E$).  Note, then, that the wedge product of two $E$-valued forms is not $E$-valued, but in fact $E \otimes E$ valued!
Logically, however, an element of $E \otimes E$ can be thought of as a pair of elements of $E$.  So if we happen to have some binary product $\cdot: E \times E \rightarrow E$, then we can adopt the convention that any $(E \otimes E)$-valued 2-form "automatically" gets mapped to the $E$-valued 2-form
$$ \alpha \wedge \beta(u,v) = \alpha(u) \cdot \beta(v) - \alpha(v) \cdot \beta(u). $$
In the case of the vector area, the most (only?) natural binary product on $\mathbb{R}^3$ is given by the cross product, since it takes a pair of 3-vectors to a 3-vector.  However (and here's the confusion), standard notation in this case would be
$$\alpha \wedge \beta(u,v) = \alpha(u) \times \beta(v) - \alpha(v) \times \beta(u),$$
and not
$$\alpha \times \beta(u,v) = \alpha(u) \times \beta(v) - \alpha(v) \times \beta(u).$$
(Personally I prefer to stick with the $\wedge$ notation because it is the only thing that gives the reader a clue that antisymmetrization occurs.)
So, with all that pedantic garbage out of the way, we can write a really pedantic proof of the statement from the question:
Consider the algebra bundle $E = (\mathbb{R}^3, \times)$ over a disk-like region $M \subset \mathbb{R}^2$ and let $f$ be an $E$-valued 0-form representing a smooth immersion; let $N$ be the unit normal field on $M$ induced by $f$ (also viewed as an $E$-valued 0-form).  Finally, let $dA$ be the standard volume form on $M$.  Noting that
$$df \wedge df(u,v) = df(u) \times df(v) - df(v) \times df(u) = 2 df(u) \times df(v) = 2 N dA(u,v)$$
for all pairs of vectors $u,v \in \mathbb{R}^2$ (since the cross product of two tangents points in the normal direction!), we have
$$ \int_M N dA = \frac{1}{2} \int_M df \wedge df = \frac{1}{2} \int_M d(f \wedge df) = \frac{1}{2} \int_{\partial M} f \wedge df,$$
as desired (up to a change in notation).


TITLE: Peculiarities of Cardano's Formula
QUESTION [2 upvotes]: It is easy to convert the general cubic into the form $$y^3+py+q=0$$
via the Tschirnhaus transformation. However, afterwards we can derive with a lot of work Cardano's formula: $$y=\sqrt[3]{-\frac{q}{2}+\sqrt{\frac{q^2}{4}+\frac{p^3}{27}}}+\sqrt[3]{-\frac{q}{2}-\sqrt{\frac{q^2}{4}+\frac{p^3}{27}}}$$
In my book on Galois theory by Stewart, he recalls every nonzero complex number has $3$ cube roots. If one of those cube roots is $\alpha$, then the other two are $\omega\alpha$ and $\omega^2\alpha$ where $\omega=-\frac{1}{2}+i\frac{\sqrt{3}}{2}$. This is fine. However then, Stewart says directly proceeding this: "The expression for $y$ therefore appears to lead to nine solutions, of the form $$\alpha+\beta,\space\space \alpha+\omega\beta,\space\space \alpha+\omega^2\beta$$ $$\omega\alpha+\beta,\space\space \omega\alpha+\omega\beta,\space\space \omega\alpha+\omega^2\beta$$ $$\omega^2\alpha+\beta,\space\space \omega^2\alpha+\omega\beta,\space\space \omega^2\alpha+\omega^2\beta$$
where $\alpha,\beta$ are specific choices of cube roots."
I am confused how he makes the jump from Cardano's formula and a reasoning about every $z\in\mathbb{C}\setminus\{0\}$ having $n$ distinct $n\text{-th}$ roots to this permutation argument.
Stewart further argues that not all of the above expressions are zeros (by the Fundamental Theorem of Algebra) and because we let $\space3\sqrt[3]{u}\sqrt[3]{v}+p=0$ in our derivation of Cardano's formula, we should choose $\alpha,\beta$ so that $3\alpha\beta+p=0,$ then the solutions are $$\alpha+\beta,\space\space \omega\alpha+\omega^2\beta, \space\space\omega^2\alpha+\omega\beta$$
How does Stewart arrive at this permuted list of $9$ solutions and further, how does he narrow them down to the three above?
I hope since these two questions are so very closely related in this section of the book I'll be able to ask them both in the same question here. Thanks!

REPLY [1 votes]: We obtain the roots by noting for each root $y$ some $u,\,v\in\Bbb C$ satisfy$$u+v=y,\,uv=-\frac{p}{3},$$from which we can prove$$u^3+v^3=-q,\,u^3v^3=-\frac{p^3}{27}.$$Then $u^3,\,v^3$ are the roots of$$t^2+qt-\frac{p^3}{27}=0,$$and the first root you listed is$$\sqrt[3]{u^3}+\sqrt[3]{v^3}.$$Every root is some $u$ + some $v$, i.e. of the form$$\sqrt[3]{u^3}\omega^j+\sqrt[3]{u^3}\omega^k.$$Without loss of generality we can set $j,\,k\in\{0,\,1,\,2\}$, but there are $3$ roots instead of $9$ (no surprise for a cubic!) because $uv=-\frac{p}{3}$ ensures the choice of $j$ fixes the choice of $k$.


text:
\begin{document}

\iftitle{A study of blow-ups in the Keller-Segel model of chemotaxis}

\ifauthor{Ibrahim Fatkullin}{Department of Mathematics, University of Arizona}

\ifabstract{We study the Keller-Segel model of chemotaxis and develop a composite particle-grid numerical method with adaptive time stepping which allows us to resolve and propagate singular solutions. We compare the numerical findings (in two dimensions) with analytical predictions regarding formation and interaction of singularities obtained through analysis of the stochastic differential equations associated with the model.}

\paragraph{Keywords:} Patlack, Keller-Segel, McKean-Vlasov, chemotaxis, particle method, blow-ups.

\vspace{2ex}
{\hfill\small\today}

\vspace{-2ex}
\section{Introduction}

In this work we discuss a numerical method for studying singularities (blow-ups) in the Keller-Segel model (K-S) of chemotaxis \cite{Kell70} and in the systems with similar mathematical structure, such as the McKean-Vlasov (McK-V) model of charged interacting particles \cite{Vill03}. (Patlack derived similar equations nearly two decades prior to Keller and Segel \cite{Pat53}.) We present illustrative numerical examples and provide analytical insights into formation and interaction of singularities from the point of view of the associated stochastic processes. We concentrate on the details of numerical modeling and several specific aspects of mathematical analysis. An interested reader should consult the exhaustive reviews for discussions of physical and biochemical phenomena involved in chemotaxis \cite{Hor03,Hor04,Tin08i,Tin08ii}.

Numerical treatment of equations with singular solutions is challenging because it requires a reliable and efficient approximation of singular functions or even distributions. For example, the conventional finite-difference or spectral methods are not fit for dealing with $\delta$-functions arising in the K-S model. While most finite-element and discontinuous Galerkin methods are efficient and accurate as long as the solutions remain sufficiently smooth, in the vicinity of singularities they struggle with positivity issues, unphysical oscillations, and related problems \cite{Fil06,Epsh08,Epsh09,Epsh09a,Stre10a,Stre10b}. (A conservative upwind finite-element method introduced by Saito for the elliptic K-S model overcomes some of these issues \cite{Sai07}.) Generally, unless some regularizations are introduced, these methods cannot propagate solutions past the moment of blow-up altogether.

The method presented in this work avoids such difficulties because the singular field is modeled by an ensemble of interacting particles. The singularities manifest themselves as particle aggregates, harmless from the numerical standpoint. Essentially, this approach to chemotaxis reverses the continuous PDE description of ensembles of interacting particles.

Ha\v skovec and Schmeiser have recently introduced a method utilizing similar ideas \cite{Has10,Has11}. Their method, however, is exclusively particle-based and is limited to one specific variety of the elliptic K-S model in the entire plane. Our method employs the grid and particle representations of the fields simultaneously and is suitable for both elliptic and parabolic K-S equations in arbitrary domains. It is also applicable to a wider range of kinetic problems. (Other technical differences are discussed where appropriate.) Conceptually, our method is more in the spirit of the Particle In Cell (PIC) methods widely used, e.g., in plasma physics \cite{Iva80,Cot86,Bird91,Vic91,Allen94,Hav96,Wall00}.

Another interesting numerical method has been recently implemented for the one-dimensional McKean-Vlasov system \cite{Bla08}. It uses that the McK-V equation describes steepest descent dynamics in the Wasserstein-2 space of probability measures and is suitable for studying singular (measure-valued) solutions.

The usefulness of the particle approach extends beyond numerical methods. Analysis of the stochastic processes underlying the K-S PDEs yields insights into formation and interaction of singularities. In particular, we derive an expression for the critical mass required to create or sustain a singularity and relate the non-uniqueness of solutions of the K-S equations to the underlying diffusion process. Even though our treatment is rather informal, it sheds new light on the mechanism of blow-ups.

Some of the results presented in this paper have been proven using different methods and under various assumptions. Much effort has been directed towards obtaining the exact value of the critical mass \cite{Dol03,Bla06,Calv06,Chav07,Bil09}. Vel\'azquez introduced measure-valued solutions and extended the K-S equations beyond the point of blow-up \cite{Vel04a,Vel04b}, see also \cite{Dol09}. These works together with the study of the K-S equations as a hydrodynamic limit of interacting particle systems by Stevens \cite{Stev00} link the particle and PDE descriptions of chemotaxis.

We also mention a rigorous analysis of phase transitions in the McK-V system \cite{Chay09} (it only covers sufficiently regular interaction kernels not relevant to the K-S model) and an investigation of singularity formation in several related aggregation models \cite{Bert09,Bert10,Carri10} (these models do not have diffusive terms and correspond to one special case of the K-S model).

\paragraph{The Keller-Segel model} is prescribed by the following system of PDEs:

\begin{subequations}\label{eq:KS}
\begin{align}
	\pd_t\rho(\vv x,t)\,&=\,\grd\cdot(\mu\grd\rho\,-\,\chi\,\rho\grd c)\label{eq:KSrho},\\
  	\alpha\,\pd_tc(\vv x,t)\,&=\,\Delta c\,-\,k^2c\,+\,\rho.\label{eq:KSc}
\end{align}
\end{subequations}

The function $\rho(\vv x)$ is the density of active particles (bacteria), $c(\vv x)$ is the concentration of chemoattractant. For numerical simulations we use Neumann (no flux) boundary conditions,

\begin{equation}
	\pd_{\vv n}c(\vv x)=\pd_{\vv n}\rho(\vv x)=0\quad\text{for}\quad\vv x\in\pd\dom
\end{equation}

in a two-dimensional square domain $\dom=(0,L)^2$.

We assume that $\rho(\vv x)$ integrates to $M$ over the entire domain; $M$ is the total mass of the bacteria. By rescaling the equations, we can always achieve that $\alpha=0$ or $\alpha=1$. In the former case the model is called {\em elliptic,} in the latter --- {\em parabolic.} The parameters $\mu$, $\chi$, and $k$ are constants. (The numerical method, however, may be readily extended to a more general class of equations, see Section~\ref{sec:num} for details.)

Physically, these equations describe the following phenomenon. The chemoattractant spreads diffusively and decays with rate $k^2$; it is also produced by the bacteria with rate $1$. (In the elliptic case these rates are infinite, i.e., the chemoattractant ``thermalizes'' infinitely fast.) The bacteria diffuse with mobility $\mu$ and also drift in the direction of the gradient of concentration of the chemoattractant with velocity $\chi|\grd c|$; $\chi$ is called chemosensitivity.

In this work we are interested in singular solutions to the K-S equations. A typical situation is illustrated in Figure~\ref{fig:blowups} where a snapshot of the concentration field $c(\vv x)$ is displayed. The peaks correspond to $\delta$-function-type singularities of the particle density $\rho(\vv x)$. (Displaying bacteria density itself is not very illustrative due to its highly singular nature.) As mentioned above, the principal challenge for numerical simulations is in approximating the fields $\rho(\vv x)$ and $c(\vv x)$ as they become unbounded. Our idea is to use a particle method for the evolution of $\rho(\vv x)$, so that its singularities manifest themselves as harmless particle aggregates. For the propagation of $c(\vv x)$, we use a second order implicit finite-difference scheme because of its simplicity and excellent stability properties. Even though $c(\vv x)$ also becomes unbounded analytically (it remains large but bounded for a given discretization), its weak logarithmic singularities do not cause any problems for sufficiently stable schemes. In general, one may use a different solver for propagation of the concentration field (including a particle-based solver as well).

\begin{figure}[h]
\scalebox{0.4}{\includegraphics{blowups}}\hfill
\parbox[b]{2.25in}{\caption{\label{fig:blowups}\small A snapshot ($t=10$) of the concentration field $c(\vv x)$ illustrating multiple blow-ups in the parabolic model. The parameters are as follows: $\alpha=k=1$, $\chi=0.1$, $\mu=0.005$, $M=25$. The spatial domain, $3.2\times3.2$, is discretized using the mesh size $\Delta x=0.05$; the time step for propagation of the concentration field is $\Delta t=0.1$; 4096 particles are used to approximate $\rho(\vv x)$ (see Section~\ref{sec:num} for details of the numerical scheme). Some random initial data is used at $t=0$.}}
\end{figure}

\paragraph{The particle-based description} utilizes that the PDE (\ref{eq:KSrho}) is a Fokker-Planck equation for the evolution of the probability density of the bacteria. Suppose first, that we are given some time-dependent concentration field $c(\vv x,t)$ and want to solve equation (\ref{eq:KSrho}). Let $\vv X^{(n)}_t$ be i.~i.~d. random variables distributed according to the initial condition $\rho(\vv x,0)$ at time $t=0$ and satisfying the following stochastic differential equations (with reflecting boundary at $\pd\dom$) for $t>0$:

\begin{equation}\label{eq:KS_SDE}
	\md \vv X^{(n)}_t\,=\,\chi\grd c\big(\vv X^{(n)}_t,t\big)\md t\,+\,\sqrt{2\mu}\,\md\vv W^{(n)}_t.
\end{equation}

Then, using the strong law of large numbers and that (\ref{eq:KSrho}) is a Fokker-Planck equation for the SDEs (\ref{eq:KS_SDE}), we can see that as $N\to\infty$, the empirical probability densities,

\begin{equation}\label{eq:ePDF}
	P_N(\vv x,t)\,=\,\frac{M}{N}\sum_{n=1}^N\delta\big(\vv x-\vv X^{(n)}_t\big),
\end{equation}

converge almost surely to $\rho(\vv x,t)$. (This convergence is in the sense of measures, point-wise in $t$; a stronger convergence can be established with more effort.)

Thus we can simulate equations (\ref{eq:KS_SDE}) for the particles, and approximate the density $\rho(\vv x,t)$ using $P_N(\vv x,t)$ whenever required.

Rigorously justifying this approach to the full K-S system (\ref{eq:KS}) is more challenging because in the full system the concentration field $c(\vv x,t)$ is a functional of $\rho(\vv x, s)$, $s\leq t$, rather than some a priori prescribed field. This implies that equations (\ref{eq:KS_SDE}) become coupled and we end up with a system of stochastic interacting particles (with memory in the parabolic case). Some progress in this direction, nevertheless, has been made.  Stevens derived the K-S equations as limit dynamics for interacting stochastic particle systems via smoothing and rescaling the interaction potentials in a particular fashion \cite{Stev00}. Ha\v skovec and Schmeiser proved convergence of their particle method under reasonable assumptions via analysis of the associated BBGKY hierarchy \cite{Has11}.

\section{Numerical method}\label{sec:num}

As mentioned in the introduction, the principal idea behind our numerical scheme is to employ a particle method for the evolution of the particle density field $\rho(\vv x,t)$. The current implementation utilizes an implicit second order finite-difference scheme for simulating equation (\ref{eq:KSc}) and an explicit Euler-Maruyama scheme for the stochastic particle dynamics (\ref{eq:KS_SDE}). We discretize the computational (rectangular) domain with grid size $\Delta x$ and propagate the concentration field $C_{ij}$ using the time step $\Delta t$.  The particle density field $P_{ij}$ is reconstructed from the particle locations. We evolve the particles using adaptive time steps which may be smaller than $\Delta t$; this is needed for stability reasons (see below).

\paragraph{The particle dynamics}

is simulated using the forward Euler-Maruyama scheme,

\begin{equation}
  \vv X^{(n)}(t+\Delta\tau)\,=\,\vv X^{(n)}(t)\,+\,\chi\grd c\big(\vv X^{(n)}(t)\big)\Delta\tau\,+\,\sqrt{2\mu\,\Delta\tau}\,\vv N(0,1),
\end{equation}

where $N(0,1)$ is a standard Gaussian random variable with mean 0 and variance 1. The gradient field $\grd c(\vv x)$ is approximated in two steps. First, we construct the gradient fields $CX$ and $CY$ using the second order approximation,

\begin{equation}
  CX_{ij}\,=\,\frac{1}{2\Delta x}\Big[C_{i+1,j}(t)\,-\,C_{i-1,j}(t)\Big],\qquad
  CY_{ij}\,=\,\frac{1}{2\Delta x}\Big[C_{i,j+1}(t)\,-\,C_{i,j-1}(t)\Big].
\end{equation}

Then we approximate $\grd c(\vv x)$ via bilinear interpolation using the values of $CX$ and $CY$ at the four nearest grid points. Because $|\grd c(\vv x)|$ becomes unbounded (very large numerically) in the vicinity of blow-ups, it is essential to choose the time step $\Delta\tau$ adaptively, see Figure~\ref{fig:adapt} and its description for details. Consequently, if needed, we subdivide the timestep $\Delta t$ into smaller intervals $\Delta\tau^{(n)}$ so that $|\grd c\big(\vv X^{(n)}\big)|\Delta\tau^{(n)}<\Delta x$. (Each particle is simulated independently of others, i.e, $\Delta\tau^{(n)}$ are chosen for each particle individually.) A failure to satisfy this condition leads to overshooting in the vicinity of blow-ups, which in turn causes artificial oscillations and damping. Finally, we enforce the no-flux Neumann  boundary conditions for $\rho(\vv x)$ by reflecting the particles escaping from the spatial domain back into it.

\begin{figure}
\scalebox{0.7}{\includegraphics{jump}}\hfill
\parbox[b]{3.75in}{\caption{\label{fig:adapt}\small Adaptive time stepping for particle dynamics. In this example a blow-up of the concentration field occurs near a grid point with index $j$. While some particular time step is appropriate for propagating the particles away from the blow-up (i), it may lead to overshooting in its vicinity (ii, iii). Thus the particles incorrectly accumulate near the grid points $j\pm1$ rather than $j$. One can avoid this phenomenon reducing the time step for the particle evolution whenever the gradient of the concentration field becomes too large. Practically, it is sufficient to adjust the time step so that the expected length of the particle jump does not exceed the mesh size $\Delta x$.}}
\end{figure}

\paragraph{The particle density field $P_{ij}(t)$} is reconstructed from the particle locations as explained in the description to Figure~\ref{fig:scheme}: each particle contributes fractions of its weight to the four nearest grid points according to the bilinear interpolation rules.  This approximation preserves the first moment of the corresponding distribution. It is important to use this kind of moment-preserving approximations, e.g., a simple bin-counting (assigning all weight of a given particle to the nearest grid point) is unsatisfactory. The reason is that this creates an artificial flux towards the grid points: the particles feel their own potentials which become artificially aligned to the nearest grid points. This flux is sufficient to pin the singularities to grid points and disrupt such phenomena as logarithmically-weak interaction of singularities with each other and with the boundary of the domain.

\paragraph{The concentration field} is propagated according to the implicit second order finite-difference\linebreak scheme,

\begin{equation}
  \frac{\alpha}{\Delta t}\,\Big[C_{ij}(t+\Delta t)\,-\,C_{ij}(t)\Big]\,=\,\frac{1}{\Delta x^2}\,D^{(2)}_{ij}C(t+\Delta t)\,-\,k^2C_{ij}(t+\Delta t)\,+\,P_{ij}(t),
\end{equation}

where $D^{(2)}$ is the standard second difference operator. It is essential to use an implicit scheme for stability reasons because the concentration field acquires logarithmic singularities. This scheme is quite efficient: even though it requires solving a system of linear equations, the matrix of the finite difference operator is symmetric and banded, and may be Cholesky-factorized  before the actual computations. The extra computational cost of solving this linear system is then alleviated by relaxation of the $\Delta t\sim\Delta x^2$ stability constraint of explicit schemes.

\begin{figure}
\scalebox{0.7}{\includegraphics{mesh}}\hfill
\parbox[b]{3.5in}{\caption{\label{fig:scheme}\small Schematic representation of the numerical scheme. The concentration field $c(\vv x)$ is sampled on a uniform grid ($\blacksquare$), while the particles ($\bigstar$) move without any restraints. The particle density $\rho(\vv x)$ is reconstructed via bilinear interpolation: each particle contributes fractions of its mass to the four nearest grid points proportionally to the relative distances from the latter, e.g., a particle located at $(x,y)$ within a unit-square grid cell contributes weights proportional to $xy$, $x(1-y)$, $y(1-x)$, and $(1-x)(1-y)$ to the respective grid points.}\rule{0pt}{2ex}}
\end{figure}

\newpage
The presented numerical method may be generalized to other similar systems of the form

\begin{subequations}\label{eq:KSgen}
\begin{align}
  \pd_t\rho(\vv x,t)\,&=\,\grd\cdot\big(\mu(\vv x,\rho,c)\grd\rho\,-\,\rho\,\vv v(\vv x,\rho,c)\big)\label{eq:KSrhogen},\\
  \alpha\,\pd_tc(\vv x,t)\,&=\,\grd\cdot\big(G(\vv x)\grd c\big)\,+\,F(\vv x,\rho,c),\label{eq:KScgen}
\end{align}
\end{subequations}

and to more complicated domains. This makes it a useful tool for studying various aggregation phenomena for which the conventional finite-difference or finite-elements schemes are inefficient. From the numerical analysis point of view, we do not currently have any estimates of the convergence rates for this scheme. Partially, this problem is complicated due to the singular nature of solutions to the K-S equations. In particular, once the blow-ups are formed, one must interpret the K-S PDEs in some proper sense, and while our numerical scheme is able to propagate the solutions for all times (i.e., it regularizes the PDEs in some fashion), relation of such numerical solutions to a particular analytical regularization must be carefully investigated.

On the other hand, we could speculate that our method exhibits the typical for this kind of schemes errors of the order $\bO(\sqrt{\Delta t})$, $\bO(1/\sqrt{N})$, and $\bO(\Delta x^2)$ for regular solutions. Such analysis, however, is not the purpose of this paper and we verify our numerical results by ``standard means,'' i.e., compare the results with simulations employing half the mesh size, double the number of particles, etc.

\section{Formation and interaction of singularities}

We now concentrate on the elliptic ($\alpha=0$) case, though some of the reasoning is equally applicable to the parabolic case as well. Generally, one can distinguish the {\em soft\,} and {\em hard\,} blow-ups: in the first case the fields $\rho(\vv x)$ and $c(\vv x)$ become unbounded but $\rho(\vv x)$ does not acquire atomic components, while in the second case the particle density field becomes a bona fide distribution. Here we concentrate on the hard blow-ups which have a more apparent physical meaning: atomic components of the particle density correspond to accumulation  of particles at some point locations on the given physical length-scale.

It is convenient to study the elliptic model in a slightly different form, as a McKean-Vlasov system. First, we solve equation (\ref{eq:KSc}) for $c(\vv x)$:

\begin{equation}\label{eq:cconv}
  c(\vv x)\,=\,-\int_\dom V_d(\vv x,\vv y)\,\rho(\vv y)\md\vv y\,=\,-\,V_d*\rho(\vv x).
\end{equation}

Here $V_d(\vv x,\vv y)$ is the Green's function for $\Delta-k^2$ in $\dom$, which is assumed to be $d$-dimensional. The entire-space Green's functions (the fundamental solutions)  whose singular part is identical to those in bounded domains are known explicitly:

\begin{subequations}\label{eq:FS}
\begin{align}
	V_1(\vv x,\vv y)\,=&\,-\frac{\me^{-k|\vv x-\vv y|}}{2k}&(d=1);\\
	V_2(\vv x,\vv y)\,=&\,-\frac{\kbes_0(k|\vv x-\vv y|)}{2\pi}&(d=2);\label{eq:fs_2d}\\
	V_3(\vv x,\vv y)\,=&\,-\frac{\me^{-k|\vv x-\vv y|}}{4\pi|\vv x-\vv y|}&(d=3).
\end{align}
\end{subequations}

Here $\kbes_0(\cdot)$ is the modified Bessel function of second kind. Substituting expression (\ref{eq:cconv}) for the concentration field $c(\vv x)$ into equation (\ref{eq:KSrho}) we get a closed integro-differential equation for the particle density field,

\begin{equation}\label{eq:KS_ell_rho}
  \pd_t\rho(\vv x,t)\,=\,\grd\cdot(\mu\grd\rho\,+\,\chi\,\rho\grd V_d*\rho)\,=\,\grd\cdot\left[\rho\grd\VD{\eNrg}{\rho}\right].
\end{equation}

Here the energy functional is given by

\begin{equation}\label{eq:nrg}
  \eNrg(\rho)\,=\,\mu\int_\dom\rho(\vv x)\ln\rho(\vv x)\md\vv x\,+\,\frac{\chi}{2}\iint_\dom\rho(\vv x)V_d(\vv x,\vv y)\rho(\vv y)\md\vv x\md\vv y.
\end{equation}

One of the features of the dynamics (\ref{eq:KS_ell_rho}) is that the functional $\eNrg(\rho)$ is non-increasing:

\begin{equation}
  \PD{\eNrg}{t}\,=\,-\int_{\dom}\left|\grd\VD{\eNrg}{\rho}\right|^2\rho(\vv x,t)\md\vv x\,\leq\,0.
\end{equation}

When $\rho(x)$ becomes singular with respect to Lebesgue measure, the first ({\em entropic}) term in (\ref{eq:nrg}) tends to infinity, which must be compensated by the second ({\em interaction}) term. If $d=1$, the Green's function $V_1(\vv x)$ is bounded below and so is the interaction term, thus the blow-ups are not possible. This is no longer the case if $d\geq 2$ and so the blow-ups are permissible. A more delicate analysis is needed to understand how exactly they are formed.

From standpoint of the stochastic process underlying the Fokker-Planck PDE (\ref{eq:KSrho}), the atomic components of $\rho(\vv x)$ correspond to particles aggregating at point locations. Thus the traps  which do not allow the particles to escape must be created. Let us first investigate when a particle aggregate creates a trap for a single particle diffusing in its field. Suppose $\rho(\vv x)=M\delta(\vv x)$, the corresponding concentration field is then proportional to the fundamental solution (\ref{eq:FS}).  Equation for the radial component of the location of a single particle diffusing in the field created by $\rho$ is equivalent in law to the following SDE:

\begin{equation}
  \md r_t\,=\,\Big[\chi c^\prime\big(r_t\big)\,+\,(d-1)\frac{\mu}{r_t}\Big]\md t\,+\,\sqrt{2\mu}\,\md W_t.
\end{equation}

For $d=1$ we get

\begin{equation}
  \md r_t\,=\,-\frac{\chi M}{2}\exp(-kr_t)\md t\,+\,\sqrt{2\mu}\,\md W_t.
\end{equation}

This equation is sufficiently regular and does not allow for existence of a trap for any value of $M$, so consistently with energy considerations, no blow-ups are possible in this case. For $d\geq 3$, a similar reasoning shows that a trap of an arbitrarily small mass may exist. So the interesting case is $d=2$.

For $d=2$, using that in the leading order, as $r\to0$, $\kbes_0(r)\sim-\ln(r)$, we obtain the following equation:

\begin{equation}\label{eq:bessel}
  \md r_t\,=\,\Big(\mu\,-\,\frac{\chi M}{2\pi}\Big)\frac{\md t}{r_t}\,+\,\sqrt{2\mu}\,\md W_t.
\end{equation}

This is a well-known Bessel process \cite{Kar81}, it's behavior near the boundary at $r=0$ depends strongly on the value of $M$. In particular, when

\begin{equation}
  M\,\geq\,\frac{4\pi\mu}{\chi}\,=\,M_c^*,
\end{equation}

the origin is the so-called {\em exit boundary,} a trap, i.e., the particles reach it in finite time, and may never escape back to the domain $r>0$. For our problem the implication is that the total mass of at least $M_c^*$ is needed for existence of a stable singularity in the particle density field. For $M\in(0,M^*_c)$, the SDE (\ref{eq:bessel}) has the so-called {\em regular boundary} at $r=0$. In this case the particles may reach the origin, but are not necessarily trapped there, i.e., may also leave it according to some rules which must be prescribed in addition to the SDE itself. This implies that the original K-S PDEs alone are not sufficient to describe the blow-up dynamics and the exchange of mass between the regular and singular components of the particle density when the mass of the singularity is less than $M_c^*$.

The critical mass $M_c^*$ is smallest mass such that a singularity does not shed its mass but only absorbs particles from the smooth component. It is, however, twice less than the mass required to create a singularity from smooth initial data. Indeed, consider the $N$-particle SDEs approximating the elliptic K-S equations:

\begin{equation}\label{eq:KS_SDE_n}
  \md \vv X^{(n)}_t\,=\,-\frac{\chi M}{N}\Pd{}{\vv X^{(n)}}\sum_{i=1,\;i\neq n}^N
  V_d\big(\vv X^{(n)}_t\!\!,\,\vv X^{(i)}_t\big)\md t\,+\,\sqrt{2\mu}\,\md\vv W^{(n)}_t.
\end{equation}

To study the aggregation dynamics of the particle system, consider the variance of the empirical probability density of the bacteria (\ref{eq:ePDF}),

\begin{equation}\label{eq:R}
	R_t^2\,=\,\frac{1}{\,2N^2}\sum_{i,\,j=1}^{N}\big|\vv X^{(i)}_t-\,\vv X^{(j)}_t\big|^2;
\end{equation}

$R_t$ will be referred to as the\, {\em radius\,} of the system. We concentrate on the problem in the entire space with $k=0$; this corresponds to $V_2(\vv x,\vv y)=1/(2\pi)\ln|\vv x-\vv y|$. Analysis for the other values of $k$, or in the bounded domains is similar (although more technical) because the singularities of $V_2(\vv x,\vv y)$  are the same in the leading order. In this case, however, there exists an explicit closed equation governing the evolution of $R_t$:

\begin{equation}\label{eq:SDE_N}
	\md R_t\;=\,\Bigg[2\mu\,\bigg(1-\frac{3}{2N}\bigg)-\,\frac{\chi M}{4\pi}\bigg(1-\frac{1}{N}\bigg)\Bigg]\;\frac{\md t}{R_t}\;+\;\sqrt{\frac{2\mu}{N}\;}\md W_t.
\end{equation}

The driving stochastic process $W_t$ obeys the following SDE:

\begin{equation}
	\md W_t\,=\,\frac{1}{N^{3/2}\,R_t}\;\sum_{i,\,j=1}^{N}\Big(\vv X^{(i)}_t-\,\vv X^{(j)}_t\Big)\cdot\,\md\vv W^{(i)}_t,
\end{equation}

and can be shown to be a Wiener process using the L\'evy characterization of Brownian motion. We can immediately see that in the limit as $N\to\infty$, if

\begin{equation}
 	M\,\neq\,\frac{8\pi\mu}{\chi}\,=\,M_c,
\end{equation}

the evolution of the radius obeys a well-known deterministic equation,

\begin{equation}
	\dot R_t\,=\,-\frac{\gamma}{2R_t},\qquad
	\gamma\,=\,4\mu(M/M_c-1).
\end{equation}

Hence $R^2_t=R^2_0-\gamma t\,$, and for supercritical masses ($M>M_c$) the blow-up occurs at time $T=R_0^2/\gamma$. It is independent of the initial distribution of mass except through $R_0$. If $M=M_c$, after the time rescaling, equation (\ref{eq:SDE_N}) may be rewritten as

\begin{equation}\label{eq:SDE_N_crit}
	\md R_{t^\prime}\,=\,\frac{1}{2}\frac{\md t^\prime}{R_{t^\prime}}\,+\,\md W_{t^\prime},
	\qquad t^\prime=\frac{2\mu t}{N}.
\end{equation}

This equation describes a critical Bessel process with {\em entrance boundary} at $R=0$, in this case $R_{t^\prime}$ a.s. never hits zero, although it visits its arbitrarily small neighborhood. In this limit the particle system remains stochastic even as $N\to\infty$ and does not approximate the deterministic K-S equations. As $t^\prime\to\infty$, the particle system approaches blow-up infinitely often, but never actually coalesces into a single particle, cf \cite{Bla08b}.

\begin{figure}[h]
\scalebox{0.4}{\includegraphics{sub}}\hfill\scalebox{0.4}{\includegraphics{super}}\\
\caption{\label{fig:supsub}\small Snapshots ($t=10000$) of the concentration field $c(\vv x)-M/(kL)^2$ demonstrating sub- and supercritical behaviors in the elliptic ($\alpha=0$) model. (The constant is subtracted because $c(\vv x)$ itself becomes unbounded in a finite domain as $k\to0$.) The other parameters are as follows: $k=0.01$, $\chi=0.1$, $\mu=0.005$; $\Delta t=0.01$, $\Delta x=0.05$, $4\cdot10^3$ particles. Initially, all particles are distributed randomly with a slight bias towards $(0,0)$. In the supercritical case ($M=0.35$, right plot), a trap is created, while in the subcritical case ($M=0.34$, left plot), all particles eventually spread uniformly over the entire domain. The theoretical critical mass for the elliptic model  with $k=0$ in the entire space is $M_c=2\pi/5\approx1.26$. Notice that only a quarter of the critical mass is needed for a trap in the corner.}
\end{figure}

\paragraph{Formation of singularities.} Our first numerical experiment is designed to test how well our numerical method predicts the value of the critical mass $M_c$ required to develop a singularity from smooth initial data. Notice, first of all, that due to reflection principle for the problem with Neumann boundary conditions, singularity in a corner of the domain, $\Omega=(0,L)^2$, only requires a quarter of the critical mass.  Since the singularity inside of the domain eventually migrates into a corner, we perform our experiment in the corner to begin with. We scatter all particles over $\Omega$ with a slight bias towards $(0,0)$ and vary their total mass observing whether the particles remain aggregated, or spread uniformly over the entire domain. The results are presented in Figure~\ref{fig:supsub}. For the specified values of the parameters, the calculated numerical value of the critical mass lies between 0.34 and 0.35, while the theoretical prediction for elliptic model with $k=0$ in the entire plane is $M_c/4=\pi/10$ --- slightly smaller. The mismatch is a consequence of the simulation in a finite domain, the limit $L\to\infty$, $k\to0$ is non-trivial because equation (\ref{eq:KSc}) with Neumann boundary conditions is not well-posed in $\Omega$ when $\alpha=k=0$. In particular, the boundary of the domain far from the corner where the singularity is formed pulls the particles away from the singularity, effectively increasing the critical mass.

Interesting phenomena occur when a singularity with mass $M\in(0,M^*_c)$ already exists in the system. In this case the trap absorbs the particles, though the latter may still escape back into the regular component of the particle density. The diffusion process underlying the K-S equations is not uniquely defined by its generator and additional rules for behavior of the particles at singular points must be specified. These rules are not inherently encoded in the K-S equations and are related to non-uniqueness of K-S regularizations.  From the modeling perspective, the exchange of mass between the regular and singular components of the particle density in this regime strongly depends on particular details of the numerical method. Note, however, that if the initial particle density is regular, the smallest possible mass for singularity is greater than $M_c=2M_c^*$,  so this scenario does not occur.

\rule{0pt}{4ex}
\begin{figure}[h]
\center\scalebox{0.225}{\includegraphics{evolv}}\\
{\caption{\label{fig:evolv}\small Snapshots of the concentration field $c(\vv x)$ illustrating the motion and interaction of singularities in the elliptic ($\alpha=0$) model. The other parameters are as follows: $k=1$, $\chi=0.1$, $\mu=0.005$, $M=4$; $\Delta t=0.01$, $4\cdot10^3$ particles. Initially, the particles are placed around four distinct locations; as time goes on, the singularities are formed, then they merge and travel towards the boundary of the domain, and finally stabilize in the corners. (Notice the scale change for the $z$-axis.)}\rule{0ex}{1ex}}
\end{figure}

\paragraph{Interaction of singularities} is illustrated in Figure~\ref{fig:evolv}. The singularities attract each other and are also attracted by the boundary of the domain. Equations governing their dynamics may be derived following \cite{Vel04a,Vel04b,Dol09}. If the density $\rho(\vv x)$ is purely atomic, we obtain

\begin{equation}\label{eq:sing_ODE}
  \dot{\vv x}_i(t)\,=\,-\chi\,\frac{\pd}{\pd\!\vv x_i}\sum_{j\neq i}m_j\,V_2(\vv x_i,\vv x_j).
\end{equation}

Here $m_j$-s are the masses of singularities located at $\vv x_j$-s, $V_2(\vv x,\vv y)$ is the Green's function. Note that unlike in the method of Ha\v skovec and Schmeiser \cite{Has10}, in our method, the dynamics of singularities is not imposed by the numerical method explicitly, i.e., it is a natural consequence of the stochastic particle dynamics.

Equations (\ref{eq:sing_ODE}) are identical to equations (\ref{eq:KS_SDE}) of the K-S model with particle diffusivity $\mu$ set to zero. The naive reasoning is that once the traps are formed, the particles' diffusion is dominated by the (infinitely strong) drift and ceases to contribute into dynamics: the stochastic dynamics self-averages and the martingale component does not contribute into the mean drift. In order to verify how well our method approximates equations (\ref{eq:sing_ODE}), we perform a numerical experiment illustrated in Figure~\ref{fig:coll}. In this experiment we create two singularities far from the boundaries in a sufficiently large domain, so that the interaction potential is well-approximated by the fundamental solution (\ref{eq:fs_2d}). We run the K-S solver and track locations of these singularities. Finally, we compare them with locations of the particles evolving according to equation (\ref{eq:sing_ODE}) from the same initial data. A perfect match indicates that the scheme is very successful in dealing with this kind of phenomena. Similarly, in a bounded domain, the particle interaction potential, $V_2(\vv x_i-\vv x_j)$, should be replaced by the Green's function which also encodes interaction with the boundary. In particular, this explains attraction of the singularities to the boundary and the corners of the domain.

\begin{figure}[h]
\rule{0ex}{4ex}\scalebox{0.25}{\includegraphics{coll}}\hfill
\parbox[b]{3.25in}{\caption{\label{fig:coll}\small Interaction of singularities in the elliptic ($\alpha=0$) model. The particle density $\rho(\vv x)$ is initialized with two delta-functions with masses 6.25 and 18.75. These singularities attract each other and eventually merge; their locations are plotted as functions of time $t$. The thin black line corresponds to the simulation of equation (\ref{eq:sing_ODE}), while the thick grey line --- to simulation of the elliptic Keller-Segel model (the lines are barely distinguishable); $\Delta t=0.01$, $4\cdot10^5$ particles are used. The computational domain is chosen sufficiently large so that the boundary effects are negligible on the presented time scale.}}
\end{figure}

\section{Discussion}

Our investigation of the K-S equations was motivated by the underlying stochastic particle dynamics. A numerical method based on similar ideas has been recently presented by Ha\v skovec and Schmeiser \cite{Has10}.  The greatest advantage of the H-S method is that the singularities are represented as deterministic particles rather than as a tight cluster of ``elemetary'' stochastic particles. This, however, requires an explicit knowledge of evolution equations governing the interaction of singularities (at the very least the knowledge of Green's function in a given domain) and is only applicable to the elliptic K-S model. Moreover, the H-S method requires a direct simulation of an ensemble of interacting particles and thus its computational cost scales as $N^2$ with respect to the total number of particles. By combining the particle dynamics and the PDE dynamics of the concentration field which mediates the particle interaction (the PIC ideas) we are able to significantly reduce this computational cost and deal with a much larger number of interacting particles. We can also treat both parabolic and elliptic models with equal ease. The ultimate numerical method should combine ideas presented in this work, the deterministic particles of the H-S method, and also high order conventional schemes in the regions where the solution remains regular. Such a multi-model implementation remains a subject for the future work.

Even though the non-uniqueness issue is well-known, the details of how the exit-entrance conditions for the underlying diffusion process are related to the exchange of mass between the singularity and the regular part  are not yet well understood. From the modeling perspective, the implication is that a particular numerical scheme provides a regularization of some sort which affects the critical mass $M_c^*$

Another interesting class of stochastic particle models related to our studies of the K-S equations arises in the context of the so-called self-gravitating Brownian particles, see e.g., \cite{Chav07,Chav07b,Sir10}. We suggest a model which fits into this class and bears an intrinsic relation to the K-S model. Consider an ensemble of particles in $d$ dimensions characterized by their masses $m_n$ and evolving according to the stochastic differential equations

\begin{equation}\label{eq:aggr}
  \md\vv X^{(n)}_t\,=\,-\chi\,\frac{\pd}{\pd\!\vv X^{(n)}}\sum_{k\neq n}m_kV_d\big(\vv X^{(n)}_t\!,\,\vv X^{(k)}_t\big)\,\md t\,+\,\sqrt{2\mu M/N_0m_n}\,\md\vv W^{(n)}_t.
\end{equation}

Here $N_0$ is the initial number of particles and $M=\sum m_n$ is their total mass. If initially $m_n=M/N_0$, equation (\ref{eq:aggr}) is precisely (\ref{eq:KS_SDE_n}). The  interaction allows for particle collisions, thus equation (\ref{eq:aggr}) is only valid until the first collision, at which point the colliding particles coalesce into a single particle which acquires their combined mass. Dynamics is then restarted with the remaining particles. In the limit as $N_0\to\infty$, whenever a particle accumulates an $\bO(1)$ mass, its diffusivity becomes zero, i.e., it obeys the deterministic equation (\ref{eq:sing_ODE}) which governs evolution of a point singularity in the K-S model. Therefore we conjecture that in a proper hydrodynamic limit this system is equivalent to the elliptic K-S model. Investigation of this system will be conducted elsewhere.

\bibliographystyle{plain}

\bibliography{../bibliography/bibl}

\end{document}

config: arxiv


text:
\begin{document}
\maketitle

\section*{Introduction}
Let $C$ be a non-singular curve of genus $g$ over an algebraically
closed field ${\bf k}$ of characteristic zero. Let $K$ be the canonical
sheaf on $C$. If $C$ is not hyperelliptic, the map associated to the
 complete canonical series $|K|$
$$C\rightarrow {\bf P}^{g-1}$$
is an embedding and the image curve is projectively normal. If the
 curve is neither trigonal nor a plane quintic, the ideal of $C$
is known to be generated by quadrics. Continuing in this vein,
Mark Green made the conjecture that the resolution of the ideal
of $C$ in ${\bf P}^{g-1}$ should depend on the linear series that $C$
has.

To make this precise, one defines property $N_p$.
Take a minimal resolution of the ideal sheaf of $C$ in ${\bf P}^{g-1}$.
Then one says that  property $N_0$ holds if $C$ is projectively
normal, $N_1$ means that the ideal of the curve is generated
by quadrics, $N_2$ means that in addition the syzygies among these quadrics
are generated by linear relations...In general $N_p$ means that $N_{p-1}$
holds and the $p^{th}$ syzygies are generated by linear relations.

Define the Clifford index of $C$ by
$$Cliff(C)=min\{ degL-2(h^0(C,L)-1)|L\in Pic C, h^0(C,L)\ge 2,h^1(C,L)\ge 2\}$$
In particular, for a curve that is generic in the sense of moduli, the
Clifford index of $C$ is given by $[(g-1)/2]$ while the most special curves
from the point of view of Clifford index are hyperelliptic curves that have
$Cliff(C)=0$.
Then
\begin{Conj}
\label{conj}
[Green](cf. \cite{G} 5.1)
 The curve $C$ has property $N_p$ if and only if
$Cliff (C)>p$.
\end{Conj}

The only if part of \ref{conj} was proved by Green and Lazarsfeld
 (cf. \cite{G} Appendix),
The conjecture has been proved for $g\le 8$ (cf \cite{S1}) and
for $p=2$ (i.e. $C$ quatrigonal) (cf \cite{V,S2,HPR}).

In a slightly more modest vein
\begin{Conj}
\label{conjgen}[Generic Green's conjecture]
(cf. \cite{G} 5.6)
 The generic curve $C$ of genus $g$ satisfies $N_{[(g-3)/2]}$.
\end{Conj}

Few effective results seem to be available with respect to \ref{conjgen}.
The only published work that we are aware of is \cite{E}
where L.Ein  showed that the first two steps of
the resolution of a generic canonical curve are as expected.
The relevance of \ref{conjgen} is enhaced by the following result of
Hirschowitz and Ramanan (cf. \cite{HR} Theorem 1.1).

\begin{Thm}[Hirschowitz-Ramanan]
\label{HR}
For odd $g=2k+1\ge 5$, Green's Conjecture holds for the generic
curve if and only if it holds for all curves of (maximal)
Clifford index $k$.
\end{Thm}

The purpose of this paper is to prove the generic Green's Conjecture.

\begin{Thm}
\label{teorema}
Let $C$ be a generic curve of genus $g$. Then, Green's conjecture on
 the syzygies of  the canonical curve $C$ holds.
\end{Thm}

As a consequence of Hirschowitz and Ramanan's Theorem,
 one obtains a very concrete open dense subset of ${\cal M}_g$
for odd $g$ where Green's Conjecture actually holds. In the language of
Hirschowitz and Ramanan, this is the specific Green's Conjecture

\begin{Cor}
If $g=2k+1\ge 5$, then Green's Conjecture holds for all curves of
maximal Clifford index $k$ (i.e. $C$ satisfies $N_{k-1}$ if and only if
$Cliff(C)=k$).
\end{Cor}

The proof of \ref{teorema} is as follows: Define a vector bundle $E$
 as the dual of the kernel of the evaluation map of the canonical
linear series. Namely, $E^*$ is defined by the exact
sequence
\begin{Def}
\label{E}
$$0\rightarrow E^*\rightarrow (H^0(C,K_C))^*
\otimes {\cal O}_C\rightarrow K \rightarrow 0.$$
\end{Def}

From a result of Paranjape-Ramanan (cf. \cite{PR} Remark 2.8, p.507),
Green's Conjecture would follow from the surjectivity of the maps
$$\wedge ^r(H^0(C,K_C))^* \rightarrow H^0(\wedge ^r E), r\le Cliff(C).$$
 We want to prove that
 this is the case for $C$ generic. To this end consider a hyperelliptic
curve $C_0$. Notice that the map
$$(H^0(C,K_{C_0}))^*\otimes {\cal O}_{C_0}\rightarrow E \rightarrow 0$$
identifies $(H^0(C,K_{C_0}))^*  $ to a subspace $W$ of $H^0(C_0,E)$.
 As $C_0$ is
 hyperelliptic, this subspace is proper. We shall start by computing
$W$ and the image $W^r\subset H^0(C_0,\wedge ^rE)$ of its exterior
 powers $\wedge ^r W$. Every infinitessimal
 defformation of the curve $C_0$
 determines a unique  infinitessimal defformation of $E$
that preserves it as the dual of the kernel of the canonical evaluation
map. Consider the defformation of $\wedge ^r E$ that this induces.
 We then see that the only sections of
$H^0(\wedge ^rE)$ that give rise to sections of the infinitessimal
defformation of $\wedge ^rE$ are those in $W^r$. This will conclude
the proof.

Acnowledgments: I would like to thank the following people and institutions
that contributed in different ways to the present work.
My interest on Green's Conjecture came from a conversation with S.Ramanan.
 Mark
Green read a first version of the paper and made several
suggestions to improve the presentation.
I am visiting the Pure Mathematics Department at the University of
Cambridge and I benefitted from conversations with Tony Wasserman and
Nick Shepherd-Barrow as well as e-mails with Loring Tu.
 I am a member of the Europroj group
Vector Bundles on Algebraic Curves.

\section{Identification of the vector bundle $E$ and its space of sections}
{\bf Notations} In this section and section 3, $C=C_0$ denotes a hyperelliptic
 curve  and $L$ the unique line bundle of degree
 two with two sections.

If $F$ is a vector bundle on a curve $C$, we write $H^i(F)$ for $H^i(C,F)$
if there is no danger of confusion.
 If several curves are involved or if we are
considering sections on an open set only, we shall make this clear.

As in the previous section $E$ will denote the vector bundle defined
in \ref{E}, $W=W^1$ will be the image of $H^0(K)^*\subset H^0(E)$ and
$W^r$ the image of the exterior powers of $W$ in $H^0(\wedge^r E)$.

\begin{Prop}
\label{W}
 Let $\pi$ be the map $C_0\rightarrow {\bf P}^1$ associated to
$L$. Let $H$ be the space of sections of $L$. Then,
$$E\cong S^{g-2}H^*\otimes \wedge ^2H^*\otimes L$$
Moreover, the map $(H^0(K))^*\rightarrow H^0(E)$ can be identified to the
natural inclusion
$$\varphi _1 :S^{g-1}H^*\rightarrow S^{g-2}H^*\otimes \wedge ^2H^*\otimes H$$
\end{Prop}
\begin{pf}
Note that the canonical sheaf $K_{C}$ on $C$ is of the form
$$K=L^{\otimes g-1}=\pi ^*({\cal O}_{\bf P^1}(g-1))$$
and
 $$H^0(C,K_C)=S^{g-1}(H)
=\pi ^*(H^0({\bf P}^1,{\cal O}_{\bf P^1}(g-1)))$$
Therefore, the exact sequence defining $E^*$ is the pull-back of
 the exact sequence in ${\bf P}^1$
\begin{equation}
\tag{*}0\rightarrow Ker \rightarrow
H^0({\bf P}^1,{\cal O}_{\bf P^1}(g-1))\otimes {\cal O}_{{\bf P}^1}
\rightarrow {\cal O}_{\bf P^1}(g-1))\rightarrow 0 \end{equation}

Denote by $\bar H$ the space of sections $H^0({\bf P^1},
 {\cal O}_{{\bf P^1}}(1))$.
Tensoring (*) with ${\cal O}(1)$ and taking global sections,
one obtains
$$0\rightarrow H^0(Ker\otimes {\cal O}(1)) \rightarrow S^{g-1}\bar H\otimes
 \bar H \rightarrow S^g \bar H\rightarrow 0$$
Hence, $H^0(Ker\otimes {\cal O}(1))\equiv S^{g-2}\bar H\otimes
 \wedge ^2\bar H$ (cf. for example \cite{FH}  p.224 15.20). Using
the exact sequence (*), one checks that $Ker$ is a direct sum of line bundles
of degree $-1$. Hence $Ker= S^{g-2}\bar H\otimes \wedge ^2\bar H\otimes
{\cal O}(-1)$. As $E^*=\pi^*(Ker), H=\pi ^*\bar H$, the result follows.

\end{pf}

\begin{Prop}
\label{W^r}
Let the notations be as above. Choose $r\le g-1$
Then,
 $$\wedge ^rE=\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes L^{\otimes r}$$
Moreover, the map
$$\wedge^rH^0(E) \rightarrow H^0(\wedge ^r E)$$
can be identified to
$$\varphi _r :  \wedge^r (S^{g-1}H^*)\rightarrow \wedge ^r(S^{g-2}H^*
\otimes \wedge ^2H^*)\otimes S^rH$$
and is an immersion.
 In particular, its image $W^r$ has dimension
 $\begin{pmatrix} g\\r\end{pmatrix}$

\end{Prop}

\begin{pf}
We can give a coordinate description of the map
$$\varphi_1: S^{g-1}H^*\rightarrow S^{g-2}H^*\otimes \wedge ^2H^*\otimes H$$
as follows. For any positive integer $k$, identify $S^kH$ with the
 space of polynomials of degree at most $k$ in a variable $x$
and $S^kH^*$ with the space of polynomials of degree at most $k$ in one
variable $y=x^*$. Identify $\wedge ^2H^*\rightarrow {\bf k}$
by the isomorphism $y\wedge 1\rightarrow 1$.
Then, the map  $\varphi_1$ acts by
$$\varphi_1(y^k)=y^k\otimes x-y^{k-1} \otimes 1$$
with the convention that on the right hand side $y^{g-1}=0, y^{-1}=0$.
By taking the wedge products of this map, one then sees that
$$\varphi_r (y^{k_1}\wedge ...\wedge y^{k_r})=\sum_{0\le \epsilon_i\le 1}
(-1)^{\epsilon _1+...+\epsilon _r}y^{k_1-\epsilon_1}\wedge...
\wedge y^{k_r-\epsilon_r}
\otimes x^{r-(\epsilon _1+...+\epsilon _r)}$$
 again with the convention that on the right hand side $y^{g-1}=0,
 y^{-1}=0$.

A basis of $\wedge^r (S^{g-1}H^*)$ consists of the elements
$y^{k_1}\wedge ...\wedge y^{k_r}$, $0\le k_1<...<k_r\le g-1$. The images
 of these elements are linearly independent as their leading terms (i.e.
the term of highest degrees on x and y jointly) obviously
are. Therefore, the dimension of $W^r$ is as stated.
\end{pf}

\section{General set up for infinitessimal deformations}
\begin{Rec}
\label{def}
We recall the basic set up for deformations of a curve, a vector bundle
and its space of sections.
\end{Rec}
Write ${\bf k}[t]/t^2={\bf  k}_{\epsilon }$. By an infinitessimal
 defformation of the curve $C$ we mean a curve ${\cal C}_{\epsilon}$
over $Spec{\bf  k}_{\epsilon }$ with central fiber $C$. Similarly,
by an infinitessimal deformation of a vector bundle $F$ we mean
a vector bundle over $C\times Spec{\bf  k}_{\epsilon }$ with central
fiber $F$. By
an infinitessimal deformation of the pair, we mean a curve
${\cal C}_{\epsilon}$ and a vector bundle ${\cal E}_{\epsilon}$
over ${\cal C}_{\epsilon}$.

Recall that the set of infinitessimal deformations of the curve $C$
can be parametrised by $H^1(C,T_C)$, the set of infinitessimal
deformations of the vector bundle $F$ can be parametrised by
$H^1(F^*\otimes F)$ while the set of infinitessimal deformations of a pair
consisting of a curve $C$ and a vector bundle $F$ on $C$ can be parametrised
by $H^1(\sum _F)$ where $\Sigma _F$ denotes the sheaf of first order
differential operators acting on $F$.

We describe next the correspondence between these objects (cf.\cite{W}
proof of Prop. 1.2 and also \cite{BR} proof of 2.3). Assume given an element
$\nu \in H^1(T_C)$. We think of the
 sections of the sheaf  $T_C$ over an open set $U$ as the set of
(${\bf k}$-linear)-maps ${\cal O}_U\rightarrow {\cal O}_U$ satisfying
$\nu (fg)=\nu (f)g+f\nu (g)$.
Take an affine open cover $C=\cup U_i$. Write $U_{ij}$ for $U_i\cap U_j$.
 Represent $\nu $ by  a
cocycle $\nu =(\nu _{ij}), \nu _{ij}\in H^0(U_{ij},T_C)$. We associate
to $\nu$ the following deformation of $C$: Consider the trivial
deformations of the $U_i$, namely $U_i\times Spec {\bf k}_{\epsilon }$.
Glue them along the intersections $U_{ij}\times Spec {\bf k}_{\epsilon }$
using the matrices
$$\begin{pmatrix} Id&0\\
 \nu_{ij}&Id
\end{pmatrix}.$$
The correspondence
$\nu _{ij}\rightarrow {\cal C}_{\epsilon}$
 obtained in this way is a bijection.

Assume now given an element $\varphi \in H^1(F^*\otimes F)$.
 Represent it by a cocycle $(\varphi _{ij})$ with
$\varphi _{ij} \in H^0(U_{ij},Hom(F,F))$. Consider the trivial extension
of $F$ to $U_i\times Spec {\bf k}_{\epsilon }$, namely
$F_{U_i}\oplus \epsilon F_{U_i}$. Take gluings on $U_{ij}$ given
by
$$\begin{pmatrix} Id&0\\
 \varphi _{ij}&Id
\end{pmatrix}.$$
This gives the correspondence between $H^1(F^*\otimes F)$
 and deformations of $F$.

 Assume now that a section $s$ of $F$
can be extended to a section $s_{\epsilon}$ of the deformation.
There exist then local sections $s'_i \in H^0(U_i, F_{|U_i})$
such that $(s_{|U_i},s_i)$ define a section of $F_{\epsilon}$.
By construction of $F_{\epsilon}$ this means that
$$\begin{pmatrix} Id&0\\
 \varphi_{ij}&Id
\end{pmatrix}
\begin{pmatrix} s_{|U_i}\\
s'_i
\end{pmatrix}=
\begin{pmatrix} s_{|U_j}\\
s'_j
\end{pmatrix}.$$

This can be written as $\varphi_{ij}(s)=s'_j-s'_i$.
Equivalently,
$$\begin{matrix} \varphi _{ij}\in Ker (& H^1(F^*\otimes F)&
\rightarrow &H^1(F))\\
 &\nu _{ij}&\rightarrow &\nu_{ij}(s) \end{matrix}$$
This result can be formulated using the language of Brill-Noether Theory:
the set of infinitessimal deformations of the vector bundle $F$ that
 have sections deforming  a certain subspace $V\subset H^0(F)$
 consists of the orthogonal to the image  of the Petri map
$$P_V: V\otimes H^0(K\otimes F^*)\rightarrow H^0(K\otimes F\otimes F^*).
\leqno (\ref{def}.1)$$

Assume now given an element $\sigma \in H^1(\Sigma  _F)$.
We think of  $\Sigma _F(U)$ as the set
 of additive morphisms $\sigma :F(U)\rightarrow F(U)$ such that for a
suitable element $\nu _{\sigma} \in T_C$,
$\sigma (fs)= \nu_{\sigma } (f)s+f\sigma (s)$.
Represent $\sigma$ by a cocycle $\sigma =(\sigma _{ij}),\  \sigma_{ij}\in
H^1(U_{ij},\Sigma _F)$. Consider the associated element $(\nu _{ij})
\in H^1(T_C)$ and the corresponding deformation ${\cal C}_{\epsilon}$
of $C$. Take then the vector bundle on ${\cal C}_{\epsilon}$ obtained
by gluing the trivial extensions of $F$ on $U_i$ by means of
the matrices
$$\begin{pmatrix} Id&0\\
 \sigma_{ij}&Id
\end{pmatrix}.$$

As in the case of deforming the line bundle alone, deformation of sections
is easy to interpret:
the set of infinitessimal deformations of the pair $(C,F)$ that have sections
deforming  a certain subspace $V\subset H^0(F)$
consists of the orthogonal to the image  of the Petri map

$$\bar P_V:V\otimes H^0(K\otimes F^*)\rightarrow H^0(K\otimes \Sigma _F^*)
\leqno (\ref{def}.2)$$

defined as the dual of the natural cup-product map

$$H^1(\Sigma _F)\rightarrow Hom(V,H^1(F)).$$

Consider the exact sequence
$$
0\rightarrow F^*\otimes F\rightarrow \Sigma_F\rightarrow T_C\rightarrow 0.$$
The canonical map $\pi :\Sigma _F \rightarrow T_C$ is defined
by $\pi (\sigma)=\nu _{\sigma}$. The map $i:F^*\otimes F\rightarrow \Sigma_F$
sends an element of $F^*\otimes F$ (considered as an endomorphism
of $F$) to itself.
One obtains a commutative diagram
$$\begin{array}{cccccccccc}
 &0&\rightarrow &H^1(F^*\otimes F)&\rightarrow &H^1(\Sigma _F)
&\rightarrow&H^1(T_C)&\rightarrow &0\cr
(\ref{def}.3)& & &\downarrow P_V^*& &\downarrow \bar P_V^* &
&\downarrow P_V^{'*} & & \cr
 &0&\rightarrow &ImP^*&\rightarrow &Hom(V,H^1(F))&\rightarrow&
(Ker P)^*&\rightarrow &0\cr
\end{array}$$
and its dual (cf. \cite{AC}p.18 )
$$\begin{array}{cccccccccc}
 &0&\leftarrow&H^0(K\otimes F\otimes F^*)&\leftarrow
&H^0(K\otimes (\Sigma_F)^*)&\leftarrow &H^0(2K)&\leftarrow &0\cr
(\ref{def}.4)& & &\uparrow & &\uparrow \bar P_V& &\uparrow P'_V& & \cr
 &0&\leftarrow &ImP& \leftarrow &V \otimes H^0(K\otimes F^*))&
\leftarrow & Ker P&\leftarrow &0\cr
\end{array}$$
When $V=H^0(F)$, we shall write $P_F$ instead of $P_V$. When $V$
and $F$ are clear, we shall suppress them from the notations.

We shall later use the following result. Its proof appears in \cite{T}
Lemma 2.12.

\begin{Lem}
\label{1x-x1}
Let $M$ be a line bundle on a curve $C$ with two independent sections $s_0,s_1$
and such that $|K\otimes M^{-2}|$ has a section $t$. Denote by $D_1$
 the fixed part of the series determined by $s_0,s_1$,
 denote by $D_2$ the divisor corresponding to the section
$t$. Denote by $R$
the ramification divisor of the map $C\rightarrow {\bf P}^1$
 associated to the series $<s_0,s_1>$.
Then, $P'_M(s_0\otimes ts_1-s_1\otimes ts_0)$ corresponds to the divisor
$2D_1+D_2+R$. In particular it is non-zero.
\end{Lem}

\begin{Lem}
\label{F=sum}
Assume that $F=\oplus _{i=1}^n F_i$ is a direct sum of vector bundles.
Then $$\Sigma _F=\oplus _{T_C}\Sigma_{F_i}\oplus [\oplus_{i\not= j}
F_i^*\otimes F_j]$$
Here $\oplus _{T_C}\Sigma_{F_i}$ denotes the fibered product over $T_C$
of the $\Sigma_{F_i}$.
\end{Lem}
\begin{pf}
 Consider an open set $U$. Let $\sigma :F(U)\rightarrow F(U)$ be
a first order differential operator acting
on $F$. Using the decomposition of $F$ as a direct sum, $\sigma$ admits
 a representation as a matrix $(\sigma _{ij})$ where
$\sigma _{ij}:F_i\rightarrow F_j$. Take a local section $s_k\in F_k(U)$.
One then checks that
 $$\sigma ((0...0,fs_k,0...0))=(\sigma_{1k}(fs_k),...,\sigma_{nk}(fs_k))$$
Using that $\sigma ((0...0,fs_k,0...0))=f\sigma ((0...0,s_k,0...0))
+\nu (f)(0...0,s_k,0...0)$, we find that $\sigma _{kk}$ is a first order
differential operator corresponding to the same $\nu$ as $\sigma$
while $\sigma _{ij}$ is ${\cal O}_C$-linear if $i\not= j$.
\end{pf}

\begin{Lem}
\label{F=sum'}
Assume that $F=\oplus _{i=1}^n F_i$ is a direct sum of vector bundles.
 Denote by $P_i,\bar P_i,P'_i$ the Petri maps corresponding to the
$F_i$. Then $P'_F$ can be obtained as the composition
$$Ker P\rightarrow \oplus _{i=1}^n KerP_i\rightarrow H^0(2K)$$
where the first map is the projection and the second map is ${1\over n}
\oplus P'_i$
 \end{Lem}
\begin{pf}
Consider the right hand square in \ref{def}.3 for each one of the $F_i$.
As in \ref{F=sum}, consider the fibered product of the $\Sigma_{F_i}$ over
$T_C$. One then has a  commutative square
$$\begin{array}{ccc}
\oplus_{H^1(T_C)} H^1(\Sigma _{F_i})&\rightarrow&H^1(T_C)\cr
\downarrow \oplus_{i=1}^n\bar P_i^*& &\downarrow \oplus_{i=1}^n  P_i^{*} \cr
\oplus_{i=1}^nHom(H^0(F_i),H^1(F_i))&\rightarrow&\oplus_{i=1}^n coker P_i^{'*}
\cr
\end{array}$$
 Take also the corresponding square for $F$
$$\begin{array}{ccc}
 H^1(\Sigma _F)&\rightarrow&H^1(T_C)\cr
\downarrow \bar P* & &\downarrow P^{'*} \cr
Hom(H^0(F),H^1(F))&\rightarrow&coker P^*\cr
\end{array}$$
Consider the cube that has these diagrams as back and front faces
respectively. We define four maps in the side edges.
The maps
$$\oplus_{H^1(T_C)} H^1(\Sigma _{F_i})\rightarrow H^1(\Sigma _F) $$
and
$$\oplus_{i=1}^nHom(H^0(F_i),H^1(F_i)) \rightarrow   Hom(H^0(F),H^1(F))$$
are natural diagonal injections (with zeroes on the terms corresponding
to a pair $F_i,F_j,i\not= j$). Notice that from \ref{F=sum}, the first map
is well defined. With these definitions, the left hand square commutes.

The map
$$\oplus_{i=1}^n (Ker P_i)^* \rightarrow (Ker P)*$$
is defined as the dual of the natural projection. By dualisation, one can then
check that the bottom face commutes.

If we take as the fourth map the homotethy
$$\times n :H^1(T_C) \rightarrow H^1(T_C)$$
 then, the top face commutes too.
 As $H^1(\Sigma _F)\rightarrow H^1(T_C)$ is onto,
 this shows that the right hand square commutes. Dualising this square,
one obtains the result in the Lemma.

\end{pf}

\section{Deformations of E}

We apply the set up of the previous section to the
 hyperelliptic curve $C_0$ and
the vector bundles $\wedge ^rE$.  As in section 1, $L$ denotes the
 hyperelliptic line bundle (i.e. the line bundle on $C$ of degree
two with two sections).
 We shall assume in all that follows
that $r\le g-1-r$

\begin{Prop}
\label{P}
 The Petri map $P_{H^0(\wedge ^r E)}$ (cf. (\ref{def}.2))
associated to the vector bundle
 $\wedge ^r (E)$ gives by restriction an isomorphism
$$P_{W^r}: W^r\otimes H^0(K\otimes (\wedge ^r E)^*)\rightarrow
 H^0(K\otimes \wedge^r E\otimes (\wedge ^r E)^*)$$
\end{Prop}
\begin{pf}
From \ref{W^r}
$$\wedge ^rE=\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes L^r$$
 so one has
 $$\wedge ^rE^*=\wedge ^r(S^{g-2}H\otimes \wedge ^2H)\otimes L^{-r}$$
Hence,
$$H^0(\wedge ^rE)=\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes
S^rH.$$
and
 $$H^0(K\otimes\wedge ^rE^*)=\wedge ^r(S^{g-2}H\otimes \wedge ^2H)\otimes
S^{g-1-r}H.$$
Then, the Petri map $P_{H^0(E)}$ and its restriction $P_{W^r}$
to $W^r\otimes H^0(K\otimes (\wedge ^r E)^*) $ can be written as the tensor
product with the vector space $ \wedge ^r(S^{g-2}H\otimes \wedge ^2H)$
of the diagram
$$\begin{matrix}
\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes S^rH\otimes S^{g-1-r}H&
\rightarrow &\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes S^{g-1}H\cr
 \uparrow & &\uparrow \cr
W^r\otimes S^{g-1-r}H& \rightarrow &\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)
\otimes S^{g-1}H\cr \end{matrix}$$
Therefore, it is enough to show that the map $p_{W^r}$ in the lower row of
this diagram is an isomorphism.
Recall that $\varphi _r$ identifies $\wedge ^rS^{g-1}H^*$ to its image $W^r$
(where $\varphi _r$ is the natural immersion
$\wedge ^rS^{g-1}H^*\rightarrow \wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)
\otimes S^r$ defined in \ref{W^r}). The natural product map $S^rH\otimes
S^{g-1-r}H\rightarrow S^{g-1}H$ can be identified to $P_{L^r}$.
Then, $p_{W^r}=(Id\otimes P_{L^r})o(\varphi_r\otimes Id)$.
We shall show that $p_{W^r}$ is an isomorphism by exhibiting its inverse
$q_{W^r}$.
Using the notations in \ref{W^r}, and the identification
$\wedge ^2H^*\cong {\bf k}$, one obtains
$$p_{W^r}(y^{k_1}\wedge ...\wedge y^{k_r}\otimes x^a)=
(Id\otimes P_{L^r})o(\varphi_r\otimes Id)
(y^{k_1}\wedge ...\wedge y^{k_r}\otimes x^a)=$$
$$=Id\otimes P_{L^r}
(\sum _{0\le \epsilon_i\le 1}(-1)^{\epsilon_1+...+\epsilon _r}
y^{k_1-\epsilon_1}\wedge ...\wedge
y^{k_r-\epsilon_r}\otimes x^{r-(\epsilon _1+...+\epsilon _r)}\otimes x^a)=$$
$$=\sum _{0\le \epsilon_i\le 1}(-1)^{\epsilon_1+...+\epsilon _r}
y^{k_1-\epsilon_1}\wedge ...\wedge
y^{k_r-\epsilon_r}\otimes x^{a+r-(\epsilon _1+...+\epsilon _r)}$$
with the convention that on the right hand side $y^{-1}=0,\ y^{g-1}=0$.
We describe $q_{W^r}$ as follows. Assume given integers  $0\le j_1<...<j_r
\le g-2,0\le b\le g-1$. There is then a value $l$ with $0\le l\le r$ such
that $j_l+1\le b\le j_{l+1}$. Define then
$$q_{W^r}(y^{j_1}\wedge ...\wedge y^{j_r}\otimes x^b)=
\sum _{0\le t_i\le j_i-j_{i-1}-1,1\le s_i\le j_{i+1}-j_i}(-1)^{r-l}
y^{j_1-t_1}\wedge ...\wedge y^{j_l-t_l}\wedge y^{j_{l+1}+s_{l+1}}\wedge ...
y^{j_r+s_r}$$
$$\otimes x^{b-r-t_1-...-t_l+s_{l+1}+...+s_r}$$
with the conventions $j_{r+1}=g-1, j_{-1}=-1$.
Notice that the map $q_{W^r}$ is well defined as
$0\le b-r-t_1-...-t_l+s_{l+1}+...+s_r\le g-1-r$ and
$0\le j_1-t_1<...<j_l-t_l<j_{l+1}+s_{l+1}\le ...\le j_r+s_r\le g-1$.
It is a slightly tedious but straightforward computation to show that
the composition
$p_{W^r}oq_{W^r}=Id$. As the two vector spaces involved have the same
 dimension, this suffices in order to prove the isomorphism.
\end{pf}

\begin{Cor}
\label{sigma_nu}
For any given infinitessimal deformation $\nu$ of the curve $C_0$,
there is a unique infinitessimal deformation $\sigma$ of the
pair consisting of $C_0$ and the vector bundle $\wedge ^rE$
such that $W^r$ can be extended to a space of sections of the deformation.
\end{Cor}
\begin{pf}

Consider the Petri map
 $$P_{W^r}:
W^r\otimes H^0(K\otimes F^*)\rightarrow H^0(K\otimes F\otimes F^*)$$
Consider the commutative diagram (\ref{def}.3) in case $F=\wedge ^rE, V=W^r$
$$\begin{array}{ccccccccc}
0&\rightarrow &H^1((\wedge ^rE)^*\otimes \wedge ^rE)&\rightarrow &H^1(\Sigma _
{\wedge ^rE})
&\rightarrow&H^1(T_C)&\rightarrow &0\cr
 & &\downarrow P^*& &\downarrow \bar P^* & &\downarrow P^{*'} & & \cr
0&\rightarrow &ImP^*&\rightarrow &Hom(W^r,H^1(\wedge ^r E))&\rightarrow&
(ker P)^*&\rightarrow &0\cr
\end{array}$$
From \ref{P},  $P^*$ is an isomorphism, $KerP=0$. Hence, every element in
 $H^1(T_C)$ has a
unique inverse image in $H^1(\Sigma _F)$ that belongs to the kernel
of $\bar P^*$. This proves the result.
\end {pf}

\begin{Cor}
\label{sigma_nu'}
For any given infinitessimal deformation $\nu$ of the curve $C_0$, the
unique infinitessimal deformation $\sigma$ of the vector bundle
$\wedge ^rE$ that preserves it as the  exterior power of the dual of
the kernel of the evaluation map is the $\sigma $ above.
\end{Cor}
\begin{pf} An infinitessimal deformation of $E$ as the dual of the kernel
of the evaluation map preserves the subspace $W=W^1$ as space of
sections . Hence, an infinitessimal deformation of $\wedge ^rE$ as the
exterior product of this dual preserves $W^r$ as space of sections.
By the unicity of such deformation, the result follows.
\end{pf}

\begin{Prop}
\label{noW'}
Take $\nu$ a generic infinitessimal deformation of $C_0$. Let $\sigma$
be the deformation of $\wedge ^r E$ associated to $\nu$ as in \ref{sigma_nu'}.
If $\hat W$ is a subspace of $H^0(\wedge ^r(E))$ that strictly contains
$W^r$, then $\hat W$ does not extend to a space of sections of
the infinitessimal deformation
of $\wedge ^r E$ corresponding to $\sigma $.
\end{Prop}
\begin{pf}
It is enough to prove the result when $\hat W$ has dimension $a=dimW^r+1$.
Denote by $S$ the Schubert cycle
of subspaces of dimension $a$ of $H^0(\wedge ^rE)$ that contain
$W^r$ ($S\subset {\bf Gr}(a, H^0(\wedge ^rE))$)  .

Consider the diagram (\ref{def}.4) for the case $F=\wedge ^rE, V=\hat W$.
 We obtain
$$\begin{array}{ccccccccc}
0&\leftarrow&H^0(K\otimes (\wedge ^rE)\otimes (\wedge ^rE)^*)&\leftarrow
&H^0(K\otimes (\Sigma_{\wedge^rE})^*)&\leftarrow &H^0(2K)&\leftarrow &0\cr
 & &\uparrow & &\uparrow \bar P_{\hat W}& &\uparrow P'_{\hat W}& & \cr
0&\leftarrow &ImP& \leftarrow &\hat W\otimes H^0(K\otimes (\wedge ^rE)^*)&
\leftarrow & Ker P_{\hat W}&\leftarrow &0\cr
\end{array}$$
Define $P'=P'_{H^0(\wedge ^r(E))}$.

We first show that the result would follow from the following statement:

(*) There is an element $\alpha \in H^0(2K)$ such that $\alpha \in P'(\hat W)$
for all $\hat W in S$.

If (*) holds, consider a direction of deformation $\nu$
in $H^1(T_C)$ not orthogonal to $\alpha$ and its corresponding
deformation $\sigma_{\nu}$. Then, $\bar P^*_{\hat W}(\sigma _{\nu })\not= 0$
 for
all $\hat W$. Therefore, no $\hat W$  extends to a space of sections
of the deformation corresponding to $\sigma _{\nu}$.

Let us see that in fact it suffices to prove (*) for the generic
 $\hat W \in S$.
Note then that from \ref{P}, $Ker P_{\hat W}$ is a subspace of
 $H^0(\wedge ^rE)\otimes H^0(K\otimes (\wedge ^rE)^*)$ of fixed dimension
 $b=h^0(K\otimes (\wedge ^rE)^*)$.  We have a well-defined map

$$\begin{array}{cccc} \varphi : &S&\rightarrow &
{\bf Gr}(b,H^0(\wedge ^rE)\otimes H^0(K\otimes (\wedge ^rE)^*)\cr
 &\hat W&\rightarrow &KerP_{\hat W}
\end{array}$$

Define $\bar P=\bar P_{H^0(\wedge ^r(E))}$. Assume that for the generic
$\hat W\in S$, $\varphi (\hat W)$ intersects
$(P')^{-1}(\alpha)$ for a fixed element $\alpha \in H^0(2K)$.
Denote by $S'$ the Schubert cycle of $ {\bf Gr}(b,H^0(\wedge ^rE)
\otimes H^0(K\otimes (\wedge ^rE)^*)$ of subspaces that intersect
$(P')^{-1}(\alpha)$. Our assumption is that for $\hat W\in S$
generic, $\varphi (\hat W)\in S'$. From the irreducibility of $S$
and  the closedness of $S'$, the same is true for every element in
$S$.

 Hence, it remains to prove (*) for generic $\hat W\in S$.
This  will follow from the Lemma below and \ref{1x-x1}.
\end{pf}
\begin{Lem}
\label{1x-x1'}
Let $\hat W$ be generic in $S$, then
 $P'_{\hat W}(Ker P_{\hat W})=P'_{L^r}(Ker P_{L^r})  $.
\end{Lem}

\begin{pf}[of \ref{1x-x1'}]
As noticed before, $dimKer P_{\hat W}$ is constant for $\hat W\in S$.
 Therefore, the map
$P'_{\hat W}$ attains its maximum rank on an open set of $S$.
From \ref{F=sum'}, the image of $P'_{\hat W}$ is contained in the
image of $P'_{L^r}$. Hence, the maximum possible rank of $P'_{\hat W}$
is $dim ImP'_{L^r}$. It then suffices
to exhibit one $\hat W$ satisfying the condition.
 We use the basis of $\wedge ^r (S^{g-2}H^*\otimes \wedge ^2H^*)$ introduced
 in \ref{W^r} and the corresponding trivialisation of $\wedge ^2H^*$.
Define
$$\hat w=\sum _{j=0}^{r-1}(-1)^j1\wedge y\wedge ...\wedge y^{j-1}\wedge y^{j+1}
\wedge ...\wedge y^r\otimes x^{r-j}$$
Let $\hat W$ be the space generated by $\hat w$ and $W^r$.
We recall from the proof of \ref{P} that the map $P_{\hat W}$ can be
 obtained by taking tensor product with
$ \wedge ^r(S^{g-2}H\otimes \wedge ^2H)$
of the map
$$p_{\hat W}:\hat W\otimes S^{g-1-r}H\rightarrow
\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes S^{g-1}H$$
 Consider the
following element of $\hat W\otimes S^{g-1-r}H$
$$z=\hat w\otimes 1+\sum _{i=r+1}^{g-1}\varphi_r(y\wedge y^2\wedge ...
\wedge y^{r-1}\wedge y^i)\otimes x^{i-r}$$
We want to see that $p_{\hat W}(z)=0$.

We can write an explicit expression for $z$ using the convention
that $y^{-1}=0,y^{g-1}=0$
$$z=\sum _{j=0}^{r-1}(-1)^j1\wedge y\wedge ...\wedge y^{j-1}\wedge y^{j+1}
\wedge ...\wedge y^r\otimes x^{r-j}\otimes 1$$
$$+\sum_{i=r+1}^{g-1}\sum_{0\le \epsilon_i\le 1}
(-1)^{\epsilon _1+...+\epsilon _r}y^{1-\epsilon_1}\wedge...
\wedge y^{r-1-\epsilon_{r-1}}\wedge y^{i-\epsilon _r}\otimes
 x^{r-(\epsilon _1+...+\epsilon _r)}\otimes x^{i-r}$$
$$=\sum _{j=0}^{r-1}(-1)^j1\wedge y\wedge ...\wedge y^{j-1}\wedge y^{j+1}
\wedge ...\wedge y^r\otimes [x^{r-j}\otimes 1-x^{r-(j+1)}\otimes x]$$
$$+\sum_{i= r+1}^{g-2}\sum_{0\le \epsilon_i\le 1}
(-1)^{\epsilon _1+...+\epsilon _{r-1}}y^{1-\epsilon_1}\wedge...
\wedge y^{r-1-\epsilon_{r-1}}\wedge y^i\otimes$$
$$[ x^{r-(\epsilon _1+...+\epsilon _{r-1})}\otimes x^{i-r}-
 x^{r-(\epsilon _1+...+\epsilon _{r-1}+1)}\otimes x^{i+1-r}]$$
It is clear from the above expression for $z$ that $z\in Kerp_{W'}$.
In particular, for any element $e\in \wedge ^r(S^{g-2}H\otimes \wedge ^2H)$,
$z\otimes e \in KerP_{W'}$.

From \ref{F=sum'}, the map $P'_{\hat W}$ is up to a constant $c'$
($c'={1\over { g-2\choose r}}$) the composition of the
following maps:
natural inclusion
$$\hat W\otimes \wedge ^r(S^{g-2}H\otimes \wedge ^2H)\otimes S^{g-1-r}H
  \rightarrow $$
$$\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes
\wedge ^r(S^{g-2}H\otimes \wedge ^2H)\otimes S^{g-1}H    $$
 the cup-product with the identity element
in $\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*) \otimes
\wedge ^r(S^{g-2}H\otimes \wedge ^2H)$
$$\wedge ^r(S^{g-2}H^*\otimes \wedge ^2H^*)\otimes
\wedge ^r(S^{g-2}H\otimes \wedge ^2H \otimes S^{g-1}H
\rightarrow S^rH\otimes S^{g-1-r}H$$
followed by the natural product map
$$P'_{L^r}:S^rH\otimes S^{g-1-r}H\rightarrow S^{g-1}H$$.

 Hence, for $j=0,...,r-1$ writting $c=(-1)^jc'$,
$$P'_{\hat W}(z\otimes 1\wedge x\wedge ...\wedge x^{j-1}\wedge x^{j+1}
\wedge ...\wedge x^r)=cP'_{L^r}[x^{r-j}\otimes 1-x^{r-(j+1)}\otimes x]$$
and for $i=r+1...g-2, 0\le \epsilon _i\le 1$ writting
 $c=(-1)^{\epsilon_1+...+\epsilon_{r-1}}c'$
$$P'_{\hat W}(z\otimes x^{1-\epsilon_1}\wedge...\wedge x^{r-1-\epsilon_{r-1}}
\wedge x^i )=
cP'_{L^r}[x^{r-(\epsilon _1+...+\epsilon _{r-1})}\otimes x^{i-r}-
 x^{r-(\epsilon _1+...+\epsilon _{r-1}+1)}\otimes x^{i+1-r}]$$
As $ker P_{L^r}$ is generated by elements of the form
$$x^i\otimes x^j-x^{i+1}
\otimes x^{j-1}, \  i=0...r-1,\ j=1...g-1-r$$
the result is proved.
\end{pf}

\begin{section}{Extending the results to the generic curve.}

\begin{Prop}
\label{inj}
Let $C$ be a generic curve of genus $g$. Let $E$ be defined as in \ref{E}.
Denote by $W$ the image of $(H^0(C,K_C))^*$ in $H^0(C,E)$. Then,
 the natural map
$$\psi _{C,r}: \wedge ^rW \rightarrow H^0(C,\wedge ^rE)$$
is injective.
\end{Prop}
Note that for $C$ non-hyperelliptic, $W=H^0(E)$.
This follows from the projective normality of $C$ (case $p=0$ of the
conjecture).

\begin{pf}
If for a given curve $C$, $\psi _{C,r}$ is injective, the same holds
for every curve in a neighborhood of $C$ in ${\cal M}_g$.
As $dim \wedge ^rW=\binom{g}{r}$, \ref{W^r} shows that $\psi _{C_0,r}$
is injective for $C_0$ hyperelliptic. Hence the result follows.
\end{pf}

The following proposition now concludes the proof of \ref{teorema}:

\begin{Prop}
Let $C$ be a generic curve of genus $g$. Then, $h^0(C,\wedge ^rE)=\binom{g}
{r}$ and $\wedge ^rW\rightarrow H^0(C,\wedge ^rE)$ is an isomorphism.
\end{Prop}
\begin{pf}
From \ref{inj}, $h^0(C,\wedge ^rE)\ge dim Im
(\wedge ^rW \rightarrow H^0(C,\wedge ^rE))=dim\wedge ^rW=\binom{g}{r}$.
From \ref{noW'}, $h^0(C,\wedge^rE)\le dim W^r_{C_0}=\binom{g}{r}$.
This concludes the proof.
\end{pf}
\end{section}


text:
\begin{document}
\title{Numerical integration of discontinuous functions in many dimensions}

\author{Vassilis Pandis}

\maketitle
\begin{abstract}
We consider the problem of numerically integrating functions with hyperplane discontinuities
over the entire Euclidean space in many dimensions. We describe a simple process through which the
Euclidean space is partitioned into simplices on which the integrand is smooth, generalising
the standard practice of dividing the interval used in one-dimensional problems. Our procedure is
combined with existing adaptive cubature algorithms to significantly reduce the necessary number
of function evaluations and memory requirements of the integrator. The method is embarrassingly
parallel and can be trivially scaled across many cores with virtually no overhead.
Our method is particularly pertinent to the integration of Green's functions, a problem directly related
to the perturbation theory of impurity models. In three spatial dimensions we observe
a speed-up of order $100$ which increases with increasing dimensionality.
\end{abstract}

\section{Introduction}\label{sec:intro}

The numerical integration of functions in many dimensions has been a central topic in numerical
analysis for a long time. Current schemes such as adaptive cubature and adaptive Monte Carlo
perform best for smooth integrand functions. However integrands with discontinuities can arise
quite naturally in a variety of contexts, such as the calculation of the fermionic self-energy
in condensed-matter physics or the study of multiphase flows in the context of computational fluid dynamics.

In one dimension discontinuities are easily accommodated within an adaptive framework simply
by dividing the region of integration into sub-regions on which the integrand is smooth.
In this paper we show how this process can be extended to higher dimensionalities. Our method is applicable
to integrals which are discontinuous on any number of hyperplanes that contain the origin, and in any number of
dimensions. We limit our attention to integrals over the entire $\RN$ --- this
is not a material limitation as integrals over a proper hyperrectangle can be straightforwardly mapped onto
$\RN$. We assume that the discontinuities in question arise from terms of the form sign$(C_{\vec{x}})$ where
$C_{\vec{x}}$ is any linear combination of the coordinates. These is precisely the form of the discontinuities
encountered in the Green's functions of fermionic systems.

We write the integral in question as
\begin{equation}
I = \int_\RN   \prod_{i=1}^{M} F_i(\vec{x})d\vec{x},
\label{eq:integral}
\end{equation}
where $F_i(\vec{x})$ is discontinuous on the hyperplane with equation $\vec{a_i}\cdot \vec{x} = 0$.
We construct the $M\times N$ \emph{discontinuity} matrix $\mat{C}$ such that $C_{ij} = (\vec{a_i})_j$.
We will here assume $M \geq N$ --- we will comment on this at end the next section.

\section{Method}

Let $S$ be the set of all $M\times M$ diagonal matrices with diagonal components $\pm 1$ ($|S| = 2^M$).
To determine the regions on which the integrand is continuous we thus have to solve the homogeneous system of simultaneous
inequalities
\begin{equation}
\mat{C}\mat{S}_i\vec{x} \geq 0,
\label{eq:system}
\end{equation}
for every $\mat{S}_i \in S$. Each inequality defines a closed-half space; the solution to the system is
the intersection of these half-spaces which can be interpreted geometrically as a convex polytope in
its half-space representation  (see \cite[p. 31]{Grunbaum67}). In the case of a homogeneous system the resultant polytope is
in fact a polyhedral (infinite) convex cone~\cite{skeleton}.

Let $P$ be the set of cones obtained by solving Eq.~\eqref{eq:system}.
A set of vectors $W_K = \left\{\vec{w}_1, \vec{w}_2, \ldots ..\vec{w}_p \right\}$ is a \emph{skeleton} of a cone $K$
if $\vec{x} = \sum_{i=1}^p \lambda_i \vec{w_i}$ belongs to $K$ for every $\lambda_i \geq 0$~\cite{skeleton}.
The duality in the representation of a cone as either a system of linear inequalities or a conical combination
of the skeleton is the essence of the well known Weyl-Minkowski theorem on cones. As we are only interested in
subspaces of $\RN$ with dimension $N$ --- lower-dimensional subspaces correspond to polyhedral facets which do
not contribute to the integral --- we can assume that $p \geq N$.
The skeleton of an acute cone is unique up to scalar multiplication of the vectors~\cite{skeleton}.
Once the normalization of the skeleton is fixed, each point in $K$ can be specified through its $\lambda$ coefficients $(\lambda_1, \ldots, \lambda_p)$.

Having achieved our goal of partitioning $\RN$ into regions where the integrand is
continuous we now have to consider how to perform the integration over a cone $K\in P$. When $p=N$ the polytope
constitutes an $N$-simplex which can be readily mapped onto the positive orthant by exploiting the bijection between
$\vec{x} \in K$ and $\vec{\lambda}$. When $p>N$ the situation is more complex, for the skeleton is linearly dependent
and there is no bijection to be exploited. To overcome this problem, each cone $K$ is decomposed into $N$-simplices
$\gamma^{K}_1, \gamma^{K}_2., \ldots$ which can then be individually mapped onto the positive orthant. The
set of all simplices $\digamma = \left\{ \gamma^K_1, \gamma^K_2, \ldots |  K \in P \right\}$ evidently partitions $\RN$; the original
integration problem has thus been broken down into multiple, separate integrations, one over each
simplex in $\digamma$. The method is inherently parallel --- barring error control considerations each region
of integration can be processed independently of the others.

To control the precision of the calculation we use an unsophisticated two-pass scheme. The first pass consists of a
crude integration over every simplex $\gamma \in \digamma$, with a relative precision of $10\%$, yielding a
result $\mu^{(1)}_\gamma$ with an associated error $\sigma^{(1)}_\gamma$. From the $\mu^{(1)}_\gamma$ we determine the
simplex which contributes the most; let $\mu_{\textrm{max}} = \textrm{max} \left\{|\mu_\gamma|, \gamma \in \digamma \right\}$.
To achieve a requested relative precision $f$ on the entire integral $I$ we then repeat the integration, now evaluating
each simplex to an \emph{absolute} precision given by $\epsilon_{\textrm{abs}} = f\mu_\textrm{max}/\sqrt\nu$, where $\nu = |\digamma|$
denotes the total number of simplices, ensuring obviously we do not re-evaluate the regions for which $\sigma^{(1)}_\gamma < \epsilon_\textrm{abs}$.
The end result $I=\sum \mu^{(2)}_\gamma$ is then associated with an absolute error
\begin{equation}
\sigma_I = \sqrt{\sum_{\gamma \in \digamma} (\sigma^{(2)}_\gamma)^2/\nu}.
\end{equation}
In practice small deviations of the resultant precision for the requested precision may occur when there are
significant cancellations. This is not a particularly grave disadvantage as the actual error is always known.

Finally, we return to the question of the number of constraints. We have been assuming that the number of rows $M$
of the constraint matrix $\mat{C}$ is larger than the dimension of the integral, ignoring the case of an integrand which
has discontinuities on fewer than $N$ planes. This is dealt with by padding the rows of $\mat{C}$ with arbitrary vectors
(so long as they are not parallel to any other vectors) until $M\geq N$. This trick has the disadvantage of causing
unnecessary divisions of the region of integration but is necessary to guarantee the existence of cones.

\section{Implementation}

The process outlined above is implemented in \verb|C++| with support for matrices provided by GSL~\cite{gsl}. The
input is the integrand and the matrix $\mat{C}$ construct as above specifying the discontinuities. The first step is
the decomposition of $\RN$ into the polyhedral cones $P$. To this end we use \verb|skeleton|~\cite{skeleton}
which implements a modified version of the Motzkin-Burger algorithm. This package is called from our
code and returns the vectors comprising the skeleton of the polyhedral cones in $P$.

To cut the polyhedral cones in $K$ into $N$-simplices we first project the vectors in $W_K$ onto the
cone's $(N-1)$-dimensional base. By `base' here we mean the subspace obtained by subtracting from
all $\vec{w} \in W_k$ their components along the axis of the cone and then expressing them as linear
combinations of $(N-1)$ orthonormal vectors. We can then construct the desired decomposition of
$K$ into $N$-simplices by triangulating the points in the $(N-1)$-dimensional base and then adding
the origin to these $(N-1)$-simplices. In general this triangulation is not unique. There are
several algorithms to handle the triangulation of the base. We use the Quickhull algorithm
implemented in Qhull~\cite{qhull}.

Each point $\vec{x}$ of the simplex can be written as a conical combination of the $(\vec{\lambda})$
and its skeleton vectors. To map the positive orthant onto the unit hypercube we use the rule
$\lambda_i = 1/u_i - 1$. Depending on the integrand other rules may be more suitable but this
was chosen for its simplicity.

The final step is the integration itself. We use \verb|HIntLib|~\cite{hintlib1,hintlib2}, a sophisticated
\verb|C++| library that among other things implements adaptive cubature with a variety of rules and a range of Monte
Carlo methods. It would be perhaps more efficient to use an adaptive code that can directly handle the
simplicial geometry, such as \verb|CUBPACK|~\cite{cubpack1,cubpack2} but for practical reasons this
approach was not followed here. The integrations are performed in parallel using $\verb|OpenMP|$ (\verb|HIntLib|'s
native parallelization is not used).

\section{Results}

We test the method with a variety of integrands and for various dimensionalities. To do so we also have
to prescribe the discontinuities. To streamline the discussion we express Eq.~\eqref{eq:integral} as
\begin{equation}
I = \int_\RN   \prod_{i=1}^{M} F(g_i(\vec{x}))d\vec{x},
\label{eq:integral2}
\end{equation}
where $\vec{g}(\vec{x}) = \mat{C}\vec{x}$. A variety of test-matrices $\mat{C}$ are considered --- they are listed in the Appendix.
All integrations are done using \verb|HIntLib|'s adaptive routines and its implementation of the embedded degree-$7$ rule of Genz and Malik~\cite{Genz80}.

We define the following integrand test-functions
\begin{align}
F_1 (u) &= \frac{1}{u - \alpha + i \beta \textrm{sign}(u)}   \label{eq:f1} \\
F_2 (u) &= \frac{1}{u^2 - \alpha + i \beta \textrm{sign}(u)} \label{eq:f2}.

\end{align}
We note that $F_1$ is actually the non-interacting Green's function for the Anderson impurity model
in the flat-band approximation~\cite{Hewson97}. Our method was developed with this integrand
in mind --- we also consider the integrand $F_{2}$ to illustrate the more general applicability of the
method. As $F_{1}, F_{2}$ are complex-valued, we consider for brevity only the real part of Eq.~\eqref{eq:integral2}.

We compare the speed-up afforded by our partitioning in terms of the number of integrand evaluations
required to achieve a given relative error. In doing so, we seemingly ignore the computational effort
required for the partitioning itself. In practice this turns out to be essentially insignificant, owing
to the large computational cost of the integrations. Nevertheless the time spent in partitioning
can be reduced by noting that the solution of Eq.~\eqref{eq:system} and subsequent triangulation can be
carried out in parallel for each $\mat{S_i} \in S$.

Unless otherwise stated all integrations are carried out to $\epsilon_{\textrm{rel}} \approx 10^{-4}$. Due to the
two-pass technique for controlling the precision of the integration it may happen that the resultant error
is (sometimes significantly) less than requested. This is to be expected when the $\mu_\gamma$ have mostly the
same sign, i.e.\ $|I| \gg \mu_{\textrm{max}}$; in such cases the target absolute error --- which is based on $\mu_{\textrm{max}}$ ---
is smaller than necessary. The resultant error may be less than the requested in
another way: To obtain an estimate over each simplex $\gamma$, \verb|HIntLib| requires a minimum number $N_\textrm{min}$
of integrand evaluations. When this yields an estimate of the integral more precise than requested, nothing
can be done to reduce the number of evaluations.

\begin{table}
\caption{Results for $F_1$, $\alpha = -0.2,\beta=0.1$. \label{table:f1}}{
\begin{tabular}{l | c | c | c | r}
$N$   &   $M$    &   $N_{p}$   &   $N_{H}$   &  $N_{p}/N$ \\
\hline
2     &   3      & $5.3 \times 10^4$   &  $2.2 \times 10^6$ &  $41.4$ \\
2     &   4      & $8.6 \times 10^4$   &  $8.3 \times 10^5$ &  $9.6$  \\
2     &   5      & $1.3 \times 10^5$   &  $1.2 \times 10^6$ &  $9.0$  \\
2     &   6      & $1.7 \times 10^5$   &  $1.0 \times 10^6$ &  $6.3$  \\
3     &   5      & $3.3 \times 10^6$   &  $>3.0\times 10^9$ &  $>906.7$ \\
3     &   6      & $7.2 \times 10^6$   &  $1.6 \times 10^9$ &  $224.9$ \\
3     &   7      & $9.6 \times 10^6$   &  $2.6 \times 10^9$ &  $265.5$ \\
3     &   8      & $1.4 \times 10^7$   &  $2.5 \times 10^9$ &  $172.9$ \\
3     &   9      & $2.0 \times 10^7$   &  $2.7 \times 10^9$ &  $132.6$ \\
4     &   7      & $5.8 \times 10^8$   &  $>3.0\times 10^9$ &  $>5.2$ \\
5     &   9      & $4.3 \times 10^{10}$   & - &  -
\end{tabular}}
\end{table}
\begin{table}
\caption{Results for $F_2$, $\alpha = -0.2, \beta = 0.1$.\label{table:f2}}{
\begin{tabular}{l | c | c | c | c | c | c|  r}
$N$   &   $M$    &   $N_{p}$           & $\epsilon^{*}_\textrm{rel}$    &   $N_{H}$   & $N^{*}_{H}$  &  $N_{p}/N$ & $N^{*}_{p}/N$   \\
\hline
2     &    3     &   $3.6\times 10^4$   & $1.8\times10^{-7}$   & $1.3  \times 10^4$        & $2.5 \times 10^5$  &  $0.35$      & $7.0$ \\
2     &    4     &   $4.8\times 10^4$   & $4.6\times10^{-6}$   & $4.0  \times 10^5$        & $8.1 \times 10^6$  &  $8.3$       & $169$ \\
2     &    5     &   $6.0\times 10^4$   & $4.6\times10^{-6}$   & $6.5  \times 10^5$        & $1.4 \times 10^7$  &  $10.8$      & $236$ \\
2     &    6     &   $7.2\times 10^4$   & $1.6\times10^{-6}$   & $1.8  \times 10^6$        & $1.2 \times 10^8$  &  $24.9$      & $161$ \\
3     &    5     &   $1.7\times 10^5$   & $8.8\times10^{-5}$   & $9.8  \times 10^8$        & $1.2 \times 10^9$  &  $5820$      & $7418$   \\
3     &    6     &   $2.5\times 10^5$   & $7.7\times10^{-5}$   & $1.7  \times 10^8$        & $2.7 \times 10^9$  &  $6750$      & $10660$  \\
3     &    7     &   $3.8\times 10^5$   & $6.1\times10^{-5}$   & $>3.0 \times 10^9$        & $>3.0 \times 10^9$ &  $7804$      & $>7804$ \\
3     &    8     &   $4.8\times 10^5$   & $8.3\times10^{-5}$   & $>3.0 \times 10^9$        & $>3.0 \times 10^9$ &  $6244$      & $>6244$ \\
3     &    9     &   $5.8\times 10^5$   & $7.0\times10^{-5}$   & $>3.0 \times 10^9$        & $>3.0 \times 10^9$ &  $>5172$     & $>5172$ \\
4     &    7     &   $2.7\times 10^6$   & $9.8\times10^{-5}$   & $>3.0 \times 10^9$        & $>3.0 \times 10^9$ &  $>1111$     & $>5172$ \\
5     &    9     &   $4.8\times 10^{7}$& $9.9\times10^{-5}$   & $-$                       &  $-$               &   $-$          & $-$
\end{tabular}}
\end{table}

Our results are presented in Tables~\ref{table:f1}, ~\ref{table:f2}. In both tables $N$ denotes the dimension of integration, $M$ the number of hyperplane discontinuities,
$N_p$ the number of function evaluations required using the partitioning technique and $N_H$ the number of function
evaluations required to achieve the requested precision without utilising our partitioning scheme. To make the comparison
meaningful we obtain $N_{H}$ using the same adaptive cubature routines in \verb|HIntLib| with the embedded degree-7 Genz-Malik rule that were employed to carry out the simplicial integrations,
with each point $\vec{x}  \in \RN$ being mapped to $\vec{t} \in [-1, 1]^N$ through $x_i = t_i/(1-t^2_i)$. We emphasise
that our proposed integration method can be used in conjunction with any integration algorithm and is not tied to this
specific Genz-Malik rule.

We note that for many of the integrations in Table~\ref{table:f2}
we were unable to reduce the precision below a certain level. Thus for each integrand we report the relative precision $\epsilon^{*}_\textrm{rel}$ \emph{actually} reached
by our partitioning scheme. We feel it is not clear whether it would be  fairer to judge the efficacy of our method method
by comparing $N_p$ to the function evaluations required to achieve the target relative precision of $10^{-4}$ or the `accidental'
precision $\epsilon^{*}_\textrm{rel}$. We thus report both quantities, the latter denoted by $N^{*}_{H}$. From our
results it's obvious that even when calculating the integral to a precision much greater than required our method
greatly reduces the samples required of the integrand.

For the integrations attempted without partitioning $\RN$ we had to impose a maximum of $3\times10^9$ integrand evaluations to prevent the integrator
from exhausting the $8$~GB of RAM we had at our disposal. Apart from requiring fewer integrand samples, our partitioning
method also drastically reduces the amount of memory required. This is because each the grid for each simplex can be discarded
after the integration is complete rather than having to concurrently store data for all previous grid refinements. We have
however refrained from trying to quantify the improvement in the memory requirements as this is sensitive to the
details of our implementation, our choice of integration routines and the benchmarking itself rather non-trivial, given the parallel nature
of the program. Nevertheless in Table~\ref{table:f2} it is evident that the integration becomes unmanageable without
our method even in only four dimensions.

As the number of simplices into which $\RN$ is partitioned increases very rapidly with $N$, the success of the method depends on whether the advantages
of a smooth integrand outweigh the cost of having to set up a new adaptive grid for each simplex, and the
`unnecessary' function evaluations due to the crudeness of our error management. It is is evident that it does;
in Table~\ref{table:f1} we see that partitioning reduces the number of required integrand evaluations
by $1-3$ orders of magnitude and in Table~\ref{table:f2} by up to $4$ orders of magnitude.

\section{Conclusions}

We have described a method enabling the numerical integration of functions featuring hyperplane discontinuities
with existing adaptive cubature schemes. We showed how to construct a set $P$ of convex polyhedral cones
that partition $\RN$. Each cone $K \in P$ is then partitioned into simplices $\gamma^{K}_1,\gamma^K_2, \ldots$
which comprise the set $\digamma$ which partitions $\RN$. Each simplex $\gamma \in \digamma$ can then
be mapped onto a hypercube and integrated using any existing multidimensional numerical integration algorithm.
We adopted a two-pass scheme to control the precision of our calculation. This allowed the essentially
complete paralellization of the integrations.

Our method can dramatically accelerate the evaluation of such multidimensional integrals. The reduction
in the number of integrand samples required to obtain an estimate for the integral becomes more pronounced
as the dimensionality $N$ increases, and is reduced by several orders of magnitude compared to a naive integration.
Memory requirements are also greatly improved, allowing the evaluation of integrands in higher dimensions
than would be otherwise possible.

The method can be improved by coupling it directly to an integrator aware of the underlying simplicial
geometry, thereby eliminating the need for a simplex-hypercube mapping. Our precision control can also be potentially replaced with a more
advanced scheme in which no integrand evaluations are discarded and the threads dynamically synchronised.
This could improve the performance of our method but coordinating the threads will require some effort
programming-wise.

The author would like to thank Alex C Hewson, Paul Carter and Ioannis Pesmazoglou for helpful discussions. This work was
generously supported by the Engineering and Physical Sciences Research Council.

\section*{APPENDIX}
We list here the test matrices $\mat{C}_{c \times N}$ pertinent to Eq.~\eqref{eq:integral2}.
\begin{equation*}
C_{3\times2}=\begin{pmatrix}
1 & 0  \\
0 & 1  \\
1 & 1
\end{pmatrix} \,
C_{4\times2}=\begin{pmatrix}
1 & 0  \\
0 & 1  \\
2 & 1  \\
1 & -1
\end{pmatrix} \,
C_{5\times2}=\begin{pmatrix}
1 & 0  \\
0 & 1  \\
2 & 1  \\
1 & -1  \\
-1 & 2
\end{pmatrix} \,
C_{6\times2}=\begin{pmatrix}
1 & 0  \\
0 & 1  \\
2 & 1  \\
1 & 1  \\
1 & -1 \\
-1 & 2
\end{pmatrix}
\end{equation*}

\begin{equation*}
C_{5\times3} = \begin{pmatrix}
1  &  0   &   0 \\
0  &  1   &   0 \\
0  &  0   &   1 \\
1  &  1   &   -1 \\
-1  &  2   &   1
\end{pmatrix} \,
C_{6\times3} = \begin{pmatrix}
1  &  0   &   0 \\
0  &  1   &   0 \\
0  &  0   &   1 \\
1  &  -1  &   1 \\
1  &  1   &  -1 \\
-1  &  2   &   1
\end{pmatrix} \,
C_{7\times3} = \begin{pmatrix}
1  &  0   &   0 \\
0  &  1   &   0 \\
0  &  0   &   1 \\
2  &  1   &  -1 \\
1  &  1   &   1 \\
1  &  1   &  -1 \\
-1 &  \hf &  2
\end{pmatrix}
\end{equation*}
\begin{equation*}
C_{8\times3} = \begin{pmatrix}
1  &  0   &   0 \\
0  &  1   &   0 \\
0  &  0   &   1 \\
2  &  1   &  -1 \\
1  &  1   &   1 \\
1  &  1   &  -1 \\
\hf & \hf &   1 \\
-1 &  \hf &   2
\end{pmatrix} \,
C_{9\times3} = \begin{pmatrix}
1  &  0   &   0 \\
0  &  1   &   0 \\
0  &  0   &   1 \\
2  &  1   &  -1 \\
\hf&  2   &   1 \\
-1 &  1   &  \hf \\
2  &  -1  &  1 \\
1  &   1  &  -1 \\
-1 &  \hf &   2
\end{pmatrix}
\end{equation*}
\begin{equation*}
C_{7\times4} = \begin{pmatrix}
1  &  0  &  0  & 0 \\
0  &  1  &  0  & 0 \\
0  &  0  &  1  & 0 \\
0  &  0  &  0  & 1 \\
1  &  1  &  1  & 1 \\
1  &  2  &  1  & 2 \\
1  &  -2 &  2  & 1
\end{pmatrix}\,
C_{9\times5} = \begin{pmatrix}
1  &  0  &  0  & 0  & 0 \\
0  &  1  &  0  & 0  & 0 \\
0  &  0  &  1  & 0  & 0 \\
0  &  0  &  0  & 1  & 0 \\
0  &  0  &  0  & 0  & 1 \\

1  &  1  &  1  & 1  & 1 \\
\hf&  1  & \hf & 1  &\hf \\
-1 & -1  & \hf & 1  &2 \\
2  &  1  &-\hf & 2  &-\hf
\end{pmatrix}
\end{equation*}

\bibliographystyle{plain}
\bibliography{biblio.bib}

\end{document}