Approximate Generalized Additive-Quadratic Functional Equations on Ternary
Banach Algebras

Sedigheh Jahedi1,∗ and Vahid Keshavarz1

1Department of Mathematics, Shiraz University of Technology
P. O. Box 71557-13876, Shiraz, Iran
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>

Abstract. In this paper, we introduce the concept of j-hom-derivation,
$j\in\\{1,2\\}$ and solve the new generalized additive-quadratic functional
equations in the sense of ternary Banach algebras. Moreover, using the fixed
point method, we prove its Hyers-Ulam stability.

Keywords: Hyers-Ulam stability, ternary Banach algebra, additive function,
quadratic function, fixed point theorem.
Mathematics Subject Classification (2010): 17A40, 39B52, 17B40, 47B47.
00footnotetext: *Corresponding author

##  1\. Introduction and preliminaries

A ternary Banach algebra $X$ with $\|.\|$ is a complex Banach algebra equipped
with a ternary product $(xyz)\rightarrow xyz$ of $X^{3}$ into $X$. This
product is $\mathbb{C}$-linear in the outer variable, conjugate
$\mathbb{C}-$linear in the middle variable associative in the sense that
$xy(zwv)=x(zyw)v=(xyz)wv$ and satisfies
$\|xyz\|\leq\|x\|.\|y\|.\|z\|,~{}\|xxx\|=\|x\|^{3}$ (see [24]). Ternary
structures and their extensions, known as n-ary algebras have many
applications in mathematical physics and photonics, such as the quark model
and Nambu mechanics [13, 17]. Today, many physical systems can be modeled as a
linear system. The principle of additivity has various applications in physics
especially in calculating the internal energy in thermodynamic and also the
meaning of the superposition principle.

The Hyers-Ulam stability problem which arises from Ulam’s question says that
for two given fixed functions $\varphi$ and $\psi,$ the functional equation
$\mathcal{E}_{1}(F)=\mathcal{E}_{2}(F)$ is stable if for a function $f$ for
which $d(\mathcal{E}_{1}(f),\mathcal{E}_{2}(f))\leq\varphi$ holds, there is a
function $g$ such that $\mathcal{E}_{1}(g)=\mathcal{E}_{2}(g)$ and
$d(f,g)\leq\psi.$ [9, 11, 21, 23]. In 1941 [11], Hyers solved the
approximately additive mappings on the setting of Banach spaces.

The functional equation $f(x+y)=f(x)+f(y)$ is an additive equation and its
solution is called an additive mapping.

######  Theorem 1.1.

[11] Let X and Y be Banach spaces. Assume that $f:X\to Y$ satisfies

$\|f(x+y)-f(x)-f(y)\|\leq\varepsilon$

for all $x,y,z\in X$ and some $\varepsilon\geq 0.$ Then there exists a unique
additive mapping $T:X\to Y$ such that

$\|f(x)-T(x)\|\leq\varepsilon$

for all $x\in X.$

First, T. Rassias [21] and Aoki [1] and then a number of authors extended this
result by considering the unbounded Cauchy differences in different spaces.
For example see [2, 3, 7, 8, 10, 12, 14, 18]. F. Skof in 1983 [22], proved the
stability problem of quadratic functional equation between normed and Banach
spaces. The functional equation $f(x+y)+f(x-y)=2f(x)+2f(y)$ is called
quadratic equation. Cholewa [4] showed that the Skof’s theorem is also true
for the mappings defined on abelian groups. Later, a lot of research appeared
with various generalization of quadratic functional equations [5, 15, 20].
Consider the generalized additive-quadratic functional equation

$\begin{split}3^{j}f(\frac{x+y+z}{3})&+f(x)+f(y)+(-1)^{j}f(z)-2^{j}f(\frac{x+y}{2})-2^{j}f(\frac{y+z}{2})-(-1)^{j}2^{j}f(\frac{x+z}{2})\\\
&=\rho[jf(x+y+z)+jf(x)-f(x+y)-f(x+z)-(j-1)f(y+z)]\end{split}$ (1.1)

where $\rho\neq 0,\pm 1$ is a complex number and $j\in\\{1,2\\}$. In this
paper, we solve (1.1) and show that for $j=1,$ a function which satisfies
(1.1) is additive and for $j=2,$ it is quadratic. We also prove its Hyers-Ulam
stability by using the fixed point method. To do this, we use the Diaz-
Margolis fixed point theorem [16].

######  Theorem 1.2.

[16] Let $(X,d)$ be a complete generalized metric space and let
$F:X\rightarrow X$ be a strictly contractive mapping with Lipschitz constant
$0<L<1$. Then for each given element $x\in X$, either

$d(F^{n}(x),F^{n+1}(x))=\infty$

for all nonnegative integers $n$ or there exists a positive integer $n_{0}$
such that

$(1)$ $d(F^{n}(x),F^{n+1}(x))<\infty,\qquad\forall n\geq n_{0}$;

$(2)$ the sequence $\\{F^{n}(x)\\}$ converges to a unique fixed point $y^{*}$
of $F$ in the set $Y=\\{y\in X\mid d(F^{n_{0}}x,y)<\infty\\}$;

$(4)$ $d(y,y^{*})\leq\frac{1}{1-L}d(y,F(y))$ for all $y\in Y$.

##  2\. Main Results

Extending the concepts of m-homomorphism, $1\leq m\leq 4,$ and hom-derivation
which has been introduced by Eshaghi Gordji et. al [6] and park et. al [19],
respectively, we have the following definitions in the case of ternary Banach
algebras. Throughout the paper, X is a ternary Banach algebra.

We say that a functional $f$ is a j-mapping, $j\in\\{1,2\\}$ if for $j=1$ it
is an additive and for $j=2,$ f is a quadratic mapping.

######  Definition 2.1.

A mapping $f:X\to X$ is a ternary j-homomorphism, $j\in\\{1,2\\},$ if $f$ is a
j-mapping and

$f(xyz)=f(x)f(y)f(z)$

for all $x,y,z\in X.$

######  Definition 2.2.

Let $h:X\to X$ be a ternary j-homomorphism, $j\in\\{1,2\\}.$ A mapping $D:X\to
X$ is called a ternary j-hom-derivation if $D$ is a j-mapping and satisfies

$D(xyz)=D(x)(h(y))^{j}(h(z))^{j}+(h(x))^{j}D(y)(h(z))^{j}+(h(x))^{j}(h(y))^{j}D(z)$

for all $x,y,z\in X$.

For $j\in\\{1,2\\}$ and a given mapping $f:X\to X$ consider the equation

$\begin{split}\mathcal{E}^{j}f(x,y,z)=3^{j}f(\frac{x+y+z}{3})&+f(x)+f(y)+(-1)^{j}f(z)-2^{j}f(\frac{x+y}{2})-2^{j}f(\frac{y+z}{2})-(-1)^{j}2^{j}f(\frac{x+z}{2})\\\
&-\rho[jf(x+y+z)+jf(x)-f(x+y)-f(x+z)-(j-1)f(y+z)].\end{split}$ (2.1)

In the following we give the solution of the functional equation (1.1).

######  Proposition 2.3.

Let $f:X\to X$ be a mapping satisfies $\mathcal{E}^{j}(f(x,y,z))=0$,
$j\in\\{1,2\\}.$ Then $f(0)=0$ and
$(i)~{}$ for $j=1$, f is an odd and additive mapping;
$(ii)~{}$ for $j=2$, f is a quadratic mapping if f is an even function.

###### Proof.

First of all, note that $\mathcal{E}^{j}f(0,0,0)=0$ implies that $f(0)=0.$
$(i)~{}$ Suppose $j=1$. Let $x=0$ and $y=0$ in (2.1).Then

$3f(\frac{z}{3})=f(z)$ (2.2)

Put $x=0$ and replacing $z$ by $y$ in (2.1) and using (2.2), we obtain
$f(2y)-2f(y)=\rho[f(2y)-2f(y)].$ Since $\rho\neq 0,\pm 1,$ then

$f(2y)=2f(y).$ (2.3)

Again put $x=0$ in (2.1), use (2.2) and (2.3) we get
$\rho[f(y+z)-f(y)-f(z)]=0.$ So

$f(y+z)=f(y)+f(z)$

i.e., $f$ is additive. To show that $f$ is an odd function put $z=-y$ in the
above equation, so $f(-y)=-f(y)$.
$(ii)~{}$ Suppose $j=2$ and $f$ is an even function satisfies
$\mathcal{E}^{2}f(x,y,z)=0.$ So from $\mathcal{E}^{2}f(x,y,-x)=0$ we have

$\displaystyle
9f(\frac{y}{3})+f(y)+2f(x)-4f(\frac{x+y}{2})-4f(\frac{y-x}{2})=\rho[2f(y)+2f(x)-f(x+y)-f(y-x)].$
(2.4)

Letting $y=0$ and $x=2x$ in (2.4), then

$f(2x)=4f(x).$ (2.5)

Again in (2.4) put $x=0$ and use (2.5). So

$9f(\frac{y}{3})=f(y).$ (2.6)

Finally, by applying (2.5) and (2.6) in $\mathcal{E}^{2}f(x,y,z)=0,$

$f(x+y)+f(x-y)=2f(x)+2f(y).$

This completes the proof. ∎

In the sequel, assume $j\in\\{1,2\\}$ and $\delta,\sigma:X^{3}\to[0,\infty)$
are two functions satisfy conditions

$\delta(\frac{x}{2},\frac{y}{2},\frac{z}{2})\leq\frac{k}{2^{j}}\delta(x,y,z),$
(2.7)
$\sigma(\frac{x}{2},\frac{y}{2},\frac{z}{2})\leq\frac{k}{2^{3j}}\sigma(x,y,z)$
(2.8)

for all $x,y,z\in X$ and some $0<k<1.$ Clearly,
$\delta(0,0,0)=\sigma(0,0,0)=0$ and by induction one can obtain that

$2^{nj}\delta(\frac{x}{2},\frac{y}{2},\frac{z}{2})\leq k^{n}\delta(x,y,z),$
${2^{3nj}}\sigma(\frac{x}{2},\frac{y}{2},\frac{z}{2})\leq k^{n}\sigma(x,y,z)$

for all $n\in\mathbb{N}.$ If $f:X\to X$ is a function such that
$\|\mathcal{E}^{j}f(x,y,z)\|\leq\delta(x,y,z)$ then we have $f(0)=0.$

To prove the following results, we need consider two cases. For $j=1$, suppose
that $f_{j}$ and $g_{j}$ are odd and in case $j=2$ assume $f_{j}$ and $g_{j}$
are even functions.

######  Theorem 2.4.

Let $j\in\\{1,2\\}$ and $f_{j}:X\to X$ be a function satisfies

$\|\mathcal{E}^{j}f_{j}(x,y,z)\|\leq\delta(x,y,z)$ (2.9)

where $\delta:X^{3}\to[0,\infty)$ fulfill (2.7). Then there exists a unique
ternary j-mapping $h_{j}:X\to X$ such that

$\|f_{1}(x)-h_{1}(x)\|\leq\frac{k}{2(1-k)}\delta(0,x,-x),$ (2.10)
$\|f_{2}(x)-h_{2}(x)\|\leq\frac{k}{2(1-k)}\delta(x,0,-x)$ (2.11)

###### Proof.

Let $\mathfrak{A}$ be the set of all functions $g:X\to X$ with $g(0)=0$.
Define the mapping $Q_{j}:\mathfrak{A}\to\mathfrak{A}$ by
$Q_{j}(g)=2^{j}g(\frac{x}{2})$ and for every $g,h\in\mathfrak{A}$ and $x\in X$
define

$d_{1}(g,h)=\inf\\{\alpha>0~{}:~{}\|g(x)-h(x)\|\leq\alpha\delta(0,x,-x)\\},$
$d_{2}(g,h)=\inf\\{\alpha>0~{}:~{}\|g(x)-h(x)\|\leq\alpha\delta(x,0,-x)\\}$

where $\inf\emptyset=+\infty.$ It is easy to show that for each
$j\in\\{1,2\\}$, $d_{j}$ is a generalized metric on $\mathfrak{A}$ and
$(\mathfrak{A},d_{j})$ is a complete generalized metric space. Let
$g,h\in\mathfrak{A}.$ Then

$\displaystyle\|Q_{1}(g)(x)-Q_{1}(h)(x)\|\leq
2d_{1}(g,h)\frac{k}{2}\delta(0,x,-x),$
$\displaystyle\|Q_{2}(g)(x)-Q_{2}(h)(x)\|\leq
2^{2}d_{2}(g,h)\frac{k}{2^{2}}\delta(x,0,-x)$

Then $d_{j}\Big{(}Q_{j}(g),Q_{j}(h)\Big{)}\leq kd_{j}(g,h),$ i.e. $Q_{j}$ is a
contraction mapping. In case $j=1,$ put $x=0$ and $z=-y$ in (2.9), then

$\|f_{1}(y)-2f_{1}(\frac{y}{2})\|\leq\frac{1}{2}\delta(0,y,-y).$

In case $j=2,$ by setting $y=0$ and $z=-x$ in (2.9) we have

$\|f_{2}(x)-2^{2}f_{2}(\frac{x}{2})\|\leq\frac{1}{2}\delta(x,0,-x).$

Above relations imply that $d_{j}(f_{j},Q_{j}(f_{j}))\leq\frac{1}{2}.$ Hence
by Theorem 1.2, there exist a positive integer $n_{0}$ and a unique fixed
point $h_{j}$ of $Q_{j}$ in set
$\Omega=\\{g\in\mathfrak{A}:~{}d(Q_{j}^{n_{0}}(f),g)<\infty\\}$ and
$\lim_{n\to\infty}Q_{j}^{n}(f_{j})(x)=h_{j}(x),$ . SO, for all $x\in X,$
$Q_{j}(h_{j})(x)=h_{j}(x)$ and
$\lim_{n\to\infty}{2^{jn}}f(\frac{x}{2^{n}})=h_{j}(x).$ Also we have
$d_{j}(f_{j},h_{j})\leq\frac{k}{2(1-k)}.$ This implies that

$\|f_{1}(x)-h_{1}(x)\|\leq\frac{k}{2(1-k)}\delta(0,x,-x),$
$\|f_{1}(x)-h_{1}(x)\|\leq\frac{k}{2(1-k)}\delta(x,0,-x).$

But for each $j\in\\{1,2\\},$

$\displaystyle\|\mathcal{E}^{j}h_{j}(x,y,z)\|$
$\displaystyle=\lim_{n\to\infty}{2^{jn}}\|\mathcal{E}^{j}f(\frac{x}{2^{n}},\frac{y}{2^{n}},\frac{z}{2^{n}})\|$
$\displaystyle\leq\lim_{n\to\infty}{2^{jn}}\delta(\frac{x}{2^{n}},\frac{y}{2^{n}},\frac{z}{2^{n}})$
$\displaystyle\leq\lim_{n\to\infty}k^{n}\delta(x,y,z)$ $\displaystyle=0.$

Hence by Proposition 2.3, for each $j\in\\{1,2\\},h_{j}$ is a j-mapping and
the proof is complete. ∎

Now, we are going to prove Hyers-Ulam stability of ternary j-hom-derivations
in ternary Banach algebras corresponding to the functional equation (1.1), by
using the fixed point method.

######  Theorem 2.5.

Assume $\delta,\sigma:X^{3}\to[0,\infty)$ are two functions which satisfy
conditions (2.7) and (2.8) for some constant $k\in(0,1).$ Suppose that $f_{j}$
and $g_{j}$ fulfill the following relations

$\|\mathcal{E}^{j}f(x,y,z)\|\leq\delta(x,y,z),$ (2.12)
$\|\mathcal{E}^{j}g(x,y,z)\|\leq\delta(x,y,z),$ (2.13)
$\|f_{j}(xyz)-f_{j}(x)f_{j}(y)f_{j}(z)\|\leq\sigma(x,y,z),$ (2.14)
$\|g_{j}(xyz)-g_{j}(x)(f_{j}(y))^{j}(f_{j}(z))^{j}-(f_{j}(x))^{j}g_{j}(y)(f_{j}(z))^{j}-(f_{j}(x))^{j}(f_{j}(y))^{j}g_{j}(z)\|\leq\sigma(x,y,z).$
(2.15)

for all $x,y,z\in X.$ Then there exist a unique ternary j-homomorphism
$H_{j}:X\to X$ and a unique ternary j-hom-derivation $D_{j}:X\to X$ such that

$\|f_{1}(x)-H_{1}(x)\|\leq\frac{k}{2(1-k)}~{}\delta(0,x,-x),$ (2.16)
$\|f_{2}(x)-H_{2}(x)\|\leq\frac{k}{2(1-k)}~{}\delta(x,0,-x),$ (2.17)
$\|g_{1}(x)-D_{1}(x)\|\leq\frac{k}{2(1-k)}~{}\delta(0,x,-x),$ (2.18)
$\|g_{2}(x)-D_{2}(x)\|\leq\frac{k}{2(1-k)}~{}\delta(x,0,-x)$ (2.19)

for all $x\in X$ and some $0<k<1.$

###### Proof.

Let $\mathfrak{A},d_{j}$ and $Q_{j}$, $j\in\\{1,2\\},$ be those as defined in
the proof of Theorem 2.4. Similar to the proof of Theorem 2.4, there exist
unique j-mappings $H_{j},D_{j}$ from X into X such that

$H_{j}(x)=\lim_{n\to\infty}Q_{j}^{n}(f_{j})(x)=\lim_{n\to\infty}{2^{jn}}f_{j}(\frac{x}{2^{n}}),$
(2.20)
$D_{j}(x)=\lim_{n\to\infty}Q_{j}^{n}(g_{j})(x)=\lim_{n\to\infty}{2^{jn}}g_{j}(\frac{x}{2^{n}})$
(2.21)

and satisfying (2.16), (2.17),(2.18) and(2.19) as desired. Mappings $H_{j}$,
$j=\\{1,2\\},$ are ternary j-homomorphism. In fact, by (2.14) and (2.20), we
have

$\displaystyle\|H_{j}(xyz)-H_{j}(x)H_{j}(y)H_{j}(z)\|$
$\displaystyle=\lim_{n\to\infty}{2^{3nj}}\|f_{j}(\frac{xyz}{2^{3n}})-f_{j}(\frac{x}{2^{n}})f_{j}(\frac{y}{2^{n}})f_{j}(\frac{z}{2^{n}})\|$
$\displaystyle\leq\lim_{n\to\infty}2^{3nj}\sigma(\frac{x}{2^{n}},\frac{y}{2^{n}},\frac{z}{2^{n}})$
$\displaystyle\leq\lim_{n\to\infty}k^{n}\sigma(x,y,z)$ $\displaystyle=0.$

Using (2.15) and (2.21) shows that $D_{j}$ is a ternary j-hom-derivation. In
fact,

$\displaystyle\|D_{j}$
$\displaystyle(xyz)-D_{j}(x)(H_{j}(y))^{j}(H_{j}(z))^{j}-(H_{j}(x))^{j}D_{j}(y)(H_{j}(z))^{j}-(H_{j}(x))^{j}(H_{j}(y))^{j}D_{j}(z)\|$
$\displaystyle=\lim_{n\to\infty}{2^{3nj}}\|g_{j}(\frac{xyz}{2^{3n}})-g_{j}(\frac{x}{2^{n}})(f_{j}(\frac{y}{2^{n}}))^{j}(f_{j}(\frac{z}{2^{n}}))^{j}-(f_{j}(\frac{x}{2^{n}}))^{j}g_{j}(\frac{y}{2^{n}})(f_{j}(\frac{z}{2^{n}}))^{j}-(f_{j}(\frac{x}{2^{n}}))^{j}(f_{j}(\frac{y}{2^{n}}))^{j}g_{j}(\frac{z}{2^{n}})\|$
$\displaystyle\leq\lim_{n\to\infty}{2^{3nj}}\sigma(\frac{x}{2^{n}},\frac{y}{2^{n}},\frac{z}{2^{n}})$
$\displaystyle\leq\lim_{n\to\infty}k^{n}\sigma(x,y,z)$ $\displaystyle=0.$

Now, the proof is complete. ∎

In Theorem 2.4 and Theorem 2.5, by taking $k=2^{1-r}$ and

$\delta(x,y,z)=\sigma(x,y,z)=s(\|x\|^{r}+\|y\|^{r}+\|z\|^{r})$

where $x,y,z\in X,~{}r\neq 1$ and s are nonnegative real numbers, we obtain
the following result.

######  Corollary 2.6.

Let $r\neq 1$ and s be two elements of $\mathbb{R}_{+}.$ Suppose that
$\delta(x,y,z)=\sigma(x,y,z)=s(\|x\|^{r}+\|y\|^{r}+\|z\|^{r})$. Assume
$f_{j},g_{j}:X\to X,~{}j\in\\{1,2\\}$ are functions satisfying in (2.12),
(2.13), (2.14) and (2.15). Then there exist ternary j-homomorphism $H_{j}$ and
ternary j-hom-derivation $D_{j}$ such that

$\left\\{\begin{array}[]{ll}\|f_{j}(x)-H_{j}(x)\|\leq\frac{2s}{2-2^{r}}\|x\|^{r}&\\\
{}{}{}{}{}&\hbox{for~{} $r>1$}\\\
\|g_{j}(x)-D_{j}(x)\|\leq\frac{2s}{2-2^{r}}\|x\|^{r}&\end{array}\right.$

and

$\left\\{\begin{array}[]{ll}\|f_{j}(x)-H_{j}(x)\|\leq\frac{2s}{2^{r}-2}\|x\|^{r}&\\\
{}{}{}{}{}&\hbox{for~{} $r<1$}\\\
\|g_{j}(x)-D_{j}(x)\|\leq\frac{2s}{2^{r}-2}\|x\|^{r}&\end{array}\right.$

Declarations of interest: none.

## References

  * [1] T. Aoki, On the stability of the linear transformation in Banach spaces. J. Math. Soc. Japan, 2, 64–66 (1950).
  * [2] G. Chen, D. Li, O.V. Gaans and S. V. Lunel, stability results for nonlinear functional differential equations using fixed point methods, Indagationes Mathematicae, 29 671–686 (2018).
  * [3] C. K. Choi, Hyperstability of logarithmic-type functional equation on restricted domains, Indagationes Mathematicae, 30 240–249 (2019).
  * [4] P. W. Cholewa, Remarks on the stability of functional equations, Aequations Mathematicae, 27, No. 1-2, 76-86 (1984).
  * [5] S. Czerwik, On the stability of the quadratic mapping in normed spaces, Abh. Math. Sem. Univ. Hamburg 62, 59–64 (1992).
  * [6] M. Eshaghi Gordji, Z. Alizadeh, H. Khodaei and C. Park, On approximate homomorphisms: a fixed point approach, Mathematical Sciences, 6 (2012), doi.org/10.1186/2251-7456-6-59.
  * [7] M. Eshaghi Gordji, H. Khodaei, Th. M. Rassias, Fixed points and generalized stability for quadratic and quartic mappings in ${C^{*}}$-algebras, J. Fixed Point Theory and Appl., 17, 703–715 (2018).
  * [8] I.-I. EL-Fassi, S. Kabbaj and A. Charifi, Hyperstability of Cauchy-Jensen functional equations, Indagationes Mathematicae, 27, 855–867 (2016).
  * [9] P. Găvruta, A generalization of the Hyers-Ulam-Rassias stability of approximately additive mappings, J. Math. Anal. Appl. 184, 431–436 (1994).
  * [10] I. Hwang and C. Park, Bihom derivations in Banach algebras, J. Fixed Point Theory Appl., (2019) doi.org/10.1007/s11784-019-0722-y.
  * [11] D. H. Hyers, On the stability of the linear functional equation, Proceedings of the National Academy of Sciences of the United States of America, 27, 222–224, (1941).
  * [12] S. Jahedi, V. Keshavarz, C. Park and S. Yun, Stability of ternary Jordan bi-derivations on C*-ternary algebras for bi-Jensen functional equation, J. Comput. Anal. Appl., 26, NO. 1, 140-145 (2019).
  * [13] Kerner R: Ternary algebraic structures and their applications in physics. Pierre et Marie Curie University, Paris; 2000.http://arxiv.org/list/math-ph/0011 Ternary algebraic structures and their applications in physics, Proc. BTLP, 23rd International Conference on Group Theoretical Methods in Physics, Dubna, Russia (2000).
  * [14] V. Keshavarz, S. Jahedi and M. Eshaghi Gordji, Ulam-Hyers stability of $C^{*}$-ternary 3-Jordan derivations, South East Asian Bulletin of Mathematics, preprint.
  * [15] H.-M. Kim and H.-Y. Shin, Refined stability of additive and quadratic functional equations in modular spaces, J Inequal Appl., 146, (2017) doi.org/10.1186/s13660-017-1422-z.
  * [16] B. Margolis, J.B. Diaz, A fixed point theorem of the alternative for contractions on the generalized complete metric space, Bull. Amer. Math. Soc., 126, 305–309 (1968).
  * [17] Y. Nambu, Generalized Hamiltonian mechanics, Phys. Rev., 7, 2405–-2412 (1973).
  * [18] I. Nikoufar, Jordan $(\theta,\phi)-$derivations on Hilbert $C^{*}$-Modules, Indagationes Mathematicae, 26, 421–430 (2015).
  * [19] C. Park, J. Rye Lee and X. Zhang, Additive s-functional inequality and hom-derivations in Banach algebras, J. Fixed Point Theory Appl., (2019) doi.org/10.1007/s11784-018-0652-0.
  * [20] C. Park, K. Tamilvanan, G. Balasubramanian, B. Noori and A. Najati, On a functional equation that has the quadratic-multiplicative property, Open Mathematics, 18, 837-–845 (2020).
  * [21] Th. M. Rassias, On the stability of the linear mapping in Banach spaces, Proc. Amer. Math. Soc., 72, 297–300 (1978).
  * [22] F. Skof, Propriet$\acute{a}$ locali e approssimazione di operatori, Rend. Sem. Mat. Fis. Milano 53, 113–-129 (1983).
  * [23] S. M. Ulam, Problems in Modern Mathematics, Chapter VI,Sscience ed. Wiley, New York, 1940.
  * [24] H. Zettl, A characterization of ternary rings of operators, Adv. Math., 48, 117–143 (1983).


# Deep Learning-Based Segmentation of Tumors in PET/CT Volumes: Benchmark of
Different Architectures and Training Strategies

Monika Górka 1,†, Daniel Jaworek 1,†, and Marek Wodzinski 1,2

1 AGH University of Krakow, Department of Measurement and Electronics
Krakow, Poland<EMAIL_ADDRESS>
2 University of Applied Sciences Western Switzerland (HES-SO Valais)
Information Systems Institute, Sierre, Switzerland
† These authors contributed equally to this work

###### Abstract

Cancer is one of the leading causes of death globally, and early diagnosis is
crucial for patient survival. Deep learning algorithms have great potential
for automatic cancer analysis. Artificial intelligence has achieved high
performance in recognizing and segmenting single lesions. However, diagnosing
multiple lesions remains a challenge. This study examines and compares various
neural network architectures and training strategies for automatically
segmentation of cancer lesions using PET/CT images from the head, neck, and
whole body. The authors analyzed datasets from the AutoPET and HECKTOR
challenges, exploring popular single-step segmentation architectures and
presenting a two-step approach. The results indicate that the V-Net and nnU-
Net models were the most effective for their respective datasets. The results
for the HECKTOR dataset ranged from 0.75 to 0.76 for the aggregated Dice
coefficient. Eliminating cancer-free cases from the AutoPET dataset was found
to improve the performance of most models. In the case of AutoPET data, the
average segmentation efficiency after training only on images containing
cancer lesions increased from 0.55 to 0.66 for the classic Dice coefficient
and from 0.65 to 0.73 for the aggregated Dice coefficient. The research
demonstrates the potential of artificial intelligence in precise oncological
diagnostics and may contribute to the development of more targeted and
effective cancer assessment techniques

Keywords: Deep Learning, Segmentation, Cancer, PET/CT, AutoPET, HECKTOR,
Benchmark

##  1 Introduction

###  1.1 Overview

Globally, cancer is acknowledged as one of the primary sources of premature
mortality, alongside cardiovascular diseases, and presents a formidable global
health challenge. It is predicted that by the end of the century, cancer will
overtake cardiovascular diseases, becoming the primary cause of premature
mortality in most nations across the world [9]. Early detection of cancer
lesions in patients is crucial for improving survival rates. The prognosis and
treatment options depend on the location and stage of the lesions [25, 10].
Computed tomography (CT) and positron emission tomography (PET) are playing a
key role in tumor analysis, providing valuable information about the location,
anatomy and stage of the tumor. Manual analysis and interpretation of medical
images is a time-consuming and error-prone process that may lead to a lack of
consistency in conclusions between different specialists [27, 29]. Therefore,
it is extremely important to continuously develop deep learning methods to
incorporate them into healthcare to accelerate and improve diagnostic
processes and analysis of biological data, which will increase the
effectiveness of detecting and understanding cancer mechanisms [34].
Currently, despite the high effectiveness of deep learning methods in
detecting and segmenting single lesions, the task of diagnosing multiple
lesions is still much more difficult due to low variability between lesions or
too large a number of lesions [19]. Another significant challenge is the
detection of small-volume tumors at an early stage of their development. In
this situation, the PET/CT imaging technique plays an invaluable role,
enabling the localization of even minimal cancer lesions. Moreover, the
integration of methods based on deep learning in the analysis of data obtained
from PET/CT imaging can significantly increase the efficiency of detection of
early-stage small-volume tumors. This, in turn, has a direct impact on
increasing the therapeutic potential and improving the life prognosis of
cancer patients.

###  1.2 Related Work

Many publications have explored the application of deep learning techniques in
automated cancer analysis. Our study focused on two issues: automatic analysis
of head and neck cancers and automatic analysis of whole-body cancers. Authors
of [35] reviewed the literature on deep learning techniques in multiomics
analysis and diagnosis of head and neck cancers, with a focus on convolutional
neural networks (CNNs). In [17] advanced methods of deep learning were
introduced by the authors, such as convolutional autoencoders, generative
adversarial networks (GANs), and transformers, to perform automatic analysis
of head and neck cancers for multiple tasks. In [20] the applications of Deep
Learning Radiomics (DLR) technology in predicting head and neck cancer were
described, and it was concluded that DLR has significant potential for
segmenting regions of interest (ROI), extracting and fusing features, and
constructing models using PET/CT volumes. Article [14] reviews the use of
artificial intelligence to enhance the radiotherapy workflow for head and neck
cancer, focusing on segmentation, planning, and delivery. Study [36] focuses
on the comparison of deep neural networks for automatic localization and
segmentation of head and neck cancers based on PET/CT images. A study [32]
aimed to assess the performance of PET and CT image fusion for gross tumor
volume segmentations of head and neck cancers utilizing conventional, deep
learning, and output-level voting-based fusions. In [39] an approach for
automatically segmenting multiple organs in three-dimensional computed
tomography (CT) images using deep learning was proposed. In study [8],
researchers examined the use of a cascade convolutional neural network
architecture to automatically identify and segment tumors in whole-body FDG-
PET images. Study [19] provided a review of deep learning algorithms that can
detect multifocal neoplastic cancer lesions in the human body, including
scenarios that involve whole-body imaging. The study [31] presents an
unconventional approach to precise initialization of organs-at-risk and
tumors. The authors proposed a framework implementing a cascade coarse-to-fine
segmentation with adaptive module for both small and large organs. The paper
[38] presents an overview of the latest deep learning-based multimodal tumor
segmentation methods. Public datasets, evaluation methods and multimodal data
processing are presented, as well as typical deep learning network structures,
techniques and methods of multimodal image fusion used in various tumor
segmentation tasks. Publication [15] summarizes the AutoPET Grand Challenge
2022 and describes the whole-body tumor segmentation methods proposed by the
participants. Publications [6, 5, 4] provide summaries of the HECKTOR
challenges from 2020, 2021, and 2020, respectively, along with a description
of the head and neck tumor segmentation methods employed by the participants.

###  1.3 Contribution

This study contributes to the field of diagnostic imaging in oncology by
analyzing, modifying, and comparing advanced deep learning algorithms for
cancer lesion segmentation from PET/CT images. The study evaluates the
performance of various neural network architectures commonly used in medical
image analysis. The careful comparison of these models enabled the
identification of the most effective architectures for specific diagnostic
applications, both for whole-body images and specific head and neck regions.
Furthermore, emphasis was placed on the two-stage segmentation methodology,
which, although less common in the literature, shows promising results
compared to traditional one-stage segmentation. Additionally, the study
examines the effects of excluding non-cancer cases from the training dataset,
which is a novel and counter-intuitive method for optimizing deep learning
models. Consequently, this research offers valuable insights into selecting
and configuring neural network models to enhance their diagnostic imaging
capabilities. These findings have significant implications for the development
of more advanced and accurate diagnostic tools in oncology. This could lead to
better clinical outcomes and increased survival rates for cancer patients.

##  2 Materials and Methods

###  2.1 Overview

To achieve satisfactory results in segmentation tasks, we focused on tools for
effectively extracting features from medical images, despite the limited
available data. The common small dataset problem arises from the challenge of
obtaining large, well-labeled, and high-quality datasets. The efficacy of deep
learning approaches is intricately connected to the size and quality of the
available dataset [35, 17]. Our research required tools capable of effectively
segmenting a variety of tumor types with different anatomical structures
located in different parts of the patient’s body. In our study, we used two
different datasets. The first one came from the HECKTOR (Head and Neck Tumor)
Grand Challenge, and the second one came from the AutoPET Grand Challenge held
at MICCAI (Medical Image Computing and Computer Assisted Intervention 2022).
We used the UNet, VNet, and UNETR models to perform segmentation on the
HECKTOR dataset, and the UNet, UNETR, and nnUnet models to perform
segmentation on the AutoPET dataset.

###  2.2 U-Net

U-Net’s fully connected layered encoder-decoder architecture has
revolutionized semantic segmentation. The skip connections between the encoder
and decoder are the essential component that enables the efficient
concatenation of feature maps. Symmetric upsampling and downsampling processes
significantly improve localization and segmentation accuracy. U-Net’s
flexibility, adaptability, and impressive performance for a variety of medical
imaging modalities with minimal training data have led to its popularity.
Consequently, this architecture has been modified for use in various
challenges. For this reason, we decided to use the U-Net architecture from the
MONAI library, which is a set of open-source, freely available collaborative
frameworks built for accelerating research and clinical collaboration in
medical imaging [19, 7, 24, 21, 12, 30]. To adapt the architecture to our
datasets we used the standard version of the U-Net and modified its parameters
without affecting the implementation of the model. Our modification included
two input channels corresponding to the two modalities used: CT and PET. The
network contained five layers with two convolutional units, with a kernel size
of 3x3x3, where each coding layer contained residual units with two
convolutional blocks and residual connections. The network was implemented
with a PReLU (Parametric Rectified Linear Unit) activation function, an
extension of ReLU, and a transposed convolution operation in the decoder,
enabling a gradual reduction in the number of channels and an increase in
spatial dimensions. For the AutoPET dataset, the network returns a three-
dimensional binary mask indicating the detected tumor areas, while for the
HECKTOR dataset, the network returns two binary masks - each mask
corresponding to one of the two types of tumors present in the images.

###  2.3 UNETR

Transformers, first introduced in 2017 for natural language processing, have
become increasingly important in medical image segmentation. Vision
Transformers (ViT) offer a critical method for dividing images into fragments,
which are processed by a transformer with a central multi-head auto-attention
(MSA) layer. Models like TransU-Net and UNETR, which use transformers, enable
more accurate image analysis by capturing long-range relationships and complex
contexts, which is important for precise segmentation of medical images. Our
study aimed to evaluate transformer effectiveness utilizing the UNETR model,
which combines the transformer architecture with U-Net elements. This model
divides the image into equal, non-overlapping fragments and projects them into
the embedding space using a linear layer. Next, this sequence is augmented
with positional embeddings and serves as input to the transformer model.
Tensors are obtained at each resolution level of this model and then connected
in a manner analogous to the standard U-net architecture. At the last stage, a
convolutional layer, alongside a softmax activation function, is used [12, 37,
33, 11]. In our study, we used the UNETR model from the MONAI library. Like
the UNet model, the network takes a two-channel input and returns one binary
mask in the case of the AutoPET set and two binary masks in the case of the
HECKTOR set.

###  2.4 V-Net

V-Net is a convolutional neural network specifically designed for processing
volumetric medical data. The network architecture was inspired by the U-Net
architecture and also consists of encoder and decoder parts that effectively
extract image features and details, enabling precise segmentation [23]. The
implementation we used from the MONAI library consists of 5 layers with
volumetric convolution, ELU (Exponential Linear Unit) activation function and
dropout in the last two layers. We applied the V-Net only to the HECKTOR
dataset. The network takes a two-channel image as input and returns two binary
masks.

###  2.5 nnU-Net

nnU-Net is currently one of the most prevalent tools for semantic segmentation
in various scientific fields. It has achieved top results in several open
challenges in medical imaging. nnU-net is based not only on the proper
preparation of the neural network, but also on the automatic configuration of
pre- and post-processing. Its decision-making process is based on a set of
parameters, rules and empirical decisions. The framework generates a ”dataset
fingerprint” by examining the provided training examples, which is then used
to develop appropriate configurations for datasets: 2d, 3d_fullres, 3d_lowres
and 3d_cascade_fullres. Due to the automatic parameter adjustment and
preprocessing, the aforementioned data preprocessing activities were not
required. During training of the nnU-Net datasets, the 5-fold cross-validation
method was employed. This method facilitated the evaluation of different
configurations and resulted in the selection of the most optimal segmentation
for a given problem [18, 1].

###  2.6 Datasets

We used two datasets in our study. The first one came from the HECKTOR (Head
and Neck Tumor) Grand Challenge, and the second one came from the AutoPET
Grand Challenge held at MICCAI (Medical Image Computing and Computer Assisted
Intervention 2022) [15, 4].

HECKTOR dataset includes 3D scans of the head and neck area from PET and CT.
For training cases, segmentation masks with marked tumor areas are
additionally included. The masks present two structures: primary Gross Tumor
Volumes (GTVp) and nodal Gross Tumor Volumes (GTVn). Data come from 9 centers
situated in Canada, Switzerland, France, and the United States. The entire
dataset contains 845 cases, divided into 524 training cases and 356 test cases
[4, 2]. Only data from the training set was used to conduct the experiments
due to the lack of access to segmentation masks for the test set.

The AutoPET dataset contains FDG-PET/CT images that cover the whole body from
the base of the skull to the mid-thigh. The dataset include 1014 images
obtained from 900 patients who took part in a prospective clinical trial
conducted between 2014 and 2018 at the University Hospital of Tübingen. Among
all the images, 501 were positive cases, indicating cancerous lesions that
displayed enhanced metabolic activity of the FDG factor. This group comprises
patients with non-small cell lung cancer (NSCLC), malignant lymphoma, or
melanoma. The remaining 513 images exhibit negative cases in which no
metabolically active cancer has been detected [16].

###  2.7 Experimental Setup

The first step was to prepare the data. Preprocessing of the HECKTOR dataset
included resampling to a common resolution and then cropping to the head and
neck area to the size of 192x192x192. Images were cropped in the sagittal and
coronal planes from the center of the image and in the transverse plane from
the top of the image. The AutoPET data has already been resampled to a common
resolution by the challenge organizers. The following preprocessing steps were
the same for both data sets. The CT image intensity was cut to the range
[-1024, 1024] Hounsfield Units, since this range covers most anatomical
structures of the human body, and then rescaled to [-1, 1]. The PET image
intensity was also rescaled to the range of [-1, 1]. The PET and CT images
were then superimposed to result in a single two-channel image. For the
HECKTOR dataset, each segmentation mask was divided into two binary masks. One
mask corresponded to the GTVp areas, and the other mask corresponded to the
GTVn areas. Cases with only one registered structure had one blank mask. In
the case of AutoPET data, the original masks were already binary. There was no
need to split them.

Data augmentations were utilized in all experiments to minimize overfitting.
Augmentation methods included random flipping, rotation, scaling, translation,
gaussian noise, and random gamma.

In the case of AutoPET data, tumors were located in different areas of the
images, and there were discrepancies in size across images. To ensure both
full resolution and high data quality while maintaining computational
optimality, we implemented a patch-based approach for the U-Net model. This
involved training the model on subsets of the full image, thus enabling
resolution to be maintained. During the experiments, we utilized a modified
Label Sampler that employed dynamic adjustment of a parameter responsible for
controlling class sampling probability in images [28]. Patch size was
standardized to 96x96x96, and 80 samples per image volume were taken. For the
UNETR architecture, we decided to adopt a standard training method for all
images due to hardware restrictions and the intricacy of the transformer
structure. As a result, all images were cropped to 192x192x224.

Another challenge with AutoPET data, as mentioned in Section 2.6, is that
approximately half of the images represent negative cases without cancerous
lesions. Consequently, the impact of these cases on the quality of
segmentation was evaluated through an experiment solely utilizing data
featuring positive controls. Non-cancerous images were removed from the
dataset. Comparative segmentation was then performed in the corresponding
configurations for U-Net and UNETR.

Figure 1: Diagram illustrating a implementation of two-step segmentation
process with three-channel aproach.

In order to improve the results of the U-Net model, an additional segmentation
step was performed, which included first identifying potential tumor areas.
Initially, a strategy was implemented to train the model on reduced dimension
and resolution data. The data underwent processing via the applied dynamic
resampling strategy, with input images being scaled to 192x192x192. As the
initial model, the U-Net architecture with the original parameter
configuration was applied. Following the model’s initial training,
segmentation mask predictions were generated across all data sets. The masks
were reinserted into the spatial context of their corresponding source images
and added as a third input channel to the U-Net model using the described
configuration.

All models and configurations described above were trained using the Adam
optimizer and Cosine Annealing with Warm Restarts learning rate scheduler,
which cyclically lowers the value of the learning rate and then resets it to
the initial value to prevent the model from getting stuck in the local minimum
gradient during training. The loss function used was DiceFocalLoss. This
function returns the weighted sum of Dice Loss and Focal Loss.

##  3 Results

The effectiveness of the models was assessed using commonly used metrics for
image segmentation tasks. The results were evaluated using the dice similarity
coefficient (DSC), which is defined by the following formula:

$DSC\left(A,B\right)=\frac{2\left|{A\cap
B}\right|}{\left|{A}\right|+\left|{B}\right|},$ (1)

where $A$ and $B$ are two sets whose similarity should be assessed. It was
decided to use the aggregated Dice coefficient as the second metric. In the
context of image segmentation, there may be two binary images representing
regions of interest. The aggregated Dice coefficient provides information on
the overall segmentation quality for all regions of interest, which allows for
a comprehensive assessment of the algorithm’s performance.

$DSC_{agg}=\frac{2\sum_{i}\left|{A_{i}\cap
B_{i}}\right|}{\sum_{i}\left|{A_{i}}\right|+\left|{B_{i}}\right|},$ (2)

where $A_{i}$ and $B_{i}$ are the ground truth and prediction for the image i
in the test set.

###  3.1 AutoPET

As part of our study, we conducted a comparative analysis of the effectiveness
of three different one-step segmentation methods. The experiment utilized
three models: U-Net architecture, UNETR, and nnU-Net, a dedicated framework
for image segmentation. The models were trained on two datasets. The first
dataset contained all images of the initial dataset, while the second dataset
was reduced by negative control cases in which there was no case of cancer.
The second part of the study was to perform two-step segmentation on the U-Net
model. The models’ effectiveness was evaluated using two indicators: the Dice
similarity coefficient and the aggregated Dice similarity coefficient. The
obtained results are presented in Table 1. The evaluation was performed on a
separate test set from the original dataset included images with cancerous
lesions. Exemplary visualizations of the obtained results are shown in Figure
2.

Table 1: Presentation of the results of various segmentation methods. |  | Full dataset | Tumor-only dataset
---|---|---|---
|  |  DSC ($\uparrow$) |  DSCagg ($\uparrow$) |  DSC ($\uparrow$) |  DSCagg ($\uparrow$)
|  U-Net |  $\mathrm{0.58}$ |  $\mathrm{0.64}$ |  0.74 |  $\mathrm{0.80}$
One-step |  nnU-Net |  0.69 |  0.80 |  0.74 |  0.84
|  UNETR |  $\mathrm{0.39}$ |  $\mathrm{0.51}$ |  $\mathrm{0.51}$ |  $\mathrm{0.54}$
Two-steps |  U-Net |  $\mathrm{0.60}$ |  $\mathrm{0.73}$ |  - |  -

During the statistical analysis, we focused on analyzing the significance of
differences between a specific model trained on the full dataset and the set
containing only cancer data. To determine the statistical significance of
differences in results for models trained on different datasets, a p-value was
calculated. The Wilcoxon test [22] was used to calculate the p-value with the
null hypothesis (H0) that the differences in results between the analyzed
models are not statistically significant at the 0.05 significance level.

The obtained p-values are presented in Table 2. P values below 0.00001 were
approximated to 0.

Table 2: Presentation of the p-values for pairs of all models. Full dataset |  Tumor-only dataset |  p-value
---|---|---
U-Net |  U-Net |  ~$0$
nnU-Net |  nnU-Net |  $\mathrm{0.0031}$
UNETR |  UNETR |  $\mathrm{0.0017}$
Figure 2: Visualization of the most accurate segmentation (DSC: 0.94) for
various nnU-Net configurations trained on the tumor-only dataset. The ground
truth tumor is marked in blue, while the model predictions are marked in red,
orange, green, and yellow.

###  3.2 HECKTOR

Three deep neural network architectures were compared on the HECKTOR dataset:
U-Net, UNETR, and V-Net. The performance of the models was evaluated on a test
set that was previously separated from the training set. Before the
evaluation, the generated predictions were put back into the full-dimensional
spatial context of the original corresponding CT image. Model performance was
assessed using the aggregated Dice similarity coefficient. The script provided
by the organizers of the HECKTOR challenge was used for evaluation [3]. The
results are shown in Table 3.

Table 3:  Segmentation results on the HECKTOR test set for models trained on cropped images. |  DSCaggGTVp ($\uparrow$) |  DSCaggGTVn ($\uparrow$) |  mean DSCagg($\uparrow$)
---|---|---|---
U-Net |  $\mathrm{0.79}$ |  $\mathrm{0.72}$ |  $\mathrm{0.76}$
V-Net |  $\mathrm{0.80}$ |  $\mathrm{0.72}$ |  $\mathrm{0.76}$
UNETR |  $\mathrm{0.79}$ |  $\mathrm{0.69}$ |  $\mathrm{0.74}$

To determine the statistical significance of the disparities in the test set
outcomes among the models, the p-value was calculated. The metric used to
evaluate the models (DSCagg) is a global metric, calculated for the entire
test set instead of for individual samples. Conventional statistical tests
require multiple observations, and global metrics cannot directly satisfy this
requirement. Therefore, bootstrap resampling [13, 26] was used. The test set
was repeatedly sampled with replacement to generate new test sets. DSCagg was
calculated for each drawn sample, allowing for the estimation of metric
distributions across models.

As with the AutoPET data, the p-value was calculated using the Wilcoxon test
with a null hypothesis that assumes no statistically significant differences
in results between the analyzed models at a significance level of 0.05. The
resultant p-values have been presented in Table 4.5. P-values below 0.00001
have been approximated to 0.

Table 4: Matrix displaying the p-values for pairs of segmentation models.
$\begin{array}[]{M{1.3cm}|M{1.3cm}|M{1.3cm}|M{1.3cm}}\@add@centering&\rotatebox[origin={c}]{90.0}{U-Net}\@add@centering&\rotatebox[origin={c}]{90.0}{V-Net}\@add@centering&\rotatebox[origin={c}]{90.0}{UNETR}\@add@centering\cr\hline\cr
U-Net\@add@centering&\text{$-$}\@add@centering&\textasciitilde{$0$}\@add@centering&\textasciitilde{$0$}\@add@centering\cr\hline\cr
V-Net\@add@centering&\textasciitilde{$0$}\@add@centering&\text{$-$}\@add@centering&\textasciitilde{$0$}\@add@centering\cr\hline\cr
UNETR\@add@centering&\textasciitilde{$0$}\@add@centering&\textasciitilde{$0$}\@add@centering&\text{$-$}\@add@centering\end{array}$

The models were evaluated using a test set extracted from the training set.
This is distinct from the data on which the results of the HECKTOR challenge
were evaluated. The results were submitted to the challenge to obtain the
results for the HECKTOR test set, in order to be able to compare the obtained
results with those of the challenge participants. Due to the limited number of
submissions, and to avoid bias towards the test set, it was decided to send
only the results of the best model. The results for the HECKTOR Grand
Challenge test set compared to the results of challenge participants are shown
in Table 5. Exemplary visualizations of the obtained results are shown in
Figure 3.

Table 5:  Comparison of the results achieved by our V-Net with the results of the challenge participants who were included in the official ranking. The results of the challenge participants are presented in order of ranking and come from [4]. The results achieved by the V-Net are highlighted in bold. Team |  DSCaggGTVp ($\uparrow$) |  DSCaggGTVn ($\uparrow$) |  mean DSCagg($\uparrow$)
---|---|---|---
NVAUTO |  $0.80066$ |  $0.77539$ |  $0.78802$
SJTU426 |  $0.77960$ |  $0.77604$ |  $0.77782$
NeuralRad |  $0.77485$ |  $0.76938$ |  $0.77212$
LITO |  $0.77700$ |  $0.76269$ |  $0.76984$
TheDLab |  $0.77447$ |  $0.75865$ |  $0.76656$
MAIA |  $0.75738$ |  $0.77114$ |  $0.76426$
AIRT |  $0.76689$ |  $0.73392$ |  $0.75040$
AIMers |  $0.73738$ |  $0.73431$ |  $0.73584$
SMIAL |  $0.68084$ |  $0.75098$ |  $0.71591$
Test |  $0.74499$ |  $0.68618$ |  $0.71559$
BDAV USYD |  $0.76136$ |  $0.65927$ |  $0.71031$
Our V-Net |  0.75491 |  0.65396 |  0.70443
junma |  $0.70906$ |  $0.69948$ |  $0.70427$
RokieLab |  $0.70131$ |  $0.70100$ |  $0.70115$
LMU |  $0.74460$ |  $0.65610$ |  $0.70035$
TECVICO Corp |  $0.74586$ |  $0.65069$ |  $0.69827$
RT UMCG |  $0.73741$ |  $0.65059$ |  $0.69400$
HPCAS |  $0.69786$ |  $0.66730$ |  $0.68258$
ALaGreca |  $0.72329$ |  $0.61341$ |  $0.66835$
Qurit |  $0.69553$ |  $0.57343$ |  $0.63448$
VokCow |  $0.59424$ |  $0.54988$ |  $0.57206$
MLC |  $0.46587$ |  $0.53574$ |  $0.50080$
M&H lab NU |  $0.51342$ |  $0.46557$ |  $0.48949$
Figure 3: Exemplary visualizations of the obtained results from the HECKTOR
dataset.

Figure 4: An example of a false negative prediction case for the U-Net one-
step segmentation where the model failed to detect the presence of a tumor
(marked in blue in the ground truth) near the bladder.

##  4 Discussion

Based on the results of the AutoPET dataset shown in Table 4.1, it can be
concluded that removing tumor-free images during training improves the overall
segmentation performance. Depending on the technique employed in training the
model for positively annotated images, the resulting Dice coefficients varied
between 0.39 and 0.69. Initially, it was considered that negative cases would
notably reduce the quality of segmentation, indicating that outcomes for data
containing cancer would be superior. The results obtained and statistical
analysis confirmed this statement. For a dataset consisting solely of images
with cancer lesions, the results ranged from 0.51 to 0.74. When working with a
dataset of solely cancer lesions, the model can process more focused
information that’s pertinent to tumor segmentation. This allows the model to
more readily grasp crucial features for identifying cancer. In the case of
two-stage segmentation, the initial segmentation was utilized to establish a
prediction, which was then incorporated into the model as an additional input
channel. Table 4.1 indicates a positive effect on the outcomes when compared
to the first segmentation phase. The results in the second step are
significantly better than those from the initial segmentation. The Dice
coefficient value for U-Net was 0.58 initially but increased to almost 0.60
after additional segmentation. A corresponding effect was observed on the
aggregated Dice, which yielded results of around 0.64 initially, and then
increased to roughly 0.73. This is an expected effect because adding the
initial segmentation results as an additional channel provides new contextual
information to the network, which can help the network focus on improving
relevant areas that were previously misclassified. The additional channel can
act as a kind of regularization. The model endeavors to adapt itself to new
information, thus fixating its focus on new data without any significant
deviation from it. A certain limitation is introduced in the range of possible
outputs that the next network can generate. Therefore, a fusion of the two
methods may enhance segmentation accuracy. An added benefit was the dual-
training of the network; the first segmentation phase was trained via
classical methods, while the second phase utilized a patch-based approach,
resulting in a minor increase of the required computational resources.

For the HECKTOR dataset, as shown in Table 2, the GTVn scores were notably
lower (about 7-10%) than the GTVp scores for all segmentation models. A
possible explanation for this is that GTVp was present in a significantly
larger number of images compared to GTVn, potentially making it more difficult
for the model to learn the features of GTVn and leading to an increased number
of false positives and false negatives.

The participants in the HECKTOR Grand Challenge 2022 achieved an aggregated
Dice Similarity Coefficient for segmentation within the range of
0.46587-0.80066 for GTVp and 0.46557-0.77604 for GTVn, as shown in Table 2.
Results obtained by the V-Net model were also within these ranges for the
HECKTOR challenge test set. Although the presented results of other models in
this thesis fall within these ranges, it is worth emphasizing that this
comparison is only indicative. It should be noted that the models’ performance
was evaluated only on a test set separate from the training set since no
ground truths were available for the test set. Hence, it is uncertain whether
the results of the models presented in this study would be equivalent if
tested on the HECKTOR test set. The V-Net model achieved slightly worse
results on the HECKTOR test set than on the test set that was separated from
the training set, on which all models were evaluated and compared. It can be
inferred that other models would behave similarly. Nearly half of the cases in
the HECKTOR test set originated from centers that were excluded from the
training set. This could have caused the model to have more trouble
generalizing to data it had not previously encountered. Each center’s data may
contain specific features that the model could not learn if it did not have
access to the data. Furthermore, since the HECKTOR test set was significantly
larger, there is an increased likelihood of errors, particularly false
positives and false negatives. Nonetheless, despite these challenges, the
model performed well since the variation in results was not very large.

For both datasets, the segmentation process encountered the greatest
difficulties with the metabolically active areas surrounding or located near
the analyzed structure. Such regions were inaccurately recognized by the
algorithm as parts of segmented tumors. There have also been cases where a
metabolically active structure near the tumor resulted in the tumor not being
detected. There were cases where the tumor was not detected due to the
presence of a metabolically active structure near the tumor. Such a scenario
is illustrated in Figure 4 where a cancerous growth was situated close to a
healthy, metabolically-active organ (bladder). Despite the lesion being
sizeable and noticeable, the model failed to correctly identify it as a tumor,
instead classifying it entirely as healthy tissue. In the case of the HECKTOR
dataset, metabolically active structures in the neck were also problematic for
the model and were misinterpreted by the algorithm as cancerous areas.
Detecting cancer in tumor areas with low metabolic activity posed another
challenge to the model. Metabolically active areas stand out especially on PET
images due to their intense brightness. Another factor that may have had an
impact on the results was the fact that the tumors in some of the images were
very small. In such situations, a difference of just a few voxels between the
prediction and the segmentation mask could significantly reduce the results.

The study primarily analyzed how the choice of neural network architecture
affected the obtained results. As such, the models were trained using similar
parameters, and the training images were similarly prepared. It is worth
noting that data preparation and training methods are also key factors in the
effectiveness of deep learning models, in addition to selecting the
appropriate model. The models utilized have the potential to yield better
results if the data is properly prepared and optimal training parameters and
techniques are chosen. This is confirmed by the results obtained by nnUNet
using the AutoPET dataset. The results are significantly better with the
changes introduces mostly to the initial preprocessing, augmentation, and
postprocessing. It shows the crucial role of data preparation, pre- and post-
processing, in contrast to the influence of the deep network architecture.

The size of the dataset also plays a crucial role. One of the main challenges
in the automatic analysis of medical images is posed by small datasets. Deep
learning models frequently perform better when trained on large datasets, thus
expanding the dataset has the potential to enhance results. In the context of
the size of the dataset, it is worth mentioning the UNETR model, which is
based on the transformer architecture. Transformers have recently gained
popularity and have become a new standard in many areas. In the context of
medical image analysis, there are also many studies using transformers models
and presenting promising results. A model utilizing the transformer
architecture was assessed in this study, however, it did not yield superior
outcomes when compared to convolutional neural networks. It is likely that the
limited size of the dataset hindered the transformer’s ability to perform
effectively. Transformers necessitate vast quantities of data to attain
optimal efficiency. Potentially a self-supervsied pretraining approach could
improve the UNETR results, however, this is within the scope of future work.

##  5 Conclusions

This study details several methods for tumor segmentation using PET/CT images.
Only a few of the many available segmentation approaches have been utilized in
the conducted experiments, suggesting that additional analysis may lead to
better results. This is crucial in the context of implementing deep learning
models in clinical practice, where precision and accuracy are indispensable
for correct diagnoses and optimal treatment planning. The obtained results are
comparable with the results reported in prior studies on the subject. While
there is room for improvement, the results at this stage demonstrate the
immense potential of deep learning models in the acceleration and automation
of complex medical analysis. Continued research in this area is necessary to
further improve the effectiveness of diagnosis. Further work should focus on
developing more advanced segmentation techniques and research on their
integration with existing diagnostic systems. Nevertheless, the key takeaways
from the results are connected with the importance of initial data
preparation, augmentation and pre-processing that are more important that the
deep architecture itself.

## Acknowledgements

We gratefully acknowledge Poland’s high-performance computing infrastructure
PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and
support within computational grant no. PLG/2023/016239.

## References

  * [1] https://github.com/MIC-DKFZ/nnUNetUNETR.  Accessed: 21.08.2023.
  * [2] https://hecktor.grand-challenge.org/Data/Accessed: 2023-12-03.
  * [3] https://github.com/voreille/hecktor/blob/master/notebooks/evaluate_segmentation2022.ipynb.  Accessed: 16.08.2023.
  * [4] V. Andrearczyk, V. Oreiller, M. Abobakr, A. Akhavanallaf, P. Balermpas, S. Boughdad, L. Capriotti, J. Castelli, C. Cheze Le Rest, P. Decazes, R. Correia, D. El-Habashy, H. Elhalawani, C. Fuller, M. Jreige, Y. Khamis, A. Greca, A. Mohamed, M. Naser, and A. Depeursinge.  Overview of the hecktor challenge at miccai 2022: Automatic head and neck tumor segmentation and outcome prediction in pet/ct.  volume 13626, pages 1–30. 03 2023.
  * [5] V. Andrearczyk, V. Oreiller, S. Boughdad, C. C. L. Rest, H. Elhalawani, M. Jreige, J. O. Prior, M. Vallières, D. Visvikis, M. Hatt, et al.  Overview of the hecktor challenge at miccai 2021: automatic head and neck tumor segmentation and outcome prediction in pet/ct images.  In 3D head and neck tumor segmentation in PET/CT challenge, pages 1–37. Springer, 2021.
  * [6] V. Andrearczyk, V. Oreiller, M. Jreige, M. Vallieres, J. Castelli, H. Elhalawani, S. Boughdad, J. O. Prior, and A. Depeursinge.  Overview of the hecktor challenge at miccai 2020: automatic head and neck tumor segmentation in pet/ct.  In Head and Neck Tumor Segmentation: First Challenge, HECKTOR 2020, Held in Conjunction with MICCAI 2020, Lima, Peru, October 4, 2020, Proceedings 1, pages 1–21. Springer, 2021.
  * [7] R. Azad, E. K. Aghdam, A. Rauland, Y. Jia, A. H. Avval, A. Bozorgpour, S. Karimijafarbigloo, J. P. Cohen, E. Adeli, and D. Merhof.  Medical image segmentation review: The success of u-net.  arXiv preprint arXiv:2211.14830, 2022.
  * [8] C. Bauer, S. Sun, W. Sun, J. Otis, A. Wallace, B. J. Smith, J. J. Sunderland, M. M. Graham, M. Sonka, J. M. Buatti, et al.  Automated measurement of uptake in cerebellum, liver, and aortic arch in full-body fdg pet/ct scans.  Medical physics, 39(6Part1):3112–3123, 2012.
  * [9] F. Bray, M. Laversanne, E. Weiderpass, and I. Soerjomataram.  The ever-increasing importance of cancer as a leading cause of premature death worldwide.  Cancer, 127(16):3029–3030, 2021.
  * [10] L. Q. Chow.  Head and neck cancer.  New England Journal of Medicine, 382(1):60–72, 2020.
  * [11] M. Contributors.  Unetr implementation in monai, 2023.  Accessed: 21.08.2023.
  * [12] P.-H. Conze, G. Andrade-Miranda, V. K. Singh, V. Jaouen, and D. Visvikis.  Current and emerging trends in medical image segmentation with deep learning.  IEEE Transactions on Radiation and Plasma Medical Sciences, 2023.
  * [13] B. Efron and R. J. Tibshirani.  An introduction to the bootstrap.  CRC press, 1994.
  * [14] C. Franzese, D. Dei, N. Lambri, M. A. Teriaca, M. Badalamenti, L. Crespi, S. Tomatis, D. Loiacono, P. Mancosu, and M. Scorsetti.  Enhancing radiotherapy workflow for head and neck cancer with artificial intelligence: A systematic review.  Journal of Personalized Medicine, 13(6):946, 2023.
  * [15] S. Gatidis, M. Früh, M. Fabritius, S. Gu, K. Nikolaou, C. La Fougère, J. Ye, J. He, Y. Peng, L. Bi, et al.  The autopet challenge: Towards fully automated lesion segmentation in oncologic pet/ct imaging.  2023\.
  * [16] S. Gatidis, T. Hepp, M. Früh, C. La Fougère, K. Nikolaou, C. Pfannenberg, B. Schölkopf, T. Küstner, C. Cyran, and D. Rubin.  A whole-body fdg-pet/ct dataset with manually annotated tumor lesions.  Scientific Data, 9(1):601, 2022.
  * [17] M. Illimoottil and D. Ginat.  Recent advances in deep learning and medical imaging for head and neck cancer treatment: Mri, ct, and pet scans.  Cancers, 15(13):3267, 2023.
  * [18] F. Isensee, P. F. Jaeger, S. A. Kohl, J. Petersen, and K. H. Maier-Hein.  nnu-net: a self-configuring method for deep learning-based biomedical image segmentation.  Nature methods, 18(2):203–211, 2021.
  * [19] H. Jiang, Z. Diao, T. Shi, Y. Zhou, F. Wang, W. Hu, X. Zhu, S. Luo, G. Tong, and Y.-D. Yao.  A review of deep learning-based multiple-lesion recognition from medical images: classification, detection and segmentation.  Computers in Biology and Medicine, page 106726, 2023.
  * [20] S. Li, J. Liu, Z. Wang, Z. Cao, Y. Yang, B. Wang, S. Xu, L. Lu, M. Iqbal Saripan, X. Zhang, et al.  Application of pet/ct-based deep learning radiomics in head and neck cancer prognosis: a systematic review.  Radiology Science, 2022.
  * [21] P. Malhotra, S. Gupta, D. Koundal, A. Zaguia, W. Enbeyle, et al.  Deep neural networks for medical image segmentation.  Journal of Healthcare Engineering, 2022, 2022.
  * [22] E. McCrum-Gardner.  Which is the correct statistical test to use?  British Journal of Oral and Maxillofacial Surgery, 46(1):38–41, 2008.
  * [23] F. Milletari, N. Navab, and S.-A. Ahmadi.  V-net: Fully convolutional neural networks for volumetric medical image segmentation.  In 2016 fourth international conference on 3D vision (3DV), pages 565–571. Ieee, 2016.
  * [24] Y. Mo, Y. Wu, X. Yang, F. Liu, and Y. Liao.  Review the state-of-the-art technologies of semantic segmentation based on deep learning.  Neurocomputing, 493:626–646, 2022.
  * [25] K. Munir, H. Elahi, A. Ayub, F. Frezza, and A. Rizzi.  Cancer diagnosis using deep learning: a bibliographic review.  Cancers, 11(9):1235, 2019.
  * [26] K. Noguchi, Y. R. Gel, and C. R. Duguay.  Bootstrap-based tests for trends in hydrological time series, with application to ice phenology data.  Journal of Hydrology, 410(3-4):150–161, 2011.
  * [27] V. Oreiller, V. Andrearczyk, M. Jreige, S. Boughdad, H. Elhalawani, J. Castelli, M. Vallieres, S. Zhu, J. Xie, Y. Peng, et al.  Head and neck tumor segmentation in pet/ct: the hecktor challenge.  Medical image analysis, 77:102336, 2022.
  * [28] F. Pérez-García, R. Sparks, and S. Ourselin.  Torchio: a python library for efficient loading, preprocessing, augmentation and patch-based sampling of medical images in deep learning.  Computer Methods and Programs in Biomedicine, 208:106236, 2021.
  * [29] L. Rebaud, T. Escobar, F. Khalid, K. Girum, and I. Buvat.  Simplicity is all you need: out-of-the-box nnunet followed by binary-weighted radiomic model for segmentation and outcome prediction in head and neck pet/ct.  In 3D Head and Neck Tumor Segmentation in PET/CT Challenge, pages 121–134. Springer, 2022.
  * [30] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * [31] F. Shi, W. Hu, J. Wu, M. Han, J. Wang, W. Zhang, Q. Zhou, J. Zhou, Y. Wei, Y. Shao, et al.  Deep learning empowered volume delineation of whole-body organs-at-risk for accelerated radiotherapy.  Nature Communications, 13(1):6566, 2022.
  * [32] I. Shiri, M. Amini, F. Yousefirizi, A. Vafaei Sadr, G. Hajianfar, Y. Salimi, Z. Mansouri, E. Jenabi, M. Maghsudi, I. Mainta, et al.  Information fusion for fully automated segmentation of head and neck tumors from pet and ct images.  Medical Physics, 51(1):319–333, 2024.
  * [33] H. Thisanke, C. Deshan, K. Chamith, S. Seneviratne, R. Vidanaarachchi, and D. Herath.  Semantic segmentation using vision transformers: A survey.  Engineering Applications of Artificial Intelligence, 126:106669, 2023.
  * [34] K. A. Tran, O. Kondrashova, A. Bradley, E. D. Williams, J. V. Pearson, and N. Waddell.  Deep learning in cancer diagnosis, prognosis and treatment selection.  Genome Medicine, 13(1):1–17, 2021.
  * [35] X. Wang and B.-b. Li.  Deep learning in head and neck tumor multiomics diagnosis and analysis: review of the literature.  Frontiers in Genetics, 12:624820, 2021.
  * [36] Y. Wang, E. Lombardo, L. Huang, M. Avanzo, G. Fanetti, G. Franchin, S. Zschaeck, J. Weingärtner, C. Belka, M. Riboldi, et al.  Comparison of deep learning networks for fully automated head and neck tumor delineation on multi-centric pet/ct images.  Radiation Oncology, 19(1):3, 2024.
  * [37] H. Xiao, L. Li, Q. Liu, X. Zhu, and Q. Zhang.  Transformers in medical image segmentation: A review.  Biomedical Signal Processing and Control, 84:104791, 2023.
  * [38] H. Xue, Y. Yao, and Y. Teng.  Multi-modal tumor segmentation methods based on deep learning: a narrative review.  Quantitative Imaging in Medicine and Surgery, 14(1):1122, 2024.
  * [39] X. Zhou, R. Takayama, S. Wang, T. Hara, and H. Fujita.  Deep learning of the sectional appearances of 3d ct images for anatomical structure segmentation based on an fcn voting method.  Medical physics, 44(10):5221–5233, 2017.


# Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”

Jae-Mo Lihm<EMAIL_ADDRESS>Department of Physics and Astronomy, Seoul
National University, Seoul 08826, Korea

###### Abstract

In a recent paper, Ibañez-Azpiroz et al. [Phys. Rev. B 97, 245143 (2018)]
derive a band-truncation-error-free formula for calculating the generalized
derivative of the interband dipole matrix using Wannier interpolation. In
practice, the denominators involving intermediate states are regularized by
introducing a finite broadening parameter. In this Comment, I show that when a
finite broadening parameter is used, a correction term must be added to the
generalized derivative to obtain results that are independent of the phase
convention for the Bloch sums.

Ibañez-Azpiroz  _et al._ (2018) derive a truncation-error-free formula for the
Wannier interpolation of the generalized derivative of the interband dipole
matrix. This result enables an efficient and accurate calculation of the
shift-current response.

In practice, to avoid numerical problems due to near-degenerate states, the
energy denominator involving intermediate states is regularized by introducing
a broadening parameter $\eta$ (see Eq. (38) of Ibañez-Azpiroz  _et al._
(2018)). Within the Wannier-interpolation scheme, the regularized version of
the generalized derivative of the interband dipole matrix element is given
from Eqs. (36) and (38) of Ibañez-Azpiroz  _et al._ (2018) by

$\displaystyle\widetilde{r}^{a;b}_{nm}(\eta)=\widetilde{r}^{a;b}_{nm}(0)+\frac{i}{\omega_{nm}}$
$\displaystyle\times\sum_{p\neq
n,m}^{M}\Bigg{\\{}\real\left[\frac{-i\eta}{\omega_{pm}(\omega_{pm}+i\eta)}\right](\mathbbm{v}_{np}^{a}+i\omega_{nm}\mathbbm{a}_{np}^{a})\mathbbm{v}_{pm}^{b}$
$\displaystyle\quad-\real\left[\frac{-i\eta}{\omega_{np}(\omega_{np}+i\eta)}\right]\mathbbm{v}_{np}^{b}(\mathbbm{v}_{pm}^{a}+i\omega_{nm}\mathbbm{a}_{pm}^{a})\Bigg{\\}},$
(1)

where $M$ is the number of Wannier functions. Here,
$\widetilde{r}^{a;b}_{nm}(0)$ refers to the case without regularization, which
equals to $r^{a;b}_{nm}$ in Eq. (36) of Ibañez-Azpiroz  _et al._ (2018). I use
tilde to distinguish the quantity obtained from the Wannier interpolation
scheme of Ibañez-Azpiroz  _et al._ (2018) from a direct calculation using
Kohn-Sham wavefunctions.

In this Comment, I show that when $\eta$ is nonzero,
$\widetilde{r}^{a;b}_{nm}(\eta)$ defined in Eq. (Comment on “Ab initio
calculation of the shift photocurrent by Wannier interpolation”), as well as
the corresponding shift-current spectrum, depends on the phase convention for
the Bloch sums. I derive a correction term that should be added to Eq.
(Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”) to obtain a phase-convention-independent result.

I consider two phase conventions, the tight-binding (TB) convention [Eq. (14)
of Ibañez-Azpiroz  _et al._ (2018)]

$\ket{u_{{\bm{\mathbf{k}}}j}^{\mathrm{(W,\,TB)}}}=\sum_{{\bm{\mathbf{R}}}}e^{-i{\bm{\mathbf{k}}}\cdot(\hat{{\bm{\mathbf{r}}}}-{\bm{\mathbf{R}}}-{\bm{\mathbf{\tau}}}_{j})}\ket{{\bm{\mathbf{R}}}j}$
(2)

and the Wannier90 (W90) convention

$\ket{u_{{\bm{\mathbf{k}}}j}^{\mathrm{(W,\,W90)}}}=\sum_{{\bm{\mathbf{R}}}}e^{-i{\bm{\mathbf{k}}}\cdot(\hat{{\bm{\mathbf{r}}}}-{\bm{\mathbf{R}}})}\ket{{\bm{\mathbf{R}}}j}.$
(3)

Here, ${\bm{\mathbf{\tau}}}_{j}$ is the position of the $j$-th atom of the
unit cell. The “internal” matrix elements that appear in Eq. (Comment on “Ab
initio calculation of the shift photocurrent by Wannier interpolation”) such
as $\mathbbm{v}$ and $\mathbbm{a}$ depend on the phase convention. In the
tight-binding (TB) convention, one finds

$\partial^{a}H_{ij}^{\rm(W,TB)}=\sum_{\bm{\mathbf{R}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R}}}+{\bm{\mathbf{\tau}}}_{j}-{\bm{\mathbf{\tau}}}_{i})}i(R^{a}+\tau^{a}_{j}-\tau^{a}_{i})H_{ij;{\bm{\mathbf{R}}}},$
(4)
$A_{ij}^{a\mathrm{(W,TB)}}=\sum_{\bm{\mathbf{R}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R}}}+{\bm{\mathbf{\tau}}}_{j}-{\bm{\mathbf{\tau}}}_{i})}r^{a}_{ij;{\bm{\mathbf{R}}}}-\tau^{a}_{i}\delta_{ij}.$
(5)

In the W90 convention, one finds

$\partial^{a}H_{ij}^{\rm(W,W90)}=\sum_{{\bm{\mathbf{R}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}iR^{a}H_{ij;{\bm{\mathbf{R}}}},$
(6)
$A_{ij}^{a\mathrm{(W,W90)}}=\sum_{\bm{\mathbf{R}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}r^{a}_{ij;{\bm{\mathbf{R}}}}.$
(7)

The unitary matrices that diagonalize the Hamiltonians in the two conventions
can be related as

$U^{(\mathrm{TB})}_{in}=e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{\tau}}}_{i}}U^{(\mathrm{W90})}_{in}.$
(8)

Comparing Eqs. (6, 7) with Eqs. (4, 5), one can show that the internal matrix
elements in the two conventions are related as

$\mathbbm{v}^{a,(\mathrm{TB})}_{nm}=\mathbbm{v}^{a,(\mathrm{W90})}_{nm}+i\omega_{nm}\bbtau^{a,(\mathrm{TB})}_{nm}$
(9)

and

$\mathbbm{a}^{a,(\mathrm{TB})}_{nm}=\mathbbm{a}^{a,(\mathrm{W90})}_{nm}-\bbtau^{a,(\mathrm{TB})}_{nm}$
(10)

where

$\bbtau^{a,(\mathrm{TB})}_{nm}\equiv\sum_{i}U^{(\mathrm{TB})\dagger}_{ni}\tau^{a}_{i}U^{(\mathrm{TB})}_{im}.$
(11)

Thus, if there is an atom that is not at the origin of the unit cell, the two
phase conventions give different internal matrix elements. Ibañez-Azpiroz  _et
al._ (2018) showed that when $\eta=0$, the additional terms involving $\bbtau$
cancel out so that the two conventions give identical shift-current spectrum.
Here, I show that such cancellation is incomplete when $\eta\neq 0$.

First, let us define the generalized derivative with finite broadening $\eta$
based on Bloch wavefunctions as follows:

$\displaystyle
r^{a;b}_{nm}(\eta)=\frac{i}{\omega_{nm}}\Bigg{\\{}\frac{v_{nm}^{a}\Delta_{nm}^{b}+v_{nm}^{b}\Delta_{nm}^{a}}{\omega_{nm}}-w^{ab}_{nm}$
$\displaystyle+\sum_{p\neq
n,m}^{\infty}\Big{[}\real\left(\frac{1}{\omega_{pm}+i\eta}\right)v_{np}^{a}v_{pm}^{b}$
$\displaystyle\quad-\real\left(\frac{1}{\omega_{np}+i\eta}\right)v_{np}^{b}v_{pm}^{a}\Big{]}\Bigg{\\}}.$
(12)

When $\eta=0$, $r^{a;b}_{nm}(\eta)$ reduces to Eq. (12) of Ibañez-Azpiroz  _et
al._ (2018). Our goal is to obtain a Wannier interpolation formula for
$r^{a;b}_{nm}(\eta)$ which holds true at finite $\eta$.

We find

$\displaystyle
r^{a;b}_{nm}(\eta)-r^{a;b}_{nm}(0)=\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{\infty}A_{nm,p}^{a,b}(\eta)$ (13)

where we defined

$\displaystyle A_{nm,p}^{a,b}(\eta)$
$\displaystyle=\real\left(\frac{1}{\omega_{pm}+i\eta}-\frac{1}{\omega_{pm}}\right)v_{np}^{a}v_{pm}^{b}$
$\displaystyle-\real\left(\frac{1}{\omega_{np}+i\eta}-\frac{1}{\omega_{np}}\right)v_{np}^{b}v_{pm}^{a}.$
(14)

In principle, Eq. (13) requires a sum over an infinite number of unoccupied
bands and thus cannot be computed using Wannier interpolation. However, note
that

$\displaystyle\real\left(\frac{1}{\omega_{pm}+i\eta}-\frac{1}{\omega_{pm}}\right)$
$\displaystyle=\real\left[\frac{-i\eta}{\omega_{pm}(\omega_{pm}+i\eta)}\right]$
$\displaystyle=-\frac{\eta^{2}}{\omega_{pm}(\omega_{pm}^{2}+\eta^{2})}$ (15)

is very small as long as $\absolutevalue{\omega_{pm}}$ is much larger than
$\eta$. Also, when computing the shift-current response, $r^{a;b}_{nm}$ is
computed only for the low-energy states $n$ and $m$ which are a pair of
occupied and unoccupied states that can be resonantly excited by an external
electric field. Thus, Eq. (13) converges rapidly with the number of states
included in the sum. We exploit this rapid convergence to approximate the
infinite sum in Eq. (13) by a finite sum over Wannier-interpolated bands

$r^{a;b}_{nm}(\eta)-r^{a;b}_{nm}(0)\approx\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{M}A_{nm,p}^{a,b}(\eta),$ (16)

where $A_{nm,p}^{a,b}(\eta)$ [Eq. (Comment on “Ab initio calculation of the
shift photocurrent by Wannier interpolation”)] is calculated using the Wannier
interpolation of energy and velocity matrix elements. Later, I demonstrate the
validity of this approximation using numerical calculation (see Fig. 3).

Now, let us study the relation between Eq. (Comment on “Ab initio calculation
of the shift photocurrent by Wannier interpolation”) and Eq. (Comment on “Ab
initio calculation of the shift photocurrent by Wannier interpolation”). The
velocity matrix element can be written as

$v_{nm}^{a}=\mathbbm{v}_{nm}^{a}+i\omega_{nm}\mathbbm{a}_{nm}^{a}.$ (17)

This expression can be obtained from Eqs. (11), (22), and (23) of Ibañez-
Azpiroz  _et al._ (2018). From Eqs. (9, 10), one can easily show that
$v_{nm}^{a}$ does not depend on the phase convention. By using Eqs. (Comment
on “Ab initio calculation of the shift photocurrent by Wannier interpolation”,
Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”), we find

$\displaystyle\frac{i}{\omega_{nm}}\sum_{p\neq n,m}^{M}A_{nm,p}^{a,b}(\eta)=$
$\displaystyle\widetilde{r}^{a;b}_{nm}(\eta)-\widetilde{r}^{a;b}_{nm}(0)+\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left[\frac{-i\eta}{\omega_{pm}(\omega_{pm}+i\eta)}\right]\left[v^{a}_{np}v^{b}_{pm}-(\mathbbm{v}_{np}^{a}+i\omega_{nm}\mathbbm{a}_{np}^{a})\mathbbm{v}_{pm}^{b}\right]$
$\displaystyle\hskip 72.00008pt-\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left[\frac{-i\eta}{\omega_{np}(\omega_{np}+i\eta)}\right]\left[v^{b}_{np}v^{a}_{pm}-\mathbbm{v}_{np}^{b}(\mathbbm{v}_{pm}^{a}+i\omega_{nm}\mathbbm{a}_{pm}^{a})\right]$
$\displaystyle=$
$\displaystyle\widetilde{r}^{a;b}_{nm}(\eta)-\widetilde{r}^{a;b}_{nm}(0)+\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left[\frac{-i\eta}{\omega_{pm}(\omega_{pm}+i\eta)}\right]\left[i\omega_{mp}\mathbbm{a}_{np}^{a}\mathbbm{v}_{pm}^{b}+v^{a}_{np}i\omega_{pm}\mathbbm{a}_{pm}^{b}\right]$
$\displaystyle\hskip 72.00008pt-\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left[\frac{-i\eta}{\omega_{np}(\omega_{np}+i\eta)}\right]\left[i\omega_{pn}\mathbbm{v}_{np}^{b}\mathbbm{a}_{pm}^{a}+i\omega_{np}\mathbbm{a}_{np}^{b}v^{a}_{pm}\right]$
$\displaystyle=$
$\displaystyle\widetilde{r}^{a;b}_{nm}(\eta)-\widetilde{r}^{a;b}_{nm}(0)+\frac{1}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left(\frac{-i\eta}{\omega_{pm}+i\eta}\right)\left[\mathbbm{a}_{np}^{a}\mathbbm{v}_{pm}^{b}-(\mathbbm{v}_{np}^{a}+i\omega_{np}\mathbbm{a}_{np}^{a})\mathbbm{a}_{pm}^{b}\right]$
$\displaystyle\hskip 72.00008pt-\frac{1}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\real\left(\frac{-i\eta}{\omega_{np}+i\eta}\right)\left[\mathbbm{v}_{np}^{b}\mathbbm{a}_{pm}^{a}-\mathbbm{a}_{np}^{b}(\mathbbm{v}_{pm}^{a}+i\omega_{pm}\mathbbm{a}_{pm}^{a})\right].$
(18)

Using Eqs. (16, Comment on “Ab initio calculation of the shift photocurrent by
Wannier interpolation”) and $r^{a;b}_{nm}(0)=\widetilde{r}^{a;b}_{nm}(0)$
which is proven in Ibañez-Azpiroz  _et al._ (2018), we find

$\displaystyle r^{a;b}_{nm}(\eta)\approx\widetilde{r}^{a;b}_{nm}(\eta)-$
$\displaystyle\frac{1}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\frac{\eta^{2}}{\omega_{pm}^{2}+\eta^{2}}\left[\mathbbm{a}_{np}^{a}\mathbbm{v}_{pm}^{b}-(\mathbbm{v}_{np}^{a}+i\omega_{np}\mathbbm{a}_{np}^{a})\mathbbm{a}_{pm}^{b}\right]$
$\displaystyle+$ $\displaystyle\frac{1}{\omega_{nm}}\sum_{p\neq
n,m}^{M}\frac{\eta^{2}}{\omega_{np}^{2}+\eta^{2}}\left[\mathbbm{v}_{np}^{b}\mathbbm{a}_{pm}^{a}-\mathbbm{a}_{np}^{b}(\mathbbm{v}_{pm}^{a}+i\omega_{pm}\mathbbm{a}_{pm}^{a})\right].$
(19)

The right-hand side of Eq. (Comment on “Ab initio calculation of the shift
photocurrent by Wannier interpolation”) is the desired phase-convention-
independent Wannier interpolation formula for the generalized derivative of
the interband dipole matrix element at finite $\eta$. This expression is
independent of the phase convention. The reason is as follows. Equation
(Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”) is phase-convention independent because $A^{a,b}_{nm,p}(\eta)$
[Eq. (Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”)] is defined only in terms of phase-convention independent
matrix elements. Also, $\widetilde{r}^{a;b}_{nm}(0)$ is phase-convention
independent because it equals $r^{a;b}_{nm}(0)$ which is defined in terms of
Kohn-Sham wavefunctions. Since the right-hand side of Eq. (Comment on “Ab
initio calculation of the shift photocurrent by Wannier interpolation”) is the
sum of Eq. (Comment on “Ab initio calculation of the shift photocurrent by
Wannier interpolation”) and $\widetilde{r}^{a;b}_{nm}(0)$, it is also
independent of the phase convention.

When $\eta=0$, the correction terms in Eq. (Comment on “Ab initio calculation
of the shift photocurrent by Wannier interpolation”) vanish as expected. When
$\eta$ is nonzero, the correction terms are finite. Also, we find that these
correction terms depend on the phase convention used for the Bloch sums. Thus,
if this correction is not applied, the calculated $r^{a;b}_{nm}(\eta)$ matrix
elements, as well as the shift-current spectrum, becomes dependent on the
phase convention.

Figure 1:  Difference of the $zzz$ component of the shift-current spectrum of
monolayer GeS computed using the W90 and TB conventions (a) without and (b)
with the finite-$\eta$ correction [Eq. (Comment on “Ab initio calculation of
the shift photocurrent by Wannier interpolation”)].  Figure 2:  Change in the
shift-current spectrum of monolayer GeS computed with the (a) W90 convention
and the (b) TB convention due to the finite-$\eta$ correction:
$\Delta\sigma=\sigma^{\rm w/corr.}-\sigma^{\rm w/o\ corr.}$.

Now, I demonstrate these findings with numerical calculations. I implemented
the correction formula Eq. (Comment on “Ab initio calculation of the shift
photocurrent by Wannier interpolation”) in Wannier90 and calculated the shift-
current spectrum of monolayer GeS. For the density functional theory
calculation, I used fully relativistic ONCV pseudopotentials Hamann (2013)
taken from the PseudoDojo library (v0.4) van Setten  _et al._ (2018) and a
wavefunction kinetic energy cutoff of 80 Ry. Other computational parameters
are the same as in Ibañez-Azpiroz  _et al._ (2018) unless otherwise noted.

Figure 1 shows the difference between the shift-current spectra calculated
using the two conventions. Without the correction term, there is a small but
nonzero discrepancy between the two spectra. The discrepancy increases with
$\eta$. When the correction term is added, the shift-current spectra
calculated using the TB and the W90 conventions perfectly agree for all values
of $\eta$. This result demonstrates the phase-convention independence of Eq.
(Comment on “Ab initio calculation of the shift photocurrent by Wannier
interpolation”).

Figure 2 shows the change in the shift-current spectrum due to the correction
term. For $\eta$=40 meV, which is used in Ibañez-Azpiroz  _et al._ (2018), the
change in the shift-current spectrum due to the correction is below 0.3
$\mathrm{\mu A/V^{2}}$. This correction is nonzero but will have little effect
on the calculated shift-current spectrum [Fig. 5(a) of Ibañez-Azpiroz  _et
al._ (2018)] whose magnitude is around 50 $\mathrm{\mu A/V^{2}}$.

Finally, I test the validity of approximating the infinite sum over bands in
Eq. (16) by a finite sum over Wannier-interpolated bands as Eq. (16). Let us
define the following approximant of Eq. (13):

$f_{n,m}^{a;b}(E_{\rm cutoff};\eta)=\frac{i}{\omega_{nm}}\sum_{p\neq
n,m}^{\varepsilon_{p}\leq E_{\rm cutoff}}A_{nm,p}^{a,b}(\eta).$ (20)

Here, I replaced the infinite sum over $p$ in Eq. (13) with a finite sum for
states with $\varepsilon_{p}\leq E_{\rm cutoff}$. Equation (13) is recovered
in the limit $E_{\rm cutoff}\rightarrow\infty$. In a Wannier interpolation,
only the bands inside the frozen window for Wannierization are accurately
interpolated. So, to use Wannier interpolation, we choose $E_{\rm cutoff}$ to
be $E_{\rm froz.}$, the upper bound of the frozen window for Wannierization.
Then, the truncation error in evaluating Eq. (13) is

$\delta r^{a;b}_{nm}(\eta)=f_{n,m}^{a;b}(\infty;\eta)-f_{n,m}^{a;b}(E_{\rm
froz.};\eta).$ (21)

Figure 3:  Truncation error in the finite-$\eta$ correction of the shift
current spectrum of monolayer GeS calculated by replacing $r^{a;b}_{nm}$ in
the shift current formula [Eqs. (5, 8) of Ibañez-Azpiroz  _et al._ (2018)]
with $\delta r^{a;b}_{nm}(\eta)$ [Eq. (21)]. We used a (a)
100$\times$100$\times$1 and a (b) 200$\times$200$\times$1 $k$-point grid.

To study the magnitude of the error in the shift-current spectra due to this
truncation, I calculated the shift-current spectrum with $r^{a;b}_{nm}$ in the
shift current formula [Eqs. (5, 8) of Ibañez-Azpiroz  _et al._ (2018)]
replaced with $\delta r^{a;b}_{nm}(\eta)$ [Eq. (21)]. Since the high-energy
bands need to be explicitly included, I performed the calculation at the DFT
level. The infinity in Eq. (21) was replaced with the conduction band minimum
energy plus 20 eV, as the contribution from bands with even higher energies
was found to be negligible. Due to the high computational cost, I used
$k$-point grids coarser than the 1000$\times$1000$\times$1 grid used in the
calculation of the shift-current spectra. Figure 3 shows the truncation error
of the shift-current spectra. The results indicate that the magnitude of the
truncation error is almost converged at a 200$\times$200$\times$1 grid. The
truncation error of the shift-current spectra is at most 0.012 $\mathrm{\mu
A/V^{2}}$, even at a large $\eta$ of 100 meV up to $\omega$=4 eV when the
upper limit of the frozen window was set to be 5 eV above the top of the
valence band Ibañez-Azpiroz  _et al._ (2018). For each $\eta$, the truncation
error is much smaller than the finite-$\eta$ correction of the shift-current
spectrum shown in Fig. 2. This analysis shows that the states with energy
above the frozen window give negligible contributions to the finite-$\eta$
correction of the shift-current spectra. Therefore, we can safely approximate
the infinite sum in Eq. (13) by a sum over Wannier-interpolated bands as in
Eq. (16).

Note that the use of $\eta$ is necessary to avoid numerical problems due to
near-degenerate states. As one uses smaller $\eta$, a larger number of $k$
points need to be used to converge the calculated shift-current spectrum. In
other words, the number of $k$ points used gives the lower bound of $\eta$ for
the calculated spectrum to remain stable.

I also note that I fixed a small inconsistency between the shift current
implementation of Wannier90 and the formalism of Ibañez-Azpiroz  _et al._
(2018). In the Wannier90 implementation, the broadening $\eta$ was included in
denominators that do not contain any intermediate states, such as
$1/\omega_{nm}$ in Eq. (Comment on “Ab initio calculation of the shift
photocurrent by Wannier interpolation”), while the formalism Ibañez-Azpiroz
_et al._ (2018) does not include $\eta$ in such cases. I changed the Wannier90
code so that it is consistent with Ibañez-Azpiroz  _et al._ (2018). Without
this change, the exact agreement between the two phase conventions was not
reached. All results in this Comment were obtained after this change.

To summarize, I showed that when using a finite broadening parameter $\eta$ to
calculate the generalized derivative of the interband dipole matrix, one needs
to add a correction term [Eq. (Comment on “Ab initio calculation of the shift
photocurrent by Wannier interpolation”)] to Eq. (36) of Ibañez-Azpiroz  _et
al._ (2018). Otherwise, the calculated shift-current spectrum becomes
dependent on the phase convention for the Bloch sums.

###### Acknowledgements.

Computational resources were provided by KISTI Supercomputing Center
(KSC-2020-INO-0078).

## References

  * Ibañez-Azpiroz  _et al._ (2018) Julen Ibañez-Azpiroz, Stepan S. Tsirkin, and Ivo Souza, “ _Ab initio_ calculation of the shift photocurrent by Wannier interpolation,” Physical Review B 97, 245143 (2018).
  * Hamann (2013) D. R. Hamann, “Optimized norm-conserving Vanderbilt pseudopotentials,” Physical Review B 88, 085117 (2013).
  * van Setten  _et al._ (2018) M.J. van Setten, M. Giantomassi, E. Bousquet, M.J. Verstraete, D.R. Hamann, X. Gonze, and G.-M. Rignanese, “The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table,” Computer Physics Communications 226, 39–54 (2018).