The Malgrange Form and Fredholm Determinants

The Malgrange Form and Fredholm Determinants

Marco BERTOLA †‡

M. Bertola

† Department of Mathematics and Statistics, Concordia University, Montréal,
Canada<EMAIL_ADDRESS>
‡ Area of Mathematics SISSA/ISAS, Trieste, Italy<EMAIL_ADDRESS>
Received March 12, 2017, in final form June 17, 2017; Published online June
22, 2017

We consider the factorization problem of matrix symbols relative to a closed
contour, i.e., a Riemann–Hilbert problem, where the symbol depends
analytically on parameters. We show how to define a function $\tau$ which is
locally analytic on the space of deformations and that is expressed as a
Fredholm determinant of an operator of “integrable” type in the sense of
Its–Izergin–Korepin–Slavnov. The construction is not unique and the non-
uniqueness highlights the fact that the tau function is really the section of
a line bundle.

Malgrange form; Fredholm determinants; tau function

35Q15; 47A53; 47A68

##  1 Introduction

We shall consider the following prototypical matrix Riemann–Hilbert problem
(RHP) on the unit circle $\Sigma$ (or any smooth closed simple contour):

$\displaystyle\Gamma_{+}(z;\mathbf{t})=\Gamma_{-}(z;\mathbf{t})M(z;\mathbf{t}),\qquad\forall\,z\in\Sigma,\qquad\Gamma(\infty)={\bf
1}.$ (1.1)

Here $\mathbf{t}$ stands for a vector of parameters which we refer to as
“deformation parameters”. The assumptions are the following;

  1. 1.

The matrix $M(z;\mathbf{t})\in{\rm GL}_{n}({\mathbb{C}})$ is jointly analytic
for $z$ in a fixed tubular neighbourhood $N(\Sigma)$ of $\Sigma$ and
$\mathbf{t}$ in an open connected domain $\mathcal{S}$, which we refer to as
the “deformation space”.

  2. 2.

The index of $\det M(z;\mathbf{t})$ around $\Sigma$ vanishes for all
$\mathbf{t}\in\mathcal{S}$.

  3. 3.

The partial indices are generically zero, i.e., the RHP (1.1) generically
admits solution.

Let us remind the reader of some facts that can be extracted from [7]

  1. –

There exists a matrix function $Y_{-}(z)$ analytic and analytically invertible
in $\operatorname{Ext}(\Sigma)\cup N(\Sigma)$ (and uniformly bounded) and
similarly a matrix function $Y_{+}(z)$ analytic and analytically invertible in
$\operatorname{Int}(\Sigma)\cup N(\Sigma)$ and $n$ integers
$k_{1},\dots,k_{n}$ (called partial indices) such that

$\displaystyle D(z)Y_{+}(z)=Y_{-}(z)M(z),\quad z\in\Sigma,\qquad
D=\operatorname{diag}\big{(}z^{k_{1}},\dots,z^{k_{n}}\big{)}.$

  2. –

The RHP (1.1) is solvable if and only if all partial indices vanish,
$k_{j}=0$, $\forall\,j=1,\dots,n$.

Note that since $\operatorname{ind}_{\Sigma}\det
M=\sum\limits_{j=1}^{n}k_{j}$, the condition (2) in our assumptions is
necessary for the solvability of (1.1).

We shall denote by $\mathbb{D}_{\pm}$ the interior ($+$) and the exterior
($-$) regions separated by $\Sigma$. Define ${\cal H}_{+}$ to be the space of
functions that are in $L^{2}(\Sigma,|\mathrm{d}z|)$ and extend to analytic
functions in the interior. We will use the notation $\vec{\cal H}_{+}={\cal
H}_{+}\otimes{\mathbb{C}}^{r}$ (i.e., vector-valued such functions). The
vectors will be thought of as row-vectors. We also introduce the Cauchy
projectors $C_{\pm}\colon L^{2}(\Sigma,|\mathrm{d}z|)\to{\cal H}_{\pm}$:

$\displaystyle
C_{\pm}[f](z)=\oint_{\Sigma}\frac{f(w)\mathrm{d}w}{(w-z)2i\pi},\qquad
z\in\mathbb{D}_{\pm}.$

It is well known [14] that the RHP (1.1) is solvable if and only if the
Toeplitz operator111Due to our choices of symbols, the matrix symbol of the
relevant Toeplitz operator is $M^{-1}$. We apologize for the inconvenience.

$\displaystyle T_{S}\colon\ \vec{\cal H}_{+}\to\vec{\cal H}_{+},$
$\displaystyle T_{S}[\vec{f}]=C_{+}[\vec{f}S],\qquad
S(z;\mathbf{t}):=M^{-1}(z;\mathbf{t})$

is invertible, in which case the inverse is given by

$\displaystyle
T_{S}^{-1}[\vec{f}]=C_{+}\big{[}\vec{f}\Gamma_{-}^{-1}\big{]}\Gamma_{+}.$

Moreover the operator is Fredholm and

$\displaystyle\dim\ker(T_{S})-\dim\operatorname{coker}(T_{S})=\operatorname{ind}_{\Sigma}\det
M=0.$

There is no reasonable way, however, to define a “determinant” of $T_{S}$ as
it stands. Such a function of $\mathbf{t}$ would desirably have the property
that the RHP (1.1) is not solvable if and only if this putative determinant is
zero.

While this is notoriously impossible in this naive form, we now propose a
proxy for the notion of determinant, in terms of a simple Fredholm
determinant.

The Malgrange one-form. As Malgrange explains [15] one can define a central
extension on the loop group $\mathcal{G}:=\\{M\colon\Sigma\to{\rm
GL}_{n}({\mathbb{C}})\colon\operatorname{ind}_{\Sigma}\det M=0\\}$ given by
$\widehat{\mathcal{G}}=\\{(M,u)\in\mathcal{G}\times{\mathbb{C}}^{\times}\\}$
with the group law222We are being a bit cavalier in this description; we
invite the reader to read pp. 1373–1374 in [15].

$\displaystyle(M,u)\cdot(\widetilde{M},\widetilde{u})=\big{(}M\widetilde{M},u\widetilde{u}c\big{(}M,\widetilde{M}\big{)}\big{)},\qquad
c(M,\widetilde{M}):=\det_{{\cal
H}_{+}}\big{(}T_{M^{-1}}T_{\widetilde{M}^{-1}}T^{-1}_{(M\widetilde{M})^{-1}}\big{)}.$

The operator in the determinant is of the form $\operatorname{Id}_{{\cal
H}_{+}}+$ (trace class) and hence the Fredholm determinant is well defined.
This group law is only valid for pairs $M$, $\widetilde{M}$ for which the
inverse of $T^{-1}_{(M\widetilde{M})^{-1}}$ exists. The left-invariant
Maurer–Cartan form of this central extension is then given by
$\big{(}S^{-1}\delta S,\frac{\mathrm{d}u}{u}+\widehat{\omega}_{M}\big{)}$,
where

$\displaystyle\widehat{\omega}_{M}:=\operatorname{Tr}_{{\cal
H}_{+}}\big{(}T_{S}^{-1}\circ T_{\delta S}-T_{S^{-1}\delta S}\big{)},\qquad
S:=M^{-1},$

and this can be written as the following integral

$\displaystyle\widehat{\omega}_{M}=\oint_{\Sigma}\operatorname{Tr}\big{(}\Gamma_{+}^{-1}\Gamma_{+}^{\prime}M^{-1}\delta
M\big{)}\frac{\mathrm{d}z}{2i\pi}.$ (1.2)

Here, and below, $\delta$ denotes the exterior total differentiation in the
deformation space $\mathcal{S}$:

$\displaystyle\delta=\sum\delta t_{j}\frac{\partial}{\partial t_{j}}.$

The Malgrange form is a logarithmic form in the sense that it has only simple
poles on a co-dimension $1$ analytic submanifold of the deformation space
$\mathcal{S}$ and with positive integer Poincaré residue along it; this
manifold is precisely the exceptional “divisor” $(\Theta)\subset\mathcal{S}$
(the Malgrange divisor) where the RHP (1.1) becomes non solvable, i.e., where
some partial indices of the Birkhoff factorization become non-zero.

Closely related to (1.2) is the following one-form, which we still name after
Malgrange:

$\displaystyle\omega_{M}:=\oint_{\Sigma}\operatorname{Tr}\big{(}\Gamma_{-}^{-1}\Gamma_{-}^{\prime}\delta
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}.$ (1.3)

It is also a logarithmic form with the same pole-divisor; indeed one verifies
that

$\displaystyle\widehat{\omega}_{M}-\omega_{M}=\oint_{\Sigma}\operatorname{Tr}\big{(}M^{\prime}M^{-1}\delta
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi},$

which is an analytic form of the deformation parameters
$\mathbf{t}\in\mathcal{S}$. In [2] the one-form (1.3) was posited as an object
of interest for general Riemann–Hilbert problems (not necessarily on closed
contours) and its exterior derivative computed (with an important correction
in [3], which is however irrelevant in the present context). It was computed
(but the computation can be traced back to Malgrange himself in this case)
that

$\displaystyle\delta\omega_{M}=\frac{1}{2}\oint_{\Sigma}\operatorname{Tr}\left(\Xi(z)\wedge\frac{\mathrm{d}}{\mathrm{d}z}\Xi(z)\right)\frac{\mathrm{d}z}{2i\pi},\qquad\Xi(z;\mathbf{t}):=\delta
M(z;\mathbf{t})M(z;\mathbf{t})^{-1}.$ (1.4)

It appears from (1.4) that this two form $\delta\omega_{M}$ is not only
closed, but also smooth on the whole of $\mathcal{S}$, including the
Malgrange-divisor. As such, it defines a line-bundle $\mathcal{L}$ over
$\mathcal{S}$ by the usual construction: one covers $\mathcal{S}$ by
appropriate open sets $U_{\alpha}$ where
$\delta\omega_{M}=\delta\theta_{\alpha}$; on the overlap $U_{\alpha}\cap
U_{\beta}$ the form $\theta_{\alpha}-\theta_{\beta}$ is also exact and one
defines then the transition functions by
$g_{\alpha\beta}(\mathbf{t})=\exp\big{(}\int\theta_{\alpha}-\theta_{\beta}\big{)}$.
Then a section of this line bundle is provided by the collection of functions
$\tau_{\alpha}\colon U_{\alpha}\to{\mathbb{C}}$ such that

$\displaystyle\tau_{\alpha}(\mathbf{t})=\exp\left[\int(\omega_{M}-\theta_{\alpha})\right].$

Since $\omega_{M}$ is a logarithmic form and each $\theta_{\alpha}$ is
analytic in the respective $U_{\alpha}$, the functions
$\tau_{\alpha}(\mathbf{t})$ have zero of finite order precisely on the
Malgrange divisor $(\Theta)\subset\mathcal{S}$ (under appropriate
transversality assumptions, the order of the zero is the dimension of
$\operatorname{Ker}T_{S}$).

Our goal is to provide an explicit construction of the $\tau_{\alpha}$’s in
terms of Fredholm determinants of simple operators of the
Its–Izergin–Korepin–Slavnov “integrable” type [13]. Their definition is
recalled in due time.

##  2 Construction of the Fredholm determinants

The construction carried out below is not unique, and also only local in the
deformation space $\mathcal{S}$; this is however not only not a problem, but
rather an interesting feature, as we will illustrate in the case of ${\rm
SL}_{2}({\mathbb{C}})$. The non-uniqueness is precisely a consequence of the
fact that we are trying to compute a section of the aforementioned line
bundle.

Preparatory step. The assumption that $M(z;\mathbf{t})\in{\rm GL}_{n}$ can be
replaced without loss of generality with $M(z;\mathbf{t})\in{\rm SL}_{n}$;
this is so because of the assumption on the index of $\det M$. Indeed we can
solve the scalar problem $y_{+}(z)=y_{-}(z)\det M(z)$, $y(\infty)=1$ and then
define a new RHP where
$\widetilde{\Gamma}_{\pm}(z):=\Gamma_{\pm}\operatorname{diag}\big{(}y_{\pm}^{-1},1,\dots\big{)}$
and hence the new matrix jump for $\widetilde{\Gamma}$ is
$\widetilde{M}(z)=\operatorname{diag}(y_{-}(z),1,\dots)M(z)\operatorname{diag}\big{(}y_{+}^{-1}(z),1,\dots\big{)}$
with $\det\widetilde{M}\equiv 1$. For this reason, from here on we assume
$M\in{\rm SL}_{n}({\mathbb{C}})$.

We define an elementary matrix (for our purposes) to be a matrix of the form
${\bf 1}+cE_{jk}$, with $j\neq k$, where $E_{jk}$ denotes the $(j,k)$-unit
matrix.

######  Lemma 2.1.

Any matrix $M\in{\rm SL}_{n}({\mathbb{C}})$ can be written as a product of
elementary matrices. The entries of the factorization are rational in the
entries of $M$ with denominators that are monomials in a suitable set of $n-1$
nested minors of $M$.

###### Proof.

We recall that given any matrix $M\in{\rm SL}_{n}$, there is a permutation
$\Pi$ of the columns such that the principal minors (the determinants of the
top left square submatrices) do not vanish, and hence we can write it as

$\displaystyle M=LDU\Pi=\widehat{M}\Pi,$

where $L$, $U$ are lower/upper triangular matrices with unit on the diagonal
and $D=\operatorname{diag}(x_{1},\dots$, $x_{n})$ is a diagonal matrix (see
for example [8, Vol. 1, Chapter II]). Denote
$q_{\ell}=\det\big{[}\widehat{M}_{j,k}\big{]}_{j,k\leq\ell}$ the principal
minors; these are the nested minors of the original matrix $M$ alluded to in
the statement. The matrices $L$, $U$ are rational in the entries of $M$ and
with denominators that are monomials in the $q_{\ell}$’s.

Now, both $L$, $U$ can clearly be written as products of elementary matrices
whose coefficients are polynomials in the entries of $L$, $U$ (respectively)
and so it remains to show that we can write $D$ as product of elementary
matrices.

To this end we observe the ’LULU’ identity (there is a similar ‘ULUL’
identity)

$\displaystyle\left[\begin{matrix}x&0\\\
0&\frac{1}{x}\end{matrix}\right]=\left[\begin{matrix}1&0\\\
\frac{1-x}{x}&1\end{matrix}\right]\left[\begin{matrix}1&1\\\
0&1\end{matrix}\right]\left[\begin{matrix}1&0\\\
x-1&1\end{matrix}\right]\left[\begin{matrix}1&-\frac{1}{x}\\\
0&1\end{matrix}\right].$

As $D\in{\rm SL}_{n}$, we can represent it in terms of product of embedded
${\rm SL}_{2}$ matrices using the root decomposition of ${\rm SL}_{n}$:

$\displaystyle
D=\operatorname{diag}\left(x_{1},\frac{1}{x}_{1},1,\dots,1\right)\operatorname{diag}\left(1,x_{2}x_{1},\frac{1}{x_{1}x_{2}},1,\dots,1\right)$
$\displaystyle\hphantom{D=}{}\times\operatorname{diag}\left(1,1,x_{1}x_{2}x_{3},\frac{1}{x_{1}x_{2}x_{3}},1,\dots,1\right)\cdots,$

and then embed the LULU identity for each factor.

Finally, also permutation matrices can be written as product of elementary
matrices embedding appropriately the simple identity;

$\displaystyle\left[\begin{matrix}0&-1\\\
1&0\end{matrix}\right]=\left[\begin{matrix}1&0\\\
1&1\end{matrix}\right]\left[\begin{matrix}1&-1\\\
0&1\end{matrix}\right]\left[\begin{matrix}1&0\\\ 1&1\end{matrix}\right].$

This concludes the proof. ∎

Let now $M(z)$ be an ${\rm SL}_{n}$ matrix valued function, analytic in a
tubular neighbourhood $N(\Sigma)$ of $\Sigma$, and let $q_{\ell}(z)$,
$\ell=1,\dots,n-1$ be the nested minors alluded to in the Lemma so that they
are not identically zero. Since the entries are analytic in $N(\Sigma)$ we can
slightly deform the contour $\Sigma$ to a contour $\widetilde{\Sigma}$ that
avoids all zeroes of every principal minor $q_{\ell}(z)$. The resulting RHP is
“equivalent” to the original in the sense that the solvability of one implies
the solvability of the other. Note also that this deformation can be done in a
piecewise constant way locally with respect to $\mathbf{t}\in\mathcal{S}$.
Thus we have

$\displaystyle M(z)=F_{1}(z)\cdots F_{R}(z),\qquad F_{\nu}(z)={\bf
1}+a_{\nu}(z)E_{j_{\nu},k_{\nu}},\qquad\nu=1,\dots,R.$

$\Sigma=\Sigma_{1}$$\mathbb{D}_{1}$$\Sigma_{2}$$\mathbb{D}_{2}$$\Sigma_{3}$$\mathbb{D}_{R-1}$$\Sigma_{R}$$\mathbb{D}_{R}=\mathbb{D}_{+}$$\mathbb{D}_{0}=\mathbb{D}_{-}$

Figure 1: An illustration of the splitting of the jump matrix into elementary
jumps. The matrix $M(z;\mathbf{t})$ is analytic in the shaded regions.

Corresponding to this factorization we can define an equivalent RHP with jumps
on $R$ contours $\Sigma_{1},\dots,\Sigma_{R}$, with $\Sigma_{1}=\Sigma$ and
$\Sigma_{j+1}$ in the interior of $\Sigma_{j}$ and all of them in the joint
domain of analyticity of the scalar functions $a_{\nu}(z)$ (see Fig. 1) which
may have poles only at the zeroes of the principal minors $q_{\ell}(z)$ of
$M(z)$. This is accomplished by “extending” the matrix $\Gamma_{-}(z)$ to the
annular regions $\mathbb{D}_{0}=\mathbb{D}_{-}$ and
$\mathbb{D}_{j}=\operatorname{Int}(\Sigma_{j})\cap\operatorname{Ext}(\Sigma_{j+1})$
as

$\displaystyle\Theta_{0}(z):=\Gamma_{-}(z),\qquad\forall\,z\in\mathbb{D}_{0},$
$\displaystyle\Theta_{\nu}(z):=\Gamma_{-}(z)\overrightarrow{\prod_{\ell=1}^{\nu}}F_{\ell}(z),\qquad\forall\,z\in\mathbb{D}_{\nu}.$
(2.1)

By doing so we obtain the following relations

$\displaystyle\Theta_{\nu}(z)=\Theta_{\nu-1}(z)F_{\nu}(z),\qquad\forall\,z\in\Sigma_{\nu},\qquad\nu=1,\dots,R.$

The piecewise analytic matrix function $\Theta(z)$ whose restriction to
$\mathbb{D}_{\nu}$ coincides with the matrices $\Theta_{\nu}$ (2.1), satisfies
a final RHP

$\displaystyle\Theta_{+}(z)=\Theta_{-}(z)F_{\nu}(z),\qquad\forall\,z\in\Sigma_{\nu},\qquad\Theta(\infty)={\bf
1}.$ (2.2)

This type of RHP is of the general type of “integrable kernels” and its
solvability can be determined by computing the Fredholm determinant of an
integral operator of
$L^{2}(\bigsqcup\Sigma_{\nu},|\mathrm{d}z|)\simeq\bigoplus_{\nu=1}^{R}L^{2}(\Sigma_{\nu},|\mathrm{d}z|)$
with kernel (we use the same symbol for the operator and its kernel)

$\displaystyle
K(z,w)=\frac{\vec{f}^{T}(z)\vec{g}(w)}{2i\pi(w-z)},\qquad\vec{f}(z)=\sum_{\nu=1}^{R}{\bf
e}_{j_{\nu}}\chi_{\nu}(z)a_{\nu}(z),\qquad\vec{g}(z)=\sum_{\nu=1}^{R}{\bf
e}_{k_{\nu}}\chi_{\nu}(z),$ (2.3)

where $\chi_{\nu}(z)$ is the projector (indicator function) on the component
$L^{2}(\Sigma_{\nu},|\mathrm{d}z|)$. Indeed, as explained in [4, 10, 13], the
Fredholm determinant $\det(\operatorname{Id}-K)$ is zero if and only if the
RHP (2.2) is non-solvable and moreover the resolvent operator
$R=K(\operatorname{Id}-K)^{-1}$ of $K$ has kernel

$\displaystyle
R(z,w)=\frac{\vec{f}^{t}(z)\Theta^{t}(z)(\Theta^{t})^{-1}(w)\vec{g}(w)}{z-w}.$

######  Theorem 2.2.

The RHP (2.2) and hence (1.1) is solvable if and only if
$\tau:=\det(\operatorname{Id}-K)\neq 0$.

######  Proposition 2.3 (see, e.g., [4, Theorem 2.1]).

Let $\partial$ be any deformation of the functions $a_{\nu}(z)$, then

$\displaystyle\partial\ln\tau=\sum_{\nu=1}^{R}\oint_{\Sigma_{\nu}}\operatorname{Tr}\big{(}\Theta_{-}^{-1}\Theta_{-}^{\prime}(z)\partial
F_{\nu}F_{\nu}^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}=\sum_{\nu=1}^{R}\oint_{\Sigma_{\nu}}\big{(}\Theta_{\nu-1}^{-1}\Theta_{\nu-1}^{\prime}(z)\big{)}_{k_{\nu},j_{\nu}}\partial
a_{\nu}(z)\frac{\mathrm{d}z}{2i\pi}.$

##  3 The $\boldsymbol{{\rm SL}_{2}}$ case

We would like to express the Malgrange one-form directly in terms of the
$\tau$ function (Fredholm determinant). Rather than obscuring the simple idea
with the general case, we consider in detail the ${\rm SL}_{2}$ case. Let
$M(z;\mathbf{t})$ be analytic in $(z,\mathbf{t})\in
N(\Sigma)\times\mathcal{S}$ and with values in ${\rm SL}_{2}({\mathbb{C}})$.
Using the general scheme above, we have the following factorizations

$\displaystyle{(1)}\ \ \left[\begin{matrix}a&b\\\
c&d\end{matrix}\right]=\underbrace{\left[\begin{matrix}1&0\\\
\frac{1+c-a}{a}&1\end{matrix}\right]}_{F_{1}(z)}\underbrace{\left[\begin{matrix}1&1\\\
0&1\end{matrix}\right]}_{F_{2}(z)}\underbrace{\left[\begin{matrix}1&0\\\
a-1&1\end{matrix}\right]}_{F_{3}(z)}\underbrace{\left[\begin{matrix}1&\frac{b-1}{a}\\\
0&1\end{matrix}\right]}_{F_{4}(z)},\qquad a\not\equiv 0,$ $\displaystyle{(2)}\
\ \left[\begin{matrix}a&b\\\
c&d\end{matrix}\right]=\underbrace{\left[\begin{matrix}1&\frac{1+b-d}{d}\\\
0&1\end{matrix}\right]}_{F_{1}(z)}\underbrace{\left[\begin{matrix}1&0\\\
1&1\end{matrix}\right]}_{F_{2}(z)}\underbrace{\left[\begin{matrix}1&d-1\\\
0&1\end{matrix}\right]}_{F_{3}(z)}\underbrace{\left[\begin{matrix}1&0\\\
\frac{c-1}{d}&1\end{matrix}\right]}_{F_{4}(z)},\qquad d\not\equiv 0,$
$\displaystyle{(3)}\ \ \left[\begin{matrix}a&b\\\
-\frac{1}{b}&0\end{matrix}\right]=\underbrace{\left[\begin{matrix}1&b-ab\\\
0&1\end{matrix}\right]}_{F_{1}(z)}\underbrace{\left[\begin{matrix}1&0\\\
-\frac{1}{b}&1\end{matrix}\right]}_{F_{2}(z)}\underbrace{\left[\begin{matrix}1&b\\\
0&1\end{matrix}\right]}_{F_{3}(z)},$ $\displaystyle{(4)}\ \
\left[\begin{matrix}0&b\\\
-\frac{1}{b}&d\end{matrix}\right]=\underbrace{\left[\begin{matrix}1&\frac{d-1}{b}\\\
0&1\end{matrix}\right]}_{F_{1}(z)}\underbrace{\left[\begin{matrix}1&0\\\
b&1\end{matrix}\right]}_{F_{2}(z)}\underbrace{\left[\begin{matrix}1&-\frac{1}{b}\\\
0&1\end{matrix}\right]}_{F_{3}(z)},$ $\displaystyle{(5)}\ \
\left[\begin{matrix}0&b\\\
-\frac{1}{b}&0\end{matrix}\right]=\underbrace{\left[\begin{matrix}1&0\\\
b&1\end{matrix}\right]}_{F_{1}(z)}\underbrace{\left[\begin{matrix}1&-\frac{1}{b}\\\
0&1\end{matrix}\right]}_{F_{2}(z)}\underbrace{\left[\begin{matrix}1&0\\\
b&1\end{matrix}\right]}_{F_{3}(z)}.$ (3.1)

###  3.1 Fredholm determinants for different factorizations

Each of the factorization (3.1) leads to an integrable operator of the form
(2.3) and hence to a corresponding Fredholm determinant; we now establish
their mutual relationships.

There are two types of questions that we address here

  1. 1.

How are the Fredholm determinants associated with the different factorizations
$(1,2)$ in (3.1) related to each other?

  2. 2.

For a fixed factorization, how does the Fredholm determinant depend on the
choice of contour $\Sigma$ (within the analyticity domain $N(\Sigma)$).

Consider the cases (3.1)(ρ), $\rho=1,\dots,5$. We compute the logarithmic
derivative of the corresponding Fredholm determinant using Proposition 2.3

$\displaystyle\partial\ln\tau_{(\rho)}=\sum_{\nu=1}^{R}\oint_{\Sigma_{\nu}}\operatorname{Tr}\big{(}\Theta_{-}^{-1}\Theta_{-}^{\prime}\partial
F_{\nu}F_{\nu}^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}.$ (3.2)

By using the relationship (2.1) between $\Gamma_{-}$ (extended to an analytic
function on $\operatorname{Ext}(\Sigma)\cup N(\Sigma)$), we can re-express it
in terms of the Malgrange one-form of the original problem (1.1); we use the
fact that we can deform the contours back to $\Sigma=\Sigma_{1}$ by Cauchy’s
theorem. Plugging (2.1) appropriately in (3.2) and using Leibnitz rule, after
a short computation we obtain

$\displaystyle\delta\ln\tau_{{}_{(\rho)}}=\oint_{\Sigma}\operatorname{Tr}\big{(}\Gamma_{-}^{-1}\Gamma_{-}^{\prime}\partial
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}$ (3.3)
$\displaystyle\hphantom{\delta\ln\tau_{{}_{(\rho)}}=}{}+\\!\oint_{\Sigma}\\!\operatorname{Tr}\big{(}F_{1}^{-1}F_{1}^{\prime}\partial
F_{2}F_{2}^{-1}+F_{12}^{-1}F_{12}^{\prime}\partial
F_{3}F_{3}^{-1}+F_{123}^{-1}F_{123}^{\prime}\partial
F_{4}F_{4}^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}=:\omega_{M}+\theta_{(\rho)},\\!$

where $F_{1\dots k}=F_{1}F_{2}F_{3}\cdots F_{k}$ (if the factorization has
only three term, then we set $F_{4}\equiv{\bf 1}$).

The cases (3,4,5). The last cases (3.1)(3,4,5) lead essentially to a RHP with
a triangular jump; it suffices to re-define $\Gamma$ by
$\Gamma\left[\begin{smallmatrix}0&1\\\ -1&0\end{smallmatrix}\right]$ for
$z\in\operatorname{Int}(\Sigma)$. If the index of $b$ is zero,
$\operatorname{ind}_{\Sigma}b=0$,then the solution can be written explicitly
in closed form and it is interesting to compute the Fredholm determinant
associated to our factorization of the matrix. We will show

######  Proposition 3.1.

In the cases (3.1)(3,4,5) and under the additional assumption that
$\operatorname{ind}_{\Sigma}b=0$, the $\tau$ function given by
$\det(\operatorname{Id}_{L^{2}(\cup\Sigma_{j})}-K)$ and $K$ as in (2.3),
equals the constant in the strong Szegö formula, given by [1, 5]

$\displaystyle\tau=\exp\bigg{[}\sum_{j>0}j\beta_{-j}\beta_{j}\bigg{]}=\det_{{\cal
H}_{+}}T_{b}T_{b^{-1}},$

where now the Toeplitz operator is for the scalar symbol $b(z)$ and ${\cal
H}_{+}$ is the Hardy space of scalar functions analytic in $\mathbb{D}_{+}$
and $\beta_{j}$ are the coefficients of $\ln b(z)$ in the Laurent expansion
centered at the origin

$\displaystyle\beta(z):=\ln b(z)=\sum_{j\in{\mathbb{Z}}}\beta_{j}z^{j}.$ (3.4)

The same applies in the case that the jump is triangular $(b\equiv 0$ and/or
$c\equiv 0)$ under the assumption $\operatorname{ind}_{\Sigma}a=0$ and
replacing $b$ with $a$ in the above formulas.

###### Proof.

From a direct computation of the term $\theta_{(\rho)}$ in (3.3), we find

$\displaystyle\theta_{{}_{(3,4,5)}}=\oint_{\Sigma}\left(\delta\beta\frac{\mathrm{d}}{\mathrm{d}z}\beta\right)\frac{\mathrm{d}z}{2i\pi},\qquad\beta(z):=\ln
b(z)$

and the solution of the RHP is explicit,

$\displaystyle\Gamma(z)=\begin{cases}\displaystyle\left[\begin{matrix}\displaystyle{\rm
e}^{B(z)}&-\displaystyle{\rm e}^{-B(z)}\oint_{\Sigma}\frac{a(w)b(w){\rm
e}^{2B_{-}(w)}}{w-z}\frac{\mathrm{d}w}{2i\pi}\\\ 0&{\rm
e}^{-B(z)}\end{matrix}\right],&z\in\mathbb{D}_{-},\\\
\displaystyle\left[\begin{matrix}\displaystyle{\rm
e}^{B(z)}&-\displaystyle{\rm e}^{-B(z)}\oint_{\Sigma}\frac{a(w)b(w){\rm
e}^{2B_{-}(w)}}{w-z}\frac{\mathrm{d}w}{2i\pi}\\\ 0&{\rm
e}^{-B(z)}\end{matrix}\right]\left[\begin{matrix}0&1\\\
-1&0\end{matrix}\right],&z\in\mathbb{D}_{+},\end{cases}$ $\displaystyle
B(z):=\oint_{\Sigma}\frac{\beta(w)}{w-z}\frac{\mathrm{d}w}{2i\pi}.$

Thus we can write explicitly the Malgrange form:

$\displaystyle\omega_{M}=\oint_{\Sigma}\operatorname{Tr}\big{(}\Gamma_{-}^{-1}\Gamma_{-}^{\prime}\delta
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}=\oint_{\Sigma}\operatorname{Tr}\left(\sigma_{3}\oint_{\Sigma}\frac{\beta(w)}{(w-z_{-})^{2}}\frac{\mathrm{d}w}{2i\pi}\delta\beta(z)\sigma_{3}\right)\frac{\mathrm{d}z}{2i\pi}$
$\displaystyle\hphantom{\omega_{M}}{}=2\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\beta(w)\delta(\beta(z))}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}.$

In this integral, $z$ is integrated on a slightly “larger” contour
$\Sigma_{-}$. The Fredholm determinant satisfies

$\displaystyle\delta\ln\tau_{(3,4,5)}=\omega_{M}+\theta_{{}_{(3,4,5)}}=2\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\beta(w)\delta(\beta(z))}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}+\oint_{\Sigma}\big{[}\beta(\delta\beta)^{\prime}\big{]}\frac{\mathrm{d}w}{2i\pi}.$

We now show that

$\displaystyle\delta\ln\tau=\delta\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\beta(w)\beta(z)}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}=\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\delta\beta(w)\beta(z)+\beta(w)\delta\beta(z)}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}.$

Indeed, the exchange of order of integration of one of the addenda (and
relabeling the variables) yields the other term plus the residue on the
diagonal,

$\displaystyle\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\delta\beta(w)\beta(z)+\beta(w)\delta\beta(z)}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}$
$\displaystyle\qquad{}=\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{2\beta(w)\delta\beta(z)}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}+\oint_{\Sigma}{\beta(z)(\delta\beta)^{\prime}}\frac{\mathrm{d}z}{2i\pi}.$

Therefore, in conclusion, we have (we fix the overall constant of $\tau$ by
requiring it to be $1$ for $b\equiv 1$)

$\displaystyle\ln\tau=\oint_{\Sigma_{-}}\oint_{\Sigma}\frac{\ln b(w)\ln
b(z)}{(w-z)^{2}}\frac{\mathrm{d}w}{2i\pi}\frac{\mathrm{d}z}{2i\pi}.$ (3.5)

If we write a Laurent expansion of $\ln b$ (3.4) the formula (3.5) gives the
explicit expression

$\displaystyle\ln\tau=\sum_{j>0}j\beta_{-j}\beta_{j},$

which is also the formula for the second Szegö limit theorem for the limit of
the Toeplitz determinants of the symbol $b(z)$, and it is known to to be the
Fredholm determinant of an operator [1, 5]

$\displaystyle\tau=\det_{{\cal H}_{+}}T_{b}T_{b^{-1}},$

where now the Toeplitz operator is for the scalar symbol $b(z)$ and ${\cal
H}_{+}$ is the Hardy space of scalar functions analytic in $\mathbb{D}_{+}$. ∎

######  Remark 3.2.

The index assumption is only necessary for the case (3.1)(5) or (3.1)(1) when
$b\equiv 0\equiv c$, because in these situations the RHP separates into two
scalar problems. However, the assumption $\operatorname{ind}_{\Sigma}b=0$ for
cases (3,4) or $\operatorname{ind}_{\Sigma}a=0$ for the triangular case is not
necessary as we now show (in the latter form). Consider the RHP

$\displaystyle Y_{+}=Y_{-}\left[\begin{matrix}1&\mu(z)\\\
0&1\end{matrix}\right],\quad z\in\Sigma,\qquad Y(z)=\big{(}{\bf
1}+\mathcal{O}\big{(}z^{-1}\big{)}\big{)}z^{n\sigma_{3}},\quad z\to\infty.$

By defining

$\displaystyle\Gamma(z)=\begin{cases}Y(z),&z\in\operatorname{Int}(\Sigma),\\\
Y(z)z^{-n\sigma_{3}},&z\in\operatorname{Ext}(\Sigma),\end{cases}$

we are lead to the RHP in standard form

$\displaystyle\Gamma_{+}=\Gamma_{-}\left[\begin{matrix}z^{n}&z^{n}\mu(z)\\\
0&z^{-n}\end{matrix}\right],\quad z\in\Sigma,\qquad\Gamma(\infty)={\bf 1},$

which is case (1) with $c\equiv 0$ (or essentially cases (3,4) up to a
multiplication by piecewise constant matrices). Even if $a(z)=z^{n}$ and thus
$\operatorname{ind}_{\Sigma}a=n$, this problem is still generically solvable
in terms of appropriate “orthogonal polynomials”
$\\{p_{\ell}(z)\\}_{\ell\in{\mathbb{N}}}$ defined by the “orthogonality”
property

$\displaystyle\oint_{\Sigma}p_{\ell}(z)p_{k}(z)\mu(z)\mathrm{d}z=h_{\ell}\delta_{\ell
k}.$

Then the solution of the $Y$-problem is written as

$\displaystyle
Y(z)=\left[\begin{matrix}p_{n}(z)&\displaystyle\oint_{\Sigma}\frac{p_{n}(w)\mu(w)\mathrm{d}w}{(w-z)2i\pi}\vspace{1mm}\\\
\displaystyle\frac{-2i\pi
p_{n-1}(z)}{h_{n-1}}&\displaystyle\oint_{\Sigma}\frac{-p_{n-1}(w)\mu(w)\mathrm{d}w}{(w-z)h_{n-1}}\end{matrix}\right],$

and the solvability depends only on the condition [6]

$\displaystyle\det\left[\oint_{\Sigma}z^{\ell+j-2}\mu(z)\mathrm{d}z\right]_{j,\ell=1}^{n}\neq
0.$

The cases (1,2). A straightforward computation using the explicit expression
(3.1)(1) yields

$\displaystyle\partial\ln\tau_{{}_{(1)}}=\oint_{\Sigma}\Big{[}\operatorname{Tr}\big{(}\Gamma_{-}^{-1}\Gamma_{-}^{\prime}\partial
MM^{-1}\big{)}+\partial\ln a(\ln a)^{\prime}(1+bc)$
$\displaystyle\hphantom{\partial\ln\tau_{{}_{(1)}}=\oint_{\Sigma}\Big{[}}{}-c\partial
b(\ln a)^{\prime}+c^{\prime}\partial b-c^{\prime}b\partial(\ln
a)\Big{]}\frac{\mathrm{d}z}{2i\pi}$
$\displaystyle\hphantom{\partial\ln\tau_{{}_{(1)}}}{}=\oint_{\Sigma}\operatorname{Tr}\big{(}\Gamma^{-1}\Gamma^{\prime}\partial
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}+\underbrace{\oint_{\Sigma}\left(\frac{(a)^{\prime}}{a}(d\partial
a-c\partial
b)+{ac^{\prime}}\partial\left(\frac{b}{a}\right)\right)\frac{\mathrm{d}z}{2i\pi}}_{\theta_{{}_{(1)}}}.$
(3.6)

Since $\partial\ln\tau$ is a closed differential, the exterior derivative of
$\theta_{{}_{(1)}}$ must be opposite to the one of the first term, which is
given by (1.4). Let us verify this directly; to this end we compute the
exterior derivative of $\theta_{{}_{(1)}}$. A straightforward computation
yields

$\displaystyle\delta\theta_{{}_{(1)}}=\oint_{\Sigma}\left(\delta
b\wedge(\delta c)^{\prime}+\frac{\delta a\wedge(\delta
a)^{\prime}}{a^{2}}(1+bc)+\delta a\wedge\left(\frac{\delta
c}{a}\right)^{\prime}{b}\right.$
$\displaystyle\left.\hphantom{\delta\theta_{{}_{(1)}}=\oint_{\Sigma}}{}-\delta
a\wedge\delta b\frac{(c)^{\prime}}{a}-\delta b\wedge(\delta
a)^{\prime}\frac{c}{a}-\delta b\wedge\delta
c\frac{a^{\prime}}{a}\right)\frac{\mathrm{d}z}{2i\pi}.$

The exterior derivative of $\omega_{M}$ is given by (1.4), in which we can
insert the explicit expression of $M(z;\mathbf{t})$; after a somewhat lengthy
but straightforward computation we find that

$\displaystyle\delta\omega_{M}+\delta\theta_{{}_{(1)}}=-\frac{1}{2}\oint_{\Sigma}\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{c}{a}\delta
a\wedge\delta b-\frac{b}{a}\delta a\wedge\delta c+\delta b\wedge\delta
c\right)\frac{\mathrm{d}z}{2i\pi}=0,$

thus confirming that the differential $\delta\ln\tau_{{}_{(1)}}$ is indeed
closed (of course it must be, since the tau function is a Fredholm
determinant!). The case ${(2)}$ is analogous with the replacements
$b\leftrightarrow c$ and $a\leftrightarrow d$.

####  3.1.1 Determinants for different choices of $\boldsymbol{\Sigma}$

The factorizations (3.1)(1,2) require that we deform the contour $\Sigma$ so
that $a(z)$ (or $d(z)$) does not have any zero on $\Sigma$. This leads to
completely equivalent RHPs of the form (1.1) but not entirely equivalent RHP
when expressed in the form (2.2).

Note that (3.1)(3,4,5) do not suffer from this ambiguity, because $b(z)$
cannot have any zeroes in $N(\Sigma)$ since we assumed analyticity of the jump
matrix $M(z;\mathbf{t})$.

Consider the factorization (3.1)(1) (with similar considerations applying to
the other factorization); in general $a(z;\mathbf{t})$ has zeroes in its
domain of analyticity, and their positions depend on $\mathbf{t}$. Therefore
it may be necessary, when considering the dependence on $\mathbf{t}$, to move
the contour so that certain zeroes are to the left or to the right of it
because a zero may sweep across $N(\Sigma)$ as we vary
$\mathbf{t}\in\mathcal{S}$.

$\widetilde{\Sigma}$$\Sigma$ Figure 2: The two contours and the zeroes of
$a(z;\mathbf{t})$ within the enclosed region. The shaded area is $N(\Sigma)$
where $M(z;\mathbf{t})$ is analytic (uniformly w.r.t.
$\mathbf{t}\in\mathcal{S}$).

So, let $\Sigma$, $\widetilde{\Sigma}$ be two contours in the common domain of
analyticity of $M(z;\mathbf{t})$ and such that $a(z;\mathbf{t})$ has no zeroes
on either one and $\Sigma\subset\operatorname{Int}(\widetilde{\Sigma})$ (see
Fig. 2).

Denote with $\tau_{{}_{(1)}}$ and $\widetilde{\tau}_{{}_{(1)}}$ the
corresponding Fredholm determinants of the operators defined as described;
following the same steps as above. Our goal is to show that

######  Theorem 3.3.

The ratio of the two Fredholm determinants is given by

$\displaystyle\widetilde{\tau}_{{}_{(1)}}={\tau_{{}_{(1)}}}\prod_{v\in\operatorname{Int}(\widetilde{\Sigma})\cap\operatorname{Ext}(\Sigma):\atop
a(v)=0}(c(v(\mathbf{t});\mathbf{t}))^{-{{\rm ord}}_{v}(a)}.$

Note that the evaluation of $c$ at the zeroes of $a$ cannot vanish because
$\det M\equiv 1$.

###### Proof.

From the formula (3.6) we get

$\displaystyle\partial\ln\frac{\widetilde{\tau}_{{}_{(1)}}}{\tau_{{}_{(1)}}}=\left(\oint_{\widetilde{\Sigma}}-\oint_{\Sigma}\right)\operatorname{Tr}\big{(}\Gamma_{-}^{-1}\Gamma_{-}^{\prime}\partial
MM^{-1}\big{)}\frac{\mathrm{d}z}{2i\pi}$
$\displaystyle\hphantom{\partial\ln\frac{\widetilde{\tau}_{{}_{(1)}}}{\tau_{{}_{(1)}}}=}{}+\left(\oint_{\widetilde{\Sigma}}-\oint_{\Sigma}\right)\left(\frac{(a)^{\prime}}{a}(d\partial
a-c\partial
b)+{ac^{\prime}}\partial\left(\frac{b}{a}\right)\right)\frac{\mathrm{d}z}{2i\pi}.$

Here $\Gamma_{-}$ means the analytic extension of the solution $\Gamma$ to the
region $\operatorname{Ext}(\Sigma)\cup N(\Sigma)$. Since the integrand of the
first term is holomorphic in the region bounded by $\widetilde{\Sigma}$ and
$\Sigma$, it yields a zero contribution by the Cauchy’s theorem and we are
left only with the second term, which is computable by the residue theorem;

$\displaystyle\partial\ln\frac{\widetilde{\tau}_{{}_{(1)}}}{\tau_{{}_{(1)}}}=\sum_{v\in\operatorname{Int}(\widetilde{\Sigma})\cap\operatorname{Ext}(\Sigma):\atop
a(v)=0}\mathop{\mathrm{res}}_{z=v}\left(\frac{(a)^{\prime}}{a}(d\partial
a-c\partial b)+{c^{\prime}}a\partial\left(\frac{b}{a}\right)\right).$

We are assuming that $M(z;\mathbf{t})$ is analytic and also that $\det M\equiv
1$. Now note that a zero $v(\mathbf{t})$ of $a(z;\mathbf{t})$ in general
depends on $\mathbf{t}$; suppose that
$a(z;\mathbf{t})=(z-v(\mathbf{t}))^{k}(C_{k}(\mathbf{t})+\mathcal{O}(z-v(\mathbf{t})))$;
then

$\displaystyle\frac{\partial
a(z;\mathbf{t})}{a(z;\mathbf{t})}=-\frac{k\partial
v(\mathbf{t})}{z-v(\mathbf{t})}+\mathcal{O}(1),$

and hence the residue evaluation at $z=v(\mathbf{t})$ yields (we use
$0=a^{\prime}\partial v+\partial a|_{v}$ and $\partial(ad)|_{v}=d\partial
a|_{v}=\partial(bc)|_{v}$)

$\displaystyle\mathop{\mathrm{res}}_{z=v}\left(\frac{(a)^{\prime}}{a}\left(d\partial
a-c\partial
b\right)+{ac^{\prime}}\partial\left(\frac{b}{a}\right)\right)=k\left({-}da^{\prime}\big{|}_{v}\partial
v-c\partial b\big{|}_{v}\right)+kbc^{\prime}\partial v$
$\displaystyle\qquad{}=k\left(d\partial a\big{|}_{v}-c\partial
b\big{|}_{v}\right)+kbc^{\prime}\partial
v=k\left(\partial(bc)\big{|}_{v}-c\partial
b\big{|}_{v}\right)+kbc^{\prime}\partial v$ $\displaystyle\qquad{}=kb\partial
c\big{|}_{v}+kbc^{\prime}\big{|}_{v}\partial
v\mathop{=}^{bc|_{v}=-1}-k\partial\ln c(v(\mathbf{t});\mathbf{t})$

(note that $c(v(\mathbf{t});\mathbf{t})$ cannot be zero because $\det M\equiv
1$ and $z=v$ is already a zero of $a$). ∎

####  3.1.2 Determinants for different factorizations

The tau functions (Fredholm determinants) defined thus far should be
understood as defining a section of a line bundle over the loop group space;
this is the line bundle associated to the two fom $\delta\omega_{M}$ (1.4).

This simply means that on the intersection of the open sets where the
factorizations (1), (2) in (3.1) can be made, we have

$\displaystyle\delta\ln\frac{\tau_{{}_{(1)}}}{\tau_{{}_{(2)}}}=\delta\ln\big{(}\Upsilon_{(1,2)}\big{)},$

where

$\displaystyle\delta\ln\big{(}\Upsilon_{(1,2)}\big{)}:=\delta\big{(}\theta_{{}_{(1)}}-\theta_{{}_{(2)}}\big{)}.$

After a short computation we obtain

$\displaystyle\delta\ln\big{(}\Upsilon_{(1,2)}\big{)}=\oint_{\Sigma}\left(\frac{c^{\prime}\delta
b-b^{\prime}\delta
c}{1+bc}\right)\frac{\mathrm{d}z}{2i\pi}=\oint_{\Sigma}\left(\delta\ln
b\frac{\mathrm{d}}{\mathrm{d}z}\ln(ad)-\delta\ln(ad)\frac{\mathrm{d}}{\mathrm{d}z}\ln
b\right)\frac{\mathrm{d}z}{2i\pi}.$

Observe the last expression; in principle the functions $b(z)$, $a(z)d(z)$ may
have nonzero index around $\Sigma$, but in any case the functions $\delta\ln
b$ and $(\ln(ad))^{\prime}$ are single-valued because the increments of the
logarithms are integer multiples of $2i\pi$ and hence locally constant in the
space of deformations. To write explicitly the transition function (or rather
a representative of the same cocycle class) we need to choose a point
$z_{0}\in\operatorname{Int}(\Sigma)$; we choose $z=0$ without loss of
generality. Let $K=\operatorname{ind}_{\Sigma}b$,
$L=\operatorname{ind}_{\Sigma}(ad)$; then we can rewrite the above expression
as follows

$\displaystyle\delta\ln\big{(}\Upsilon_{(1,2)}\big{)}=\oint_{\Sigma}\bigg{[}\delta\ln\left(\frac{b}{z^{K}}\right)\frac{\mathrm{d}}{\mathrm{d}z}\ln\left(\frac{ad}{z^{L}}\right)-\delta\ln\left(\frac{ad}{z^{L}}\right)\frac{\mathrm{d}}{\mathrm{d}z}\ln\left(\frac{b}{z^{K}}\right)$
$\displaystyle\hphantom{\delta\ln\big{(}\Upsilon_{(1,2)}\big{)}=\oint_{\Sigma}\bigg{[}}{}-\delta\ln\left(\frac{ad}{z^{L}}\right)\frac{K}{z}+\delta\ln\left(\frac{b}{z^{K}}\right)\frac{L}{z}\bigg{]}\frac{\mathrm{d}z}{2i\pi},$

and after integration by parts (which is now possible since all functions
involved are single-valued in $N(\Sigma)$)

$\displaystyle\delta\ln\left(\Upsilon_{(1,2)}\right)=\delta\oint_{\Sigma}\left[\ln\left(\frac{b}{z^{K}}\right)\frac{\mathrm{d}}{\mathrm{d}z}\ln\left(\frac{ad}{z^{L}}\right)-\ln\left(\frac{ad}{z^{L}}\right)\frac{K}{z}+\ln\left(\frac{b}{z^{K}}\right)\frac{L}{z}\right]\frac{\mathrm{d}z}{2i\pi},$

so that the transition function admits the explicit expression

$\displaystyle\Upsilon_{(1,2)}=\exp\left[\oint_{\Sigma}\left[\ln\left(\frac{b}{z^{K}}\right)\frac{\mathrm{d}}{\mathrm{d}z}\ln\left(\frac{ad}{z^{L}}\right)-\ln\left(\frac{ad}{z^{L}}\right)\frac{K}{z}+\ln\left(\frac{b}{z^{K}}\right)\frac{L}{z}\right]\frac{\mathrm{d}z}{2i\pi}\right].$

##  4 Conclusion

We conclude this short note with a few comments.

First of all the choice of $\Sigma$ as a single closed contour is not
necessary; we can have a disjoint union of closed contours or an unbounded
contour as long as $M(z)$ converges to the identity sufficiently fast as
$|z|\to\infty$ or near the endpoints. The considerations extend with trivial
modifications. Our approach is similar in spirit to the approach used in [9]
to express the tau function of a general isomonodromic system with Fuchsian
singularities in terms of an appropriate Fredholm determinant.

Much less clear to the writer is how to handle the case where $\Sigma$
contains intersections; in this case we should stipulate a local “no-
monodromy” condition at the intersection points as explained in [2, 3]. The
obstacle is not the issue of factorization but the fact that the resulting RHP
of the IIKS type leads to an operator $K$ (2.3) which is not of trace-class
(and not even Hilbert–Schmidt, which would be sufficient in order to construct
a Hilbert–Carleman determinant).

Nonetheless, the two form $\delta\omega_{M}$ defines a line bundle as
explained and therefore it is possible to compute the “transition functions”
of the line-bundle; this is precisely what is accomplished (in different
setting and in different terminology) in recent works [11, 12] and the
transition functions can be expressed in terms of explicit expressions,
analogously to what we have shown here in this general but simplified setting.

### Acknowledgements

The author wishes to thank Oleg Lisovyy for asking a very pertinent question
on the representation of the Malgrange form in terms of Fredholm determinants.
Part of the thinking was done during the author’s stay at the “Centro di
Ricerca Matematica Ennio de Giorgi” at the Scuola Normale Superiore in Pisa,
workshop on “Asymptotic and computational aspects of complex differential
equations” organized by G. Filipuk, D. Guzzetti and S. Michalik. The author
wishes to thank the organizers and the Institute for providing an opportunity
of fruitful exchange.

## References

  * [1] Basor E.L., Widom H., On a Toeplitz determinant identity of Borodin and Okounkov, Integral Equations Operator Theory 37 (2000), 397–401, math.FA/9909010.
  * [2] Bertola M., The dependence on the monodromy data of the isomonodromic tau function, Comm. Math. Phys. 294 (2010), 539–579, arXiv:0902.4716.
  * [3] Bertola M., Corrigendum: The dependence on the monodromy data of the isomonodromic tau function, arXiv:1601.04790.
  * [4] Bertola M., Cafasso M., The transition between the gap probabilities from the Pearcey to the Airy process – a Riemann–Hilbert approach, Int. Math. Res. Not. 2012 (2012), 1519–1568, arXiv:1005.4083.
  * [5] Borodin A., Okounkov A., A Fredholm determinant formula for Toeplitz determinants, Integral Equations Operator Theory 37 (2000), 386–396, math.CA/9907165.
  * [6] Chihara T.S., An introduction to orthogonal polynomials, Mathematics and its Applications, Vol. 13, Gordon and Breach Science Publishers, New York – London – Paris, 1978.
  * [7] Clancey K.F., Gohberg I., Factorization of matrix functions and singular integral operators, Operator Theory: Advances and Applications, Vol. 3, Birkhäuser Verlag, Basel – Boston, Mass., 1981.
  * [8] Gantmacher F.R., The theory of matrices. Vols. 1, 2, Chelsea Publishing Co., New York, 1959.
  * [9] Gavrylenko P., Lisovyy O., Fredholm determinant and Nekrasov sum representations of isomonodromic tau functions, arXiv:1608.00958.
  * [10] Harnad J., Its A.R., Integrable Fredholm operators and dual isomonodromic deformations, Comm. Math. Phys. 226 (2002), 497–530, solv-int/9706002.
  * [11] Its A., Lisovyy O., Prokhorov A., Monodromy dependence and connection formulae for isomonodromic tau functions, arXiv:1604.03082.
  * [12] Its A., Lisovyy O., Tykhyy Yu., Connection problem for the sine-Gordon/Painlevé III tau function and irregular conformal blocks, Int. Math. Res. Not. 2015 (2015), 8903–8924, arXiv:1403.1235.
  * [13] Its A.R., Izergin A.G., Korepin V.E., Slavnov N.A., Differential equations for quantum correlation functions, Internat. J. Modern Phys. B 4 (1990), 1003–1037.
  * [14] Malgrange B., Sur les déformations isomonodromiques. I. Singularités régulières, in Mathematics and Physics (Paris, 1979/1982), Progr. Math., Vol. 37, Birkhäuser Boston, Boston, MA, 1983, 401–426.
  * [15] Malgrange B., Déformations isomonodromiques, forme de Liouville, fonction $\tau$, Ann. Inst. Fourier (Grenoble) 54 (2004), 1371–1392.


# Simulating the nanometric track-structure of carbon ion beams in liquid
water at energies relevant for hadrontherapy

Pablo de Vera1 Stefano Simonucci2 Paolo E. Trevisanutto3,4 Isabel Abril5
Maurizio Dapor3,4 Simone Taioli3,4,6 and Rafael Garcia-Molina1 1Departamento
de Física - Centro de Investigación en Óptica y Nanofísica (CIOyN),
Universidad de Murcia, 30100 Murcia, Spain  2School of Science and Technology,
University of Camerino, 62032 Camerino, Italy, and INFN, Sezione di Perugia,
06123 Perugia, Italy  3European Centre for Theoretical Studies in Nuclear
Physics and Related Areas (ECT*-FBK), 38123 Trento, Italy  4Trento Institute
for Fundamental Physics and Applications, 38123 Trento, Italy  5Departament de
Física Aplicada, Universitat d’Alacant, 03080 Alacant, Spain  6Peter the Great
St. Petersburg Polytechnic University, Russia<EMAIL_ADDRESS>
###### Abstract

The nanometric track-structure of energetic ion beams in biological media
determines the direct physical damage to living cells, which is one of the
main responsibles of their killing or inactivation during radiotherapy
treatments or under cosmic radiation bombardment. In the present work,
detailed track-structure Monte Carlo simulations, performed with the code SEED
(Secondary Electron Energy Deposition), are presented for carbon ions in a
wide energy range in liquid water. Liquid water is the main constituent of
biological tissues, and carbon ions are one of the most promising projectiles
currently available for ion beam cancer therapy. The simulations are based on
accurate cross sections for the different elastic and inelastic events
determining the interaction of charged particles with condensed-phase
materials. The latter are derived from the ab initio calculation of the
electronic excitation spectrum of liquid water by means of time-dependent
density functional theory (TDDFT), which is then used within the dielectric
formalism to obtain inelastic electronic cross sections for both carbon ions
and secondary electrons. Both the ionisation cross sections of water by carbon
ions and the excitation and ionisation cross sections for electron impact are
obtained in very good agreement with known experimental data. The elastic
scattering cross sections for electrons in condensed-phase water are also
obtained from ab initio calculations by solving the Dirac-Hartree-Fock
equation. The detailed simulations fed with reliable cross sections allow to
assess the contribution of different physical mechanisms (electronic
excitation, ionisation and dissociative electron attachment –DEA–) to the
carbon ion-induced direct biodamage.

##  1 Introduction

Among the many applications of ion beams interaction with materials, the
irradiation of living tissues with energetic carbon ions is attracting much
attention for its benefits for treating deep-seated tumours in the technique
known as hadrontherapy [1, 2]. Due to the characteristic track-structure of
carbon ions (as compared to that of photon beams used in conventional
radiotherapy) [3], they produce large clusters of inelastic events within a
few turns of the sensitive DNA molecules, thus presenting a large relative
biological effectiveness or cell-killing ability.

The development of accurate biophysical models requires the detailed
simulation of carbon ions track-structure in liquid water (the main
constituent of biological tissue) in a wide energy range, from the MeV/u
energies of the primary therapeutic beams down to the keV/u energies at which
the stopping power and biological effects are largest. For this purpose,
probability distributions (cross sections) of electronic excitations by both
carbon ions and their secondary electrons (especially for energies below
50–100 eV, as these electrons are produced in large numbers by ion-impact
ionisation) need to be accurately known, in order to be used within Monte
Carlo simulation codes for the transport of radiation.

Reliable inelastic electronic cross sections can be obtained by means of the
dielectric formalism [4, 5] for both ions [6, 7] and electrons [8, 9] in
liquid water if its electronic excitation spectrum is accurately known. The
latter was determined experimentally for relatively wide energy and momentum
transfer ranges [10, 11], although theoretical models are useful to confirm
the experiments, interpret them, as well as to fill the gap where experiments
are missing. Phenomenological approaches such as the Mermin Energy Loss
Function – Generalised Oscillator Strengths (MELF-GOS) method have been shown
to provide reliable results [12, 13]. However, more recently, time-dependent
density functional theory (TDDFT) methods have demonstrated their ability to
achieve higher accuracies from purely ab initio calculations [14]. Still,
TDDFT is too expensive to be applied for arbitrarily large energies or
momenta, making analytical approaches such as MELF-GOS necessary to complement
the ab initio calculations [15].

In this contribution, we show how the TDDFT electronic excitation spectrum can
excellently reproduce the available experimental data up to 100 eV [10, 11]
and, from it and a high energy Aand high momenta extension by means of the
MELF-GOS method, we calculate differential and total cross sections for
electronic excitation and ionisation in very good agreement with the most
recent experimental information. These cross sections are then used within the
Monte Carlo code SEED (Secondary Electron Energy Deposition) [14, 16, 17], in
order to assess the clustering of direct damaging events produced by carbon
ions, in a wide energy range, in nanocylinders of volume equivalent to two DNA
convolutions located at different distances from the ions’ paths. From the
simulations, the role of different physical mechanisms to the clustered
nanometric biodamage can be assessed. The Monte Carlo simulations can also be
checked against experimental nanodosimetric distributions for ionisation
clusters, this interaction being the one which can be most easily measured in
the laboratory [18].

##  2 Excitation spectrum of liquid water

The energy loss function (ELF) of liquid water, Im$[-1/\epsilon(k,\omega)]$,
gives the excitation spectrum over the energy $\hbar\omega$ and momentum
$\hbar k$ transfers surface, $\epsilon(k,\omega)$ being the complex dielectric
function. The ELF is obtained here by TDDFT, using the Lanczos chains
algorithm (LCA) implemented in the turboEELS code [19] and the adiabatic PBE
(APBE) kernel [20], over the energy range $0\leq E\leq 100\mbox{ eV}$ and for
momentum transfers $0\leq\hbar k\leq 2.5\mbox{ a.u.}$, with a resolution of
0.25 a.u. These calculations were done for a particular configuration of a
periodic box formed by 32 water molecules and with liquid water density,
electronically relaxed via DFT calculations in Quantum Espresso [21] using the
PBE-GGA functional [20] and Troullier-Martins (TM) norm-conserving
pseudopotentials. Such a box was carved out from a supercell of several
thousands molecules optimised by means of molecular dynamics simulations using
the empirical TIP3P force-field [22] implemented in LAMMPS [23]. Figure 1
depicts by solid lines the TDDFT results for the ELF of liquid water for two
exemplary momentum transfers of (a) 0 a.u. and (b) 2.11 a.u., together with
the available experimental data (symbols) [10, 11]. The ab initio calculations
yield excellent results, reproducing very well the shape and intensity of the
different peaks corresponding to the excitations of the valence electrons of
liquid water.

Figure 1: Energy loss function of liquid water as a function of the energy
transfer $E=\hbar\omega$ for momentum transfers (a) $\hbar k=0$ a.u. and (b)
$\hbar k=2.11$ a.u.. Symbols correspond to experimental determinations for
liquid water [10, 11], while lines depict present TDDFT calculations.

##  3 Ionisation cross sections for carbon ions in liquid water

The space and energy characteristics of the ion track-structure are determined
by the angular and energy distributions of the secondary electrons which are
put in motion by the primary beam. Through the dielectric formalism, the ELF
over the entire ($k$, $\omega$) plane allows the calculation of the angular
and energy spectra of the electrons produced by carbon ions in liquid water,
following the model presented in Refs. [6, 7]. In short, the excitation and
ionisation of the valence shell electrons are distinguished by defining a mean
binding energy for these electrons, which for liquid water has a value of
13.71 eV [9]. The $\omega$-dependence of the ELF gives the energy spectrum of
secondary electrons [6], while the $k$-dependence, through a physically
motivated relation, provides their angular dependence [7]. The calculated
angular distributions of secondary electrons produced by 6 MeV/u C6+ ions in
liquid water are shown in Fig. 2, for several kinetic energies of the
secondary electrons $W$ (i.e., the doubly differential cross sections, DDCS),
together with the experimental measurements available for water molecules
[24].

Figure 2: Angular distributions of secondary electrons ejected from water at
different kinetic energies $W$ for the impact of 6 MeV C6+ ions. Symbols
correspond to experimental data for the water molecule [24], while lines are
calculations for liquid water based on the ELF obtained from TDDFT.

Despite the expected phase effects, in general the calculations reproduce very
well the experimental points recorded in water vapour in a wide range of
energies and angles, confirming the reliability of the model. For the lower
energies, the calculations predict a forward peaked distribution in opposition
to the experiments, a difference which may be due to phase effects. The figure
exemplifies how low energy electrons are produced in much larger quantity,
while higher energy electrons tend to be ejected at the angle defined by the
binary encounter peak.

##  4 Cross sections for secondary electrons in liquid water

The dielectric formalism also allows the calculation of electron inelastic
scattering cross sections for liquid water. For that purpose, the model needs
to be amended for indistinguishability and exchange between the primary and
excited electrons, both for electronic excitation and ionisation. Furthermore,
electronic structure effects and low-energy corrections need to be considered
in order to describe the energy-loss of secondary electrons [9].

Figure 3: Total cross sections for electron impact in water. (a) Inelastic
cross sections: Solid and dashed lines correspond to present calculations for
electronic excitation and ionisation using the ELF obtained from TDDFT, the
dotted line corresponds to the recommended DEA cross section [25], open
symbols depict experimental data for ionisation (down triangles [26], circles
[27], squares [28]) and full symbols for excitation (diamonds [29], triangles
[30], stars [31]). (b) Elastic cross section: the solid lines corresponds to
the Dirac-Hartree-Fock calculation for a cluster of six water molecules, the
dotted line to the Mott calculation for the water molecule, while symbols
correspond to experimental data for the water molecule (triangles [25],
squares [32], circles [33]).

Figure 3(a) represents by solid (dashed) lines the electronic excitation
(ionisation) cross sections calculated for liquid water. Symbols depict a
compilation of experimental data which, unfortunately, are only available for
water vapour. For ionisation, the calculations reproduce very well the
experimental points [26, 27, 28] in the entire energy range, particularly the
measurements from Refs. [26, 28]. Those from Ref. [27] are slightly larger, in
line with expected phase effects. For excitation, the experimental data
available for particular channels [29, 30, 31] are scaled in order to estimate
the total cross section [9], so as the data reproduces rather well the
calculations. The cross section for the other major inelastic channel, namely
dissociative electron attachment (DEA), is more difficult to calculate, and
here is approximated from recommendations based on experimental data [25],
shown in the figure by a dotted line.

As for the elastic scattering cross sections, they can be calculated ab initio
for water in the condensed phase. The Dirac-Hartree-Fock equation is solved by
using a projected-potential approach based on Gaussian functions [34, 35, 36]
for a cluster of six water molecules [14], in an effort to approximate the
liquid environment. The elastic cross sections are represented in Fig. 3(b).
Symbols correspond to experimental data for water vapour [25, 32, 33], while
the dotted line is that calculated for the water molecule within the Mott
theory. The solid line depicts the Dirac-Hartree-Fock calculation for a
cluster of six water molecules. While the calculation for the single water
molecule reproduces fairly well the experiments, the results for the water
cluster shows a clear deviation, which manifest the effect of the condensed
phase on the elastic scattering.

##  5 Track-structure simulations

The cross sections described above are the information needed to simulate the
secondary electron emission and transport by means of the Monte Carlo code
SEED [14, 16, 17], together with the electron-phonon and electron-polaron
interactions, which are included by means of the Fröhlich and Ganachaud-
Mokrani models [14, 16].

The program SEED is used to generate secondary electrons in liquid water along
carbon ion paths of different energies, the latter passing by at different
impact parameters from the centre of a sensitive nanocylinder, perpendicular
to the path, with dimensions similar to two convolutions (or 20 base pairs) of
a DNA molecule (2.3 nm diameter and 6.8 nm height). The collision geometry is
depicted in the inset of Fig. 4(c). Such a DNA-like target is commonly used in
radiobiological studies.

Figure 4: Clusters of inelastic events produced by carbon ions of different
energies passing by some impact parameter from the centre of a cylinder of
radius 2.3 nm and height 6.8 nm, similar to two DNA convolutions or 20 base
pair. The inset in the right hand panel depicts the collision geometry. (a)
Average size $M_{1}$ of the damage clusters produced by 0.2 MeV/u carbon ions,
as a function of the impact parameter, by combinations of particular inelastic
channels: ionisation, excitation and DEA. Open circles show the results for
excitation clusters when only 40% of them produce damage. Open squares
connected by lines represent the total size of the damage cluster when only
40% of excitations produce damage. (b) The same but for 83.26 MeV/u carbon
ions. (c) Cumulative probabilities $F_{n}^{\rm ionis}$ to produce clusters of
$n$ or more ionisations, as a function of the ionisation average cluster size
$M_{1}^{\rm ionis}$. Here lines represent the results of simulations, while
symbols correspond to experimental nanodosimetric data [18].

##  6 Results and discussion

The Monte Carlo simulations allow calculating the average size of clusters of
damaging events (ionisations, dissociative excitations, DEA) produced by
carbon ions in the DNA-like volumes placed at different distances from the
ions’ paths. The impact parameters analysed range from 1 to 100 nm, and the
carbon ion energies have been varied between 0.2 MeV/u and 83.26 MeV/u. While
initial energies of hundreds of MeV/u are typical in clinical carbon ion beams
(so they present ranges of the order of tens of centimetres in tissue), the
ions progressively lose their energy while traversing the body, going down to
energies of around hundreds of keV/u around the Bragg peak region. Thus, 0.2
MeV/u corresponds to a typical energy of carbon ions around the Bragg peak
(where their biological effects are most severe), while 83.26 MeV/u
corresponds to a situation closer to the plateau region of the depth-dose
curve. Such energetic carbon ions can also be found in cosmic radiation.

Figures 4(a) and (b) represent the average size of damage clusters, $M_{1}$,
as a function of the impact parameter, for these two limiting energies of 0.2
and 83.26 MeV/u. The clusters due to different physical mechanisms
(ionisations, excitations and DEA) are represented by different symbols. It is
known that only around 40% of the electronic excitations are able to produce
water molecule dissociation [37], so we only considered this proportion of
“dissociative excitations” as damaging events. Open circles depict the average
size of clusters produced by these dissociative excitations, while open
squares connected by lines are the average sizes of the clusters due to the
combination of ionisations, DEA and dissociative excitations. Error bars
represent the statistical uncertainties obtained from a large number of
simulations [14]. For 0.2 MeV/u, the average cluster size drops to zero for
impact parameters larger than 20 nm, as secondary electrons are not able to
reach large distances. For 83.26 MeV/u, the clusters are non negligible for
longer distances, due to the production of more energetic delta electrons.
Still, the average cluster sizes are much larger for 0.2 MeV/u than for 83.26
MeV/u, remarking the much larger relative biological effectiveness of carbon
ions around the Bragg peak region. Figures 4(a) and (b) also remarks the fact
that ionisations are the main events contributing to the size of the damage
clusters, followed by dissociative excitations, and, by far, by DEA. In fact,
ionisations make up for around 75% of the size of the clusters, while
dissociative excitations contribute around 20% and DEA only around 5% [14].

Experimental nanodosimetry allows measuring the clustering of ionisation
events in macroscopic gas-phase detectors having nanometric equivalent volumes
of liquid water. Figure 4(c) shows by symbols the distributions of the
cumulative probabilities $F_{n}^{\rm ionis}$ to produce clusters of size equal
or larger than $n$ (with $n=1$, 2 or 3) as a function of the average
ionisation cluster size $M_{1}^{\rm ionis}$ coming from a compilation of
experimental data [18]. It is worth to mention that such a curve is universal,
and independent of the dimensions and particularities of the nanodosimeters.
Calculated results coming from different ion energies and impact parameters,
depicted by lines, reproduce experiments in a wide range of conditions, what
validates the performance of the simulation code.

##  7 Summary and conclusions

As the direct damage of biological materials is determined by the clustering
of inelastic events in volumes similar to the dimensions of a few convolutions
of DNA, it is important to develop simulation tools capable to predict such
clustered damage under irradiation. Given the current interest in applying
carbon ion beams for cancer radiotherapy due to their increased relative
biological effectiveness, and considering that liquid water is the main
constituent of living tissue, we developed Monte Carlo simulations to study
the clustered direct damage produced by carbon ion beams, in a wide range of
energies relevant for hadrontherapy, in nanovolumes of liquid water.
Particular attention was paid to provide accurate ab initio descriptions of
the electronic excitation spectrum of liquid water as well as the elastic
collisions, in order to feed the Monte Carlo code with reliable cross
sections, so relevant to describe the propagation of the very low energy
secondary electrons.

Time-dependent density functional theory (TDDFT) calculations yielded the ELF
of liquid water in excellent agreement with the available experimental data.
The ELF was extended to large energy and large momentum transfers by means of
the MELF-GOS method and then used to obtain electronic cross sections within
the dielectric formalism. The energy and angular distributions of secondary
electrons produced by carbon ions reproduced well the known experiments. The
same finding applies to the electronic excitation and ionisation cross
sections for electrons when appropriate low energy corrections are
implemented. Elastic scattering cross sections for liquid water were estimated
by means of calculations for clusters of six water molecules by solving the
Dirac-Hartree-Fock equation.

Simulations with the SEED code allowed to determine the clustering of
inelastic events in cylinders of volume similar to two DNA convolutions,
showing how for low carbon ion energies the average size of clusters can be
hundreds of times larger than those for high energies far away from the Bragg
peak region. Such clusters are mainly composed, independently of the carbon
ion energy, by ionising events (up to 75–80%), while DEA has a minor
contribution. Moreover, the simulated cumulative distributions of ionisation
clusters reproduce rather well the available measurements, ionisations being
the inelastic events more easily accessible to experimental nanodosimetry.

This project received funding from the European Union’s Horizon 2020 Research
and Innovation programme under the Marie Sklodowska-Curie grant agreement no.
840752. M.D. and S.T. acknowledge the Bruno Kessler Foundation and the
National Institute of Nuclear Physics for unlimited access to their computing
facilities, and the Caritro Foundation for the grant High-Z ceramic oxide
nanosystems for mediated proton cancer therapy. This work was also supported
by the Spanish Ministerio de Ciencia e Innovación and the European Regional
Development Fund (Project PGC2018-096788-B-I00); the Fundación Séneca-Agencia
de Ciencia y Tecnología de la Región de Murcia (Project 19907/GERM/15); and
the Ministry of Science and Higher Education of the Russian Federation as part
of World-class Research Center program: Advanced Digital Technologies
(contract No. 075-15-2020-934 dated 17.11.2020).

## References

## References

  * [1] Schardt D, Elsässer T and Schulz-Ertner D 2010 Rev. Mod. Phys. 82 383–425
  * [2] Ebner D K and Kamada T 2016 Front. Oncol. 6 140
  * [3] Nikjoo H, Uehara S and Emfietzoglou D 2012  (CRC Press: Boca Raton, FL, USA)
  * [4] Lindhard J 1954 Det Kongelige Danske Videnskabernes Selskab. Matematisk-fysiske Meddelelser 28(8) 1–57
  * [5] Ritchie R H 1959 Phys. Rev. 114 644–654
  * [6] de Vera P, Garcia-Molina R, Abril I and Solov’yov A V 2013 Phys. Rev. Lett. 110 148104
  * [7] de Vera P, Garcia-Molina R and Abril I 2015 Phys. Rev. Lett. 114 018101
  * [8] de Vera P and Garcia-Molina R 2019 J. Phys. Chem. C 123 2075–2083
  * [9] de Vera P, Abril I and Garcia-Molina R 2021 Phys. Chem. Chem. Phys. 23 5079–5095
  * [10] Watanabe N, Hayashi H and Udagawa Y 1997 Bull. Chem. Soc. Jpn. 70 719–726
  * [11] Hayashi H, Watanabe N, Udagawa Y and Kao C C 2000 Proc. Natl. Acad. Sci. USA 97 6264–6266
  * [12] Heredia-Avalos S, Garcia-Molina R, Fernández-Varea J M and Abril I 2005 Phys. Rev. A 72 052902
  * [13] Garcia-Molina R, Abril I, Heredia-Avalos S, Kyriakou I and Emfietzoglou D 2011 Phy. Med. Biol. 56 6475–6493
  * [14] Taioli S, Trevisanutto P E, de Vera P, Simonucci S, Abril I, Garcia-Molina R and Dapor M 2021 J. Phys. Chem. Lett. 12 487–493
  * [15] Pedrielli A, de Vera P, Trevisanutto P E, Pugno N M, Garcia-Molina R, Abril I, Taioli S and Dapor M 2021 Phys. Chem. Chem. Phys. 23 19173–19187
  * [16] Dapor M, Abril I, de Vera P and Garcia-Molina R 2017 Phys. Rev. B 96 064113
  * [17] Dapor M 2020 Transport of Energetic Electrons in Solids. Computer Simulation with Applications to Materials Analysis and Characterization 3rd ed (Cham, Switzerland: Springer International Publishing AG)
  * [18] Conte V, Selva A, Colautti P, Hilgers G, Rabus H, Bantsar A, Pietrzak M and Pszona S 2018 Radiat. Prot. Dosim. 180 150–156
  * [19] Timrov I, Vast N, Gebauer R and Baroni S 2015 Comput. Phys. Commun. 196 460–469
  * [20] Olsen T, Patrick C E, Bates J E, Ruzsinszky A and Thygesen K S 2019 npj Comput. Mater. 5 106
  * [21] Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, Ceresoli D, Chiarotti G L, Cococcioni M, Dabo I, Corso A D, de Gironcoli S, Fabris S, Fratesi G, Gebauer R, Gerstmann U, Gougoussis C, Kokalj A, Lazzeri M, Martin-Samos L, Marzari N, Mauri F, Mazzarello R, Paolini S, Pasquarello A, Paulatto L, Sbraccia C, Scandolo S, Sclauzero G, Seitsonen A P, Smogunov A, Umari P and Wentzcovitch R M 2009 J. Phys. Condens. Matter 21 395502
  * [22] MacKerell A D, Bashford D, Bellott M, Dunbrack R L, Evanseck J D, Field M J, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau F T K, Mattos C, Michnick S, Ngo T, Nguyen D T, Prodhom B, Reiher W E, Roux B, Schlenkrich M, Smith J C, Stote R, Straub J, Watanabe M, Wiórkiewicz-Kuczera J, Yin D and Karplus M 1998 J. Phys. Chem. B 102 3586–3616
  * [23] Senn M 2020 LAMMPS Molecular Dynamics Simulator https://lammps.sandia.gov/
  * [24] Dal Cappello C, Champion C, Boudrioua O, Lekadir H, Sato Y and Ohsawa D 2009 Nucl. Instrum. Methods Phys. Res. B 267 781–790
  * [25] Itikawa Y and Mason N 2005 J. Phys. Chem. Ref. Data 34 1–22
  * [26] Schutten J, de Heer F J, Moustafa H R, Boerboom A J H and Kistemaker J 1966 J. Chem. Phys. 44 3924–3928
  * [27] Bolorizadeh M A and Rudd M E 1986 Phys. Rev. A 33 882–887
  * [28] Bull J N, Lee J W and Vallance C 2014 Phys. Chem. Chem. Phys. 16 10743–10752
  * [29] Thorn P A, Brunger M J, Teubner P J, Diakomichalis N, Maddern T, Bolorizadeh M A, Newell W R, Kato H, Hoshino M, Tanaka H, Cho H and Kim Y K 2007 J. Chem. Phys. 126 064306
  * [30] Ralphs K, Serna G, Hargreaves L R, Khakoo M A, Winstead C and McKoy V 2013 J. Phys. B 46 125201
  * [31] Matsui M, Hoshino M, Kato H, Da Silva F F, Limão-Vieira P and Tanaka H 2016 Eur. Phys. J. D 70 77
  * [32] Katase A, Ishibashi K, Matsumoto Y, Sakae T, Maezono S, Murakami E, Watanabe K and Maki H 1986 J. Phys. B 19 2715–2734
  * [33] Cho H, Park Y S, Tanaka H and Buckman S J 2004 J. Phys. B 37 625–634
  * [34] Taioli S, Simonucci S and Dapor M 2009 Comput. Sci. Discov. 2 015002
  * [35] Taioli S, Simonucci S, Calliari L and Dapor M 2010 Phys. Rep. 493 237–319
  * [36] Morresi T, Taioli S and Simonucci S 2018 Adv. Theory Simul. 1 1870030
  * [37] Thorn P A 2008 Electronic State Excitations in the Water Molecule By Collisions with Low Energy Electrons Ph.D. thesis