# Anomalies and parities for quintessential and ultra-light axions

Jihn E. Kim  Department of Physics, Kyung Hee University, 26 Gyungheedaero,
Dongdaemun-Gu, Seoul 02447, Republic of Korea

###### Abstract

We discuss the energy scales of the explicit breaking terms of the global
symmetries U(1)${}_{\rm global}\times\rm SU(2)_{\rm gauge}^{2}$ needed for the
quinessential axion (QA) and the ultra-light axion (ULA). The appropriate
scale of QA is about $1.48\times 10^{8\,}\,\textrm{GeV}$.

Quintessential axion, Dark energy, QCD axion, SU(2) anomaly.

###### pacs:

11.25.Mj, 11.30.Er, 11.25.Wx, 12.60.Jv

##  I Introduction

Light pseudoscalar particles are the lampposts to the road leading to physics
scales much above their masses. The well-known example is the pion triplet
which guided toward the SU(2)${}_{L}\times$SU(2)R chiral symmetry at the
strong interaction scale Nambu61 . Another is the very light axion which hints
the intermediate scale KSVZ1 ; KNS18 . The scale where the restored symmetry
recovered is, $E\gtrsim f\sim\Lambda^{2}/m$, for the pseudoscalar mass $m$,
where $\Lambda^{4}$ is a typical energy density contributed by such a
pseudoscalar. Before obtaining mass by the energy density perturbation
$\Lambda^{4}$, these pseudoscalar degrees are phase degrees $\theta(x)$ of
some unitary operators such that they have kinetic energy terms in quantum
field theory below the pseudoscalar defining scale $f$. Above the scale $f$,
the phase is not a dynamical field, i.e. not depending on $x$, but it still
represents a phase direction of the global symmetry. Gravity does not respect
global symmetries. Global symmetries are destined to be broken by gravity.

If the pseudoscalar mass $m$ is less than 1 MeV, it can decay to two photons
and two neutrinos. For these two rare processes, the $\pi^{0}$ decays give
$\Gamma(\pi^{0}\to 2\gamma)\simeq(\alpha_{\rm
em}^{2}/64\pi^{3})(M_{\pi^{0}}/f_{\pi})^{2}M_{\pi^{0}}$ (at leading order
Holstein02 ) and $1.2\times 10^{-5}\Gamma(\pi^{0}\to 2\gamma)$ Shrock79 ,
respectively. For a pseudoscalar to be still present in our Universe, we
require $\Gamma^{-1}>4.3\times 10^{17}{\rm s}=1/(1.53\times
10^{-33}\,\mathrm{eV})$. Using the $f-m$ relation, we get an idea on $m$ from
the condition that it survives until now. If $\Lambda$ is the QCD scale

$\begin{split}\Gamma^{-1}\simeq\frac{64\pi^{3}\,\Lambda_{\rm
QCD}^{4}}{\alpha_{\rm em}^{2}m^{4}}\frac{1}{m}>\frac{1}{1.53\times
10^{-33}\,\mathrm{eV}}\end{split}$ (1)

or $m<65\,\,\mathrm{eV}$ for $\Lambda_{\rm QCD}=380\,\,\mathrm{MeV}$.111For
the very light axion, the specific coupling reduces it further to 24 eV
KimRMP10 . If we take $\Lambda^{4}$ as the current energy density of the
Universe, $(0.003\,\,\mathrm{eV})^{4}$, then pseudoscalars with mass less than
$3.3\times 10^{-4}\,\,\mathrm{eV}$ survives until now. Therefore, very light
pseudoscalars attracted a great deal of attention in cosmology. Pseudoscalars
for $m<3.3\times 10^{-4}\,\,\mathrm{eV}$ can contribute to the current energy
density and structure formation. In this context, quinessential axion (QA)
Carroll98 ; Hill02 ; KN03  was suggested for dark energy (DE), and ultra-light
axion (ULA) KimPRD16  was suggested for the galactic scale structures Witten17
. The precursers of these pseudoscalars are ‘general quintessence’, not
pinpointing to the pseudoscalars quinessence .

For a pseudoscalar created by spontaneous breaking of a global symmetry, there
is no pseudoscalar above the defining scale $f$. Below $f$, the fundamental
scale for the pseudoscalar decay-interactions is $f$ and hence $f$ is called
the decay constant. On the other hand, if a pseudoscalar is present in the
theory, then $f$ is the defining scale of that theory. In string theory, there
are antisymmetric tensor fields which are pseudoscalars, and hence $f$ in this
case is the string scale $M_{s}$. Below the compactification scale, various
scales combine to give $f$ above KNP05  or below KimPLB17  the string scale,
depending on the compactification schemes. In this paper, we will not discuss
the ultraviolet completed theories because they depend on the details at the
Planck scale physics ConlonSven16 .

In the effective field theory approach, there are two classes for the explicit
breaking scale $\Lambda$: one class in the potential $V$ of spin-0 fields and
Yukawa couplings, and the other class the non-abelian gauge anomalies. Since
gravitational interactions always break any global symmetry Barr92 , we can
consider the explicit breaking term at $\Lambda$ due to gravity, of order
$\Lambda^{n+4}/M_{\rm P}^{n}$. For this term to be sufficiently small for QA
or ULA, theory must possess some symmetry to forbid terms up to the order
$\Lambda^{n+3}/M_{\rm P}^{n-1}$. In Ref. KimPRD16 , discrete symmetries were
used to realize this scheme. In this paper, we focus on the contribution of
non-abelian gauge anomalies such as the Peccei-Quinn(PQ) symmetry PQ77
U(1)${}_{\rm global}\times\rm SU(3)_{\rm gauge}^{2}$ and another global
symmetry U(1)${}_{\rm global}\times\rm SU(2)_{\rm gauge}^{2}$.

In Sec. II, we present the feasibility for the class of gauge anomaly,
U(1)${}_{\rm global}\times\rm SU(2)_{\rm gauge}^{2}$, to work for DE, and in
Subsec. II.2 the forbidden and allowed rare processes are discussed. In Sec.
III, gravity effects are briefly discussed. Section IV is a conclusion.

##  II Anomaly breaking

Non-abelian gauge forces break global symmetries Belavin75 . In the minimal
supersymmetric standard model (MSSM), the gauge coupling runs such that the
weak SU(2) anomaly may work for DE. We have the following SU(2) gauge coupling
at $M_{Z}$, using $\alpha_{\rm em}(M_{Z})=1/(127.940\pm 0.014)$,

$\begin{split}\frac{1}{\alpha_{2}(M_{Z})}=29.600\pm 0.010.\end{split}$ (2)

In the MSSM, there already exists an excellent estimate, based on the two-loop
MSSM running of $\alpha_{2}$ Bourilkov15 ,

Threshold correction [%] | $M_{SUSY}[\,\textrm{GeV}],$ | $Mg[\,\textrm{GeV}],$ | $1/\alpha_{GUT}$ | $\chi^{2}$
---|---|---|---|---
$+1$ | $10^{3.96\pm 0.10}$ | $10^{15.85\pm 0.03}$ | $26.74\pm 0.17$ | 8.2%
$\pm 0$ | $10^{3.45\pm 0.09}$ | $10^{16.02\pm 0.03}$ | $25.83\pm 0.16$ | 8.2%
$-1$ | $10^{3.02\pm 0.08}$ | $10^{16.16\pm 0.03}$ | $25.07\pm 0.15$ | 9.5%
$-2$ | $10^{2.78\pm 0.07}$ | $10^{16.25\pm 0.02}$ | $24.63\pm 0.13$ | 25.1%
$-3$ | $10^{2.60\pm 0.06}$ | $10^{16.31\pm 0.02}$ | $24.28\pm 0.10$ | 68.1%
$-4$ | $10^{2.42\pm 0.05}$ | $10^{16.38\pm 0.02}$ | $23.95\pm 0.09$ | 138.3%
$-5$ | $10^{2.26\pm 0.05}$ | $10^{16.44\pm 0.02}$ | $23.66\pm 0.09$ | 235.7%
Table 1: Combined fit results - CMS data CMS1 ; CMS2  and GUT unification. The
threshold correction $\varepsilon_{GUT}$ at ${M_{\rm GUT}}$ is defined as
$\alpha_{3}({M_{\rm GUT}})=\alpha_{GUT}(1+\varepsilon_{GUT})$. Here, without
the threshold correction the favored SUSY scale turns out to be 2820 GeV.

on the grand unification(GUT) scale as shown in Table 1. For SU($N$) groups,
we use the following $\beta$-function with gauge bosons and fermions Jones74 ;
Caldwell74 ,

$\begin{split}\beta=-(\frac{\alpha_{s}}{4\pi})^{2}\left(\frac{11}{3}C_{2}(G)-\frac{2}{3}\sum_{R}T(R)-\frac{\alpha_{s}}{4\pi}\Big{(}\frac{10}{3}\sum_{R}C_{2}(G)T(R)+2\sum_{R}{\rm
Cashimir}_{2}({\rm
SU}(N))T(R)-\frac{34}{3}(C_{2}(G))^{2}\Big{)}\right)\end{split}$ (3)

where $\alpha_{s}$ is the QCD coupling, and

$\begin{split}C_{2}({\rm SU}(N))=N,~{}{\rm Cashimir}_{2}({\rm
SU}(N))=\frac{N^{2}-1}{2N},~{}T(R)=\ell(R),\ell({\bf
N})=\frac{1}{2}.\end{split}$ (4)

Table 1 was obtained with the following input parameters,

$\begin{split}{\rm
At}~{}M_{Z}=91.19\,\,\textrm{GeV}:~{}\left\\{~{}\begin{array}[]{l}\sin^{2}\theta_{W}\Big{|}_{\overline{\rm
MS}}=0.23126\pm 0.00005,\\\\[3.00003pt] \alpha_{s}=0.1185\pm
0.0006,\end{array}.\right.\end{split}$ (5)

Without considering the threshold corrections at the GUT scale, we obtain
Bourilkov15 ,

$\displaystyle\textrm{MSSM}:$
$\displaystyle\left[\begin{array}[]{l}e^{-2\pi/\alpha_{2}}\Big{|}_{{M_{\rm
GUT}}}=1.69\times 10^{-81},\\\\[20.00003pt]
\left\\{~{}\begin{array}[]{l}M_{\rm
SUSY}=2820+670-540~{}\,\textrm{GeV},\\\\[3.00003pt] {M_{\rm GUT}}=(1.065\pm
0.06)\times 10^{16}\,\textrm{GeV}.\end{array}\right.\end{array}\right.$ (10)
$\displaystyle\textrm{SM}:$
$\displaystyle\left[\begin{array}[]{l}e^{-2\pi/\alpha_{2}}\Big{|}_{{M_{\rm
GUT}}}=1.69\times 10^{-131},\\\\[10.00002pt] {M_{\rm GUT}}=(1.096\pm
0.06)\times 10^{15}\,\textrm{GeV}.\end{array}\right.$ (13)

Note that the SU(2) coupling does not run with one-loop corrections in the
MSSM. In the SM, the SU(2) coupling runs and we obtain a much smaller factor
if not considering the threshold corrections. If the SU(2) gauge anomaly is
responsible for DE, we require the following order for $\Lambda$ in the MSSM,

$\begin{split}\textrm{MSSM}:~{}1.69\times
10^{-81}\Lambda^{4}=(0.003\,\,\mathrm{eV})^{4}\to\Lambda\sim 1.48\times
10^{8}\,\,\textrm{GeV},\\\ \textrm{SM}:~{}1.065\times
10^{-131}\Lambda^{4}=(0.003\,\,\mathrm{eV})^{4}\to\Lambda\sim 5.25\times
10^{20}\,\,\textrm{GeV}.\end{split}$ (14)

Thus, in the SM we cannot use the SU(2) gauge anomaly for the source of DE.
But, in the MSSM the SU(2) gauge anomaly can work for the explicit breaking
term for the QA. If it oscillated before the current epoch, then it cannot
work for DE but may work for a ULA KimPRD16 .

###  II.1 QCD axion plus QA/ULA

For DE, we need another global symmetry beyond the PQ symmetry, and the needed
discrete symmetry must be very restrictive such that many lower order terms
are not allowed. If we consider another U(1) global symmetry for DE, most
probably it has a QCD anomaly also. For it to be free of the QCD anomaly, its
quantum numbers of quarks must be chosen judiciously. In this section, we
comment on the SU(2) anomaly toward QA, and the possibility with
superpotential terms from string compactification will be interesting but
postponed to a future communication. In Eq. (14), $\Lambda$ for SU(2) anomaly
toward QA was calculated as $1.48\times 10^{8}\,\,\textrm{GeV}$. Note that,
among the terms breaking global symmetries the QCD anomaly is the largest one,
and QA/ULA must be chosen after removing the QCD anomaly part. Let the singlet
field housing the QCD axion and QA are $\sigma$ and $\sigma_{\rm
quint}\equiv\sigma_{q}$ whose VEVs are

$\begin{split}\langle\sigma\rangle=\frac{f_{a}}{\sqrt{2}}\,e^{ia/f_{a}},~{}\langle\sigma_{\rm
quint}\rangle=\frac{f_{q}}{\sqrt{2}}\,e^{ia_{q}/f_{q}}.\end{split}$ (15)

Let us introduce two complex scalar fields possessing the following quantum
numbers under U(1)PQ and U(1)q,

$\begin{split}\begin{array}[]{ccc}&~{}~{}\sigma{}{}&\sigma_{\rm quint}\\\ {\rm
U(1)_{PQ}:}&1&\Gamma_{2}\\\ {\rm
U(1)_{q}:}&\Gamma_{1}&1\end{array}\end{split}$ (16)

Then, the SU(3) and SU(2) instantons generate the following potential,

$\begin{split}V=m\Lambda^{3}_{\rm
QCD}\left(\cos(\frac{a}{f_{a}}+\Gamma_{2}\frac{a_{q}}{f_{q}})+{\rm
h.c.}\right)+f_{q}^{4}e^{-2\pi/\alpha_{2}}\left(\cos(\Gamma_{1}\frac{a}{f_{a}}+\frac{a_{q}}{f_{q}})+{\rm
h.c.}\right),\end{split}$ (17)

where $m$ is a typical chiral mass for light quarks, most probably very close
to $m_{u}$. The coefficient in the second term depends on details in the
models as shown in Yosh20 ; KimSpider09  and we simply denote it by the
parameter $f_{q}^{4}$. The mass matrix from Eq. (17) is

$\begin{split}\begin{pmatrix}\frac{m\Lambda_{\rm
QCD}^{3}+\Gamma_{1}^{2}f_{q}^{4}e^{-2\pi/\alpha_{2}}}{f_{a}^{2}},&\frac{\Gamma_{1}f_{q}^{4}e^{-2\pi/\alpha_{2}}+\Gamma_{2}m\Lambda_{\rm
QCD}^{3}}{f_{a}f_{q}}\\\\[6.99997pt]
\frac{\Gamma_{1}f_{q}^{4}e^{-2\pi/\alpha_{2}}+\Gamma_{2}m\Lambda_{\rm
QCD}^{3}}{f_{a}f_{q}},&\frac{f_{q}^{4}e^{-2\pi/\alpha_{2}}+\Gamma_{2}^{2}m\Lambda_{\rm
QCD}^{3}}{f_{q}^{2}}\end{pmatrix}.\end{split}$ (18)

The eigenvalues of (18) are

$\begin{split}m^{2}_{a,a_{q}}=&\frac{1}{2f_{a}^{2}f_{q}^{2}}\Biggl{(}f_{q}^{4}e^{\frac{-2\pi}{\alpha_{2}}}(f_{a}^{2}+f_{q}^{2}\Gamma_{1}^{2})+m\Lambda_{\rm
QCD}^{3}(f_{q}^{2}+f_{a}^{2}\Gamma_{2}^{2})\\\ \\\
&\mp\sqrt{-4(\Gamma_{1}\Gamma_{2}-1)^{2}\,m\Lambda_{\rm
QCD}^{3}{f_{a}^{2}}{f_{q}^{6}}e^{\frac{-2\pi}{\alpha_{2}}}+\left(f_{q}^{4}e^{\frac{-2\pi}{\alpha_{2}}}(f_{a}^{2}+\Gamma_{1}^{2}f_{q}^{2})+m\Lambda_{\rm
QCD}^{3}(f_{q}^{2}+f_{a}^{2}\Gamma_{2}^{2})\right)^{2}}\,\Biggr{)}.\end{split}$
(19)

If we expand in terms of $e^{-2\pi/\alpha_{2}}$ factor, which is very small,
two eigenvalues are

$\begin{split}m_{a}^{2}&=\frac{m\Lambda_{\rm
QCD}^{3}}{f_{a}^{2}}\Biggl{(}1+\frac{\Gamma_{2}^{2}f_{a}^{2}}{f_{q}^{2}}\Biggr{)}+O(e^{\frac{-2\pi}{\alpha_{2}}}),\\\
m_{a_{q}}^{2}&=(\Gamma_{1}\Gamma_{2}-1)^{2}f_{q}^{2}e^{\frac{-2\pi}{\alpha_{2}}}\frac{1}{1+\Gamma_{2}^{2}(f_{a}^{2}/f_{q}^{2})}+O(e^{\frac{-4\pi}{\alpha_{2}}})\,.\end{split}$
(20)

Note that $e^{-2\pi/\alpha_{2}}$ is almost $0.169\times 10^{-40}$. With
$f_{q}\simeq 10^{8}\,\textrm{GeV}$, we obtain $m_{q}\simeq 2\times
10^{-13}\,\textrm{GeV}\approx 0.0002\,\mathrm{eV}$ and the vacuum energy
density $f_{q}^{4}e^{-2\pi/\alpha_{2}}\approx(0.64\times
10^{-3}\,\mathrm{eV})^{4}$.

###  II.2 Breaking by potential $V$

Since we try to introduce explicit breaking terms of some U(1) global
symmetries, the conditions should come outside those global symmetries that we
considered previously, U(1)${}_{\rm global}\times\rm SU(3)_{\rm gauge}^{2}$
and U(1)${}_{\rm global}\times\rm SU(2)_{\rm gauge}^{2}$. Gravitational
interactions are known to break any global symmetry, e.g. at order $1/M_{\rm
P}^{n}$ as discussed in Barr92 ,

$\begin{split}\frac{\Lambda^{n+4}}{M_{\rm P}^{n}}.\end{split}$ (21)

For a given $\Lambda$, the lowest $n$ should be sufficiently large such that a
tiny but needed explicit breaking term arises. A possibility for Eq. (21) to
arise for a large $n$ can be a big discrete symmetry. Here we discuss in the
SUSY anti-SU(5) model because it reduces a number of Yukawa couplings of the
MSSM and still it is sufficient to present the essence in generating Eq. (21)
DKN84 ; Barr82 . In particular, it is also useful in the discussion of string
compactification LNP954 ; Kim20ijmpa . Matter fermions in anti-SU(5) are

$\begin{split}{\bf 10}_{g}(+1)=\Big{(}d^{c},q,N\Big{)}_{g},~{}{\overline{\bf
5}}_{g}(-3)=\Big{(}u^{c},\ell\Big{)}_{g},~{}{\bf
1}_{g}=(e^{+})_{g}(+5);~{}g=\\{1,2,3\\},\end{split}$ (22)

where quantum number $X=\\{1,-3,5\\}$ of anti-SU$(5)$ is shown in the
parantheses, $g$ is the generation index, and the SM doublets are $q$ in ${\bf
10}_{g}(+1)$ and $\ell$ in ${\overline{\bf 5}}_{g}(-3)$. The Higgs doublets
$H_{d}$ and $H_{u}$ of the MSSM are included in

$\begin{split}{\bf 5}_{H}(-2)\ni H_{d},~{}{\overline{\bf 5}}_{H}(+2)\ni
H_{u}.\end{split}$ (23)

We also need the GUT-breaking fields in the anti-SU(5),

$\begin{split}\Sigma\equiv{\bf
10}_{H}(+1)=\Big{(}D^{c},Q,N^{\prime}\Big{)}_{H},~{}\overline{\Sigma}\equiv{\overline{\bf
10}}_{H}(-1)=\Big{(}D,\bar{Q},\bar{N}^{\prime}\Big{)}_{H}\,.\end{split}$ (24)

We focussed on the anti-SU(5) GUT because it can be obtained at the level-1 in
the compactification of the heterotic string Kim20ijmpa ; LNP954 . With the
above anti-SU$(5)$ GUT representations, the Yukawa coupling for the top quark
mass is allowed at tree level, ${\bf 10}_{3}{\overline{\bf
5}}_{3}\,{\overline{\bf 5}}_{H}$. The Yukawa coupling for the down quark mass
is also allowed at tree level, ${\bf 10}_{3}{\bf 10}_{3}\,{\bf 5}_{H}$.
However, we must pay attention to the following rare process dimension-5
operators,

  * (i)

Forbid proton decay Sakai81 ; Hall83 :
$\\{q_{g}q_{g}q_{g}\ell_{g},d^{c}_{g}d^{c}_{g}u^{c}_{g}N_{g}\\}\in{\bf
10}\cdot{\bf 10}\cdot{\bf 10}\cdot{\overline{\bf
5}},~{}d^{c}_{g}u^{c}_{g}u^{c}_{g}e^{+}_{g}\in{\bf 10}\cdot{\overline{\bf
5}}\cdot{\overline{\bf 5}}\cdot{\bf 1}$.

  * (ii)

Allow neutrino masses Weinberg79 :
$N^{\prime}{N}^{\prime}\ell_{g}\ell_{g}H_{u}H_{u}\in{\bf 10}_{H}\cdot{\bf
10}_{H}\cdot{\overline{\bf 5}}\cdot{\overline{\bf 5}}\cdot{\overline{\bf
5}}_{H}\cdot{\overline{\bf 5}}_{H}$.

  * (iii)

Allow a TeV scale $\mu$ term KimNilles84 : $\mu\,H_{u}H_{d}$.

Item (i) has to be sufficiently suppressed but Items (ii) and (iii) must
arise. In Table 2, we considered a few discrete symmetries: matter parity
${\bf Z}_{2}$, R-parity ${\bf Z}_{2R}$ and R-parity ${\bf Z}_{4R}$. In the
Check column Y denotes the satisfied condition and N denotes the un-satisfied
condition. The item (iii) is not satisfied in all cases with the MSSM fields
plus singlets $N$ (in Eq. (22)), shown up to the ${\overline{\bf 5}}_{H}$
column, because the TeV scale of $\mu$ is not a natural one in the MSSM.

|  |  | NMSSM |  |  | anti-SU$(5)$ | PQ | Check
---|---|---|---|---|---|---|---|---
| ${\bf 10}_{g}$ |  ${\overline{\bf 5}}_{g}$ | ${\bf 1}_{g}$ | ${\bf 5}_{H}$ |  ${\overline{\bf 5}}_{H}$ |  $\Sigma_{\rm GUT}$ $\overline{\Sigma}_{\rm GUT}$ | $\sigma$ | Yes or No
${\bf Z}_{2}$ | $+1$ | $+1$ | $+1$ | $0$ | $0$ |  |  | (i) N, (ii) Y, (iii) N
${\bf Z}_{2R}$ | $+1$ | $+1$ | $+1$ | $0$ | $0$ |  |  | (i) N, (ii) Y, (iii) N
${\bf Z}_{4R}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+1$ | $+1$ |  |  | (i) Y, (ii) Y, (iii) N
${\bf Z}_{4R}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+1$ | $+1$ | $+4\quad-4$ | $+2$ | (i) Y, (ii) Y, (iii) Y
Table 2: Discrete symmetries with SUSY anti-SU(5) representations.

For item (iii) to be satisfied, a solution of the gauge hierarchy problem must
accompany the discrete symmetry. In this case, we must consider the GUT
breaking fields of Eq. (24), where
$\langle\Sigma\rangle=\langle\overline{\Sigma}\rangle=O(10^{17\,}\,\textrm{GeV})$
Ellis89 ; KimKyae07 . The column PQ $\sigma$ contains the singlet field
$\sigma$ whose VEV of O$(10^{10\,}\,\textrm{GeV})$ breaks the PQ symmetry to
create the invisible axion KSVZ1 . The GUT breaking fields allow the following
couplings Ellis89 ,

$\begin{split}{\bf 5}_{H}(-2)\Sigma(1)\Sigma(1)=\epsilon_{abdce}{\bf
5}^{a}{\bf 10}^{bc}{\bf 10}^{de},~{}{\overline{\bf
5}}_{H}(+2)\overline{\Sigma}(-1)\overline{\Sigma}(-1)=\epsilon^{abdce}{\bf
5}_{a}{\bf 10}_{bc}{\bf 10}_{de}.\end{split}$ (25)

So, the colored particle in ${\bf 5}$ is removed with the anti-colored
particle $D^{c}$ in ${\bf 10}$ by $\langle{\bf 10}^{45}\rangle$ at the GUT
scale. So does the anti-colored particle in ${\overline{\bf 5}}$ with the
colored particle $D$ in ${\overline{\bf 10}}$ by $\langle{\overline{\bf
10}}_{45}\rangle$. Then, we are left with a pair of Higgs doublets below the
scale ${M_{\rm GUT}}$, and we can consider the following superpotential,
generatibng a kind of $\mu$-term,

$\begin{split}{\bf 10}_{H}^{45}\,{\overline{\bf 5}}_{H\,4}\,{\overline{\bf
5}}_{H\,5},~{}\textrm{and }~{}{\overline{\bf 10}}^{H}_{45}\,{\overline{\bf
5}}^{H\,4}\,{\overline{\bf 5}}^{H\,5},\end{split}$ (26)

which, however, are not present if there is only one pair of the Higgs
doublets. If they are present (in case more than one pair of the Higgs
doublets are present), then together with the ${\bf 5}_{H}{\overline{\bf
5}}_{H}$ term, the following $\mu$ terms arise

$\begin{split}W_{\mu}=\mu_{1}H_{d}H_{u}+\mu_{2}H_{d}^{1}H_{d}^{2}+\mu_{3}H_{u}^{1}H_{u}^{2},\end{split}$
(27)

where $\mu_{1},\mu_{2},$ and $\mu_{3}$ are at the GUT scale, leading to the
following Higgsino mass matrix,

$\begin{split}~{}\hskip-14.22636pt&\quad
H_{d}^{1}\,H_{d}^{2}\,H_{u}^{1}\,H_{u}^{2}\\\ \begin{array}[]{c}H_{d}^{1}\\\
H_{d}^{2}\\\ H_{u}^{1}\\\ H_{u}^{2}\\\
\end{array}&\begin{pmatrix}0&\mu_{2}&0&\mu_{1}\\\ \mu_{2}&0&\mu_{1}&0\\\
0&\mu_{1}&0&\mu_{3}\\\ \mu_{1}&0&\mu_{3}&0\\\ \end{pmatrix}.\end{split}$ (28)

With the D-flat condition, leading to $\mu_{2}=\mu_{3}$, four eigenvalues turn
out to be

$\begin{split}\pm(\mu_{1}\pm\mu_{2}).\end{split}$ (29)

A pair of massless Higgsino arises in case $\mu_{2}=+\mu_{1}$ or $-\mu_{1}$.
Note, however, that $\mu_{2}=\mu_{3}=0$ in Eq. (27), as commented above, due
to the bosonic symmetry. At this level, therefore, we cannot derive a pair of
massless Higgs doublets. The only way to allow a massless pair is to have
$\mu_{1}=0$, i.e. there is no dimension-2 superpotential, as envisioned in
KimNilles84 ; Munoz93 . This case is shown in the last row of Table 2,
introducing an anti-SU(5) singlet field $\sigma$. Then, a nonrenormalizable
term allowed by the discrete symmetry is

$\begin{split}W\sim\frac{1}{M_{\rm P}}{\bf 5}_{H}{\overline{\bf
5}}_{H}\sigma\sigma+{\bf 5}_{H}(-2)\Sigma(1)\Sigma(1)+{\overline{\bf
5}}_{H}(+2)\overline{\Sigma}(-1)\overline{\Sigma}(-1)\end{split}$ (30)

The VEV, $f_{a}$, of $\sigma$ generates a $\mu=O(\frac{f_{a}^{2}}{M_{\rm
P}})$. If it breaks the PQ symmetry, then it must be at the intermediate scale
Kim84PLB  and the $\mu$ term is at the TeV scale as needed in the MSSM
phenomenology KimNilles84 . This $\sigma$ field may contribute to dark matter
in the Universe but not to DE.

##  III Gravity effects to effective potential

Discrete symmetries are the basis for the appearance of effective global
symmetries, which is manifest in the string compactification Kobayashi07 .

|  |  | NMSSM |  |  | anti-SU$(5)$ | PQ | Quintessential
---|---|---|---|---|---|---|---|---
| ${\bf 10}_{g}$ |  ${\overline{\bf 5}}_{g}$ | ${\bf 1}_{g}$ | ${\bf 5}_{H}$ |  ${\overline{\bf 5}}_{H}$ | $\Sigma_{\rm GUT},\overline{\Sigma}_{\rm GUT}$ | $\sigma$ | $\sigma_{q}$
${\bf Z}_{4R}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+1$ | $+1$ | $+4,\quad+4$ | $r_{1}=+2$ | $r_{2}=+2$
U(1)X | $+1$ | $-3$ | $+5$ | $-2$ | $+2$ | $+1,\quad-1$ | $0$ | $0$
U(1)PQ | $+1$ | $+1$ | $+1$ | $0$ | $0$ | $0,\quad~{}0$ | $+1$ | $\Gamma_{2}$
U(1)q | $+1$ | $+1$ | $+1$ | $0$ | $0$ | $0,\quad~{}0$ | $\Gamma_{1}$ | $+1$
Table 3: Quantum numbers of SUSY chiral fields. Heavy quarks in the KSVZ type
axion are not shown.

Let us illustrate with the ${\bf Z}_{4R}$ quantum numbers of the chiral fields
in the SUSY anti-SU(5) GUT of Table 2, shown again in the 2nd row of Table 3.
Let the ${\bf Z}_{4R}$ be a discrete subgroup of U(1)R. The integer quantum
numbers $r_{1}$ and $r_{2}$ are taken as the U(1)R quantum numbers. Then, an
effective superpotential term must carry the U(1)R quantum number $+2$ modulo
$4\textrm{ times integer}$.

Consider $\sigma$ first. For a term $\sigma^{n}$ to be present, $nr_{1}=2$
modulo $4\textrm{ times integer}$. For $r_{1}=1,2,3,4$, the allowed values of
$n$ are 2, none, 2, and none, respectively. So, it is appropriate to choose
$r_{1}=2$ or $4\textrm{ (times integer)}$ so that the PQ symmetry breaking is
soley through the U(1)PQ-SU(3)${}_{\rm color}^{2}$ anomaly. But, if we
consider ${\bf 5}_{H}{\overline{\bf 5}}_{H}$ also, we must choose $r_{1}=2$
such that ${\bf 5}_{H}{\overline{\bf 5}}_{H}\sigma$ is forbidden.

Next, consider $\sigma_{q}$. Again, to allow a term $\sigma_{q}^{n}$, one can
choose $nr_{2}=2$ or $4\textrm{ (times integer)}$ so that the weak-SU(2)
symmetry breaking is soley through the U(1)q-SU(2)${}_{W}^{2}$ anomaly. With a
minimal K$\ddot{\rm a}$hler potential term, a superpotential term
$W\sim\sigma_{q}^{n+3}/M_{\rm P}^{n}$ would allow the vacuum energy of order
$V\sim(n+3)|\langle\sigma_{q}\rangle|^{n+3}/M_{\rm P}^{n-1}$. If we take
$\langle\sigma_{q}\rangle=1.48\times 10^{8}\,\textrm{GeV}$ as implied by Eq.
(14), we must forbid the terms in the superpotential up to the order

$\begin{split}\frac{\sigma_{q}^{11.82}}{M_{\rm P}^{8.82}},\end{split}$ (31)

which needs the superpotential terms for $n\geq 9$. But, with the
aforementioned value for $n=2$ or 4, we cannot forbid terms with $n\leq 8$.
$\sigma_{q}$ alone cannot have a solution.

Therefore, let us consider the other fields of Table 3 also, especially ${\bf
5}_{H},{\overline{\bf 5}}_{H},\Sigma_{\rm GUT}$, and $\overline{\Sigma}_{\rm
GUT}$. If we choose $n=4$ (or 2) for $\sigma_{q}$, the superpotential ${\bf
5}_{H}{\overline{\bf 5}}_{H}\sigma_{q}$ is (or not) allowed. So, let us choose
$r_{2}=2$ for $\sigma_{q}$. In addition, inclusion of $\Sigma$ and
$\overline{\Sigma}$ is required not to allow the above result, and hence ${\bf
Z}_{4R}$ quantum numbers of $\Sigma$ and $\overline{\Sigma}$ are $\pm 4$.
Then, the GUT scale VEVs of these fields do not break the ${\bf Z}_{4R}$
symmetry, as required. These ${\bf Z}_{4R}$ quantum numbers are given in Table
3.

For $\Gamma_{1}\Gamma_{2}=1$, viz. Eq. (20), the QA mass is only through the
terms from the superpotential. The ${\bf Z}_{4R}$ allowed terms are

$\begin{split}\sigma^{2}\sigma_{q},~{}\sigma\sigma_{q}^{2},~{}\Sigma\overline{\Sigma}\sigma,~{}\Sigma\overline{\Sigma}\sigma_{q},~{}{\bf
5}_{H}{\overline{\bf 5}}_{H}\sigma^{2},{\bf 5}_{H}{\overline{\bf
5}}_{H}\sigma_{q}^{2},~{}\rm etc.\end{split}$ (32)

In particular, we need the following for the electroweak scale $\mu$-term
KimNilles84 ,

$\begin{split}\frac{1}{M_{\rm P}}{\bf 5}_{H}{\overline{\bf
5}}_{H}\sigma^{2}.\end{split}$ (33)

But the terms
$\sigma^{2}\sigma_{q},~{}\sigma\sigma_{q}^{2},~{}\Sigma\overline{\Sigma}\sigma,$
and $\Sigma\overline{\Sigma}\sigma_{q}$ should be forbidden. It can be
achieved by considering ${\bf Z}_{nR}\,(n\geq 9)$ or ${\bf Z}_{4R}\times{\bf
Z}_{\rm matter}.$ For ${\bf Z}_{4R}\times{\bf Z}_{2}$, the quantum numbers are
shown in Table 4.

|  |  | NMSSM |  |  | anti-SU$(5)$ | PQ | Quintessential
---|---|---|---|---|---|---|---|---
| ${\bf 10}_{g}$ |  ${\overline{\bf 5}}_{g}$ | ${\bf 1}_{g}$ | ${\bf 5}_{H}$ |  ${\overline{\bf 5}}_{H}$ | $\Sigma_{\rm GUT},\overline{\Sigma}_{\rm GUT}$ | $\sigma$ | $\sigma_{q}$
${\bf Z}_{4R}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+\frac{1}{2}$ | $+1$ | $+1$ | $+4,\quad+4$ | $r_{1}=+2$ | $r_{2}=+2$
${\bf Z}_{2}$ | $0$ | $0$ | $0$ | $+1$ | $+1$ | $+1,\quad+1$ | $+1$ | $+1$
Table 4: Working quantum numbers of SUSY chiral fields. Not only one but a few
heavy quarks (because the SM quarks carry PQ charges) are assumed to be added
to make the domain wall number one.

##  IV Conclusion

We discussed the possibility that the weak SU(2) gauge symmetry works for the
source of explicit breaking terms of a global symmetries suggested for a
quintessential axion. We find that it is possible in supersymmetric standard
models if the vacuum expectaion value of the quintessential axion field is
about $1.48\times 10^{8\,}\,\textrm{GeV}$.

###### Acknowledgements.

I have indebted from Junu Jeong and Myungbo Shim for checking the needed
numbers. This work is supported in part by the National Research Foundation
(NRF) grant NRF-2018R1A2A3074631.

## References

  * (1)
  * (2) Y. Nambu and G. Jona-Lasinio, _Dynamical model of elementary particles based on an analogy with superconductivity. I_ , Phys. Rev. 122 (1961) 345 [doi:10.1103/PhysRev. 122.345];
M. Gell-Mann, R. J. Oakes, and B. Renner, _Behavior of current divergences
under SU(3) $\times$SU(3)_, Phys. Rev. D 175 (1968) 2195
[doi:10.1103/PhysRev.175.2195].

  * (3) J. E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , Phys. Rev. Lett. 43 (1979) 103 [doi:10.1103/ PhysRevLett.43.103].
  * (4) For a review, see, J. E. Kim, S. Nam, and Y. Semertzidis, _Fate of global symmetries in the Universe: QCD axion, quintessential axion and trans-Planckian inflaton decay-constant_ , Int. J. Mod. Phys. A 33 (2018) 1830002 [e-print: 1712.08648 [hep-ph]].
  * (5) J. L. Goity, A. M. Bernstein, and B. R. Holstein, _Decay $\pi^{0}\to\gamma\gamma$ to next to leading order in chiral perturbation theory_, Phys. Rev. D 66 (2002) 076014 [e-print:hep-ph/0206007].
  * (6) R. E. Shrock and M. Voloshin, _Bounds on quark mixing angles from the decay $K_{L}\to\mu\bar{\mu}$_, Phys. Lett. B 87 (1979) 375 [doi:10.1016/0370-2693(79)90557-4].
  * (7) J. E. Kim and G. Carosi, _Axions and the strong CP problem_ , Rev. Mod. Phys. 82 (2010) 557 [e-print: 0807.3125 [hep-ph]].
  * (8) S. M. Carroll, _Quintessence and the rest of the world_ , Phys. Rev. Lett. 81 (1998) 3067 [e-Print: astro-ph/9806099 [astro-ph]].
  * (9) C. T. Hill, _Natural theories of ultralow mass pseudo Nambu-Goldstone bosons: Axions and quintessence_ , Phys. Rev. D 66 (2002) 075010 [e-print:hep-ph/0205237].
  * (10) J. E. Kim and H. P. Nilles, _A quintessential axion_ , Phys. Lett. B 553 (2003) 1 [e-print: hep-ph/0210402 [hep-ph]].
  * (11) J. E. Kim and D. J. E. Marsh, _An ultralight pseudoscalar boson_ , Phys. Rev. D 93 (2016) 025027 [e-Print: 1510.01701 [hep-ph]].
  * (12) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, _Ultralight scalars as cosmological dark matter_ , Phys. Rev. D 95 (2017) 043541 [e-print:1610.08297 [astro-ph.CO]].
  * (13) M. $\ddot{\rm O}$zer and M. O. Taha, _A model of the universe free of cosmological problems_ , Nucl. Phys. B 287 (1987) 776;
B. Ratra and P. J. E. Peebles, _Cosmological consequences of a rolling
homogeneous scalar field_ , Phys. Rev. D 37 (1988) 3406
[doi:10.1103/PhysRevD.37.3406];
C. Wetterich, _Cosmologies with variable Newton’s ’constant’_ , Nucl. Phys. B
302 (1988) 645 [doi:10.1016/0550-3213(88)90192-7];
H. Gies and C. Wetterich, _Renormalization flow from UV to IR degrees of
freedom_ , hep-ph/0205226;
J. A. Frieman, C. T. Hill, and R. Watkins, _Late-time cosmological phase
transitions: Particle-physics models and cosmic evolution_ , Phys. Rev. D 46
(1992) 1226 [doi:10.1103/PhysRevD.46.1226];
R. Caldwell, R. Dave, and P. J. Steinhardt, _Cosmological imprint of an energy
component with general equation of state_ , Phys. Rev. Lett. 80 (1998) 1582
[doi: 10.1103/PhysRevLett.80.1582];
P. Binetruy, _Models of dynamical supersymmetry breaking and quintessence_ ,
Phys. Rev. D 60 (1999) 063502 [e-Print: hep-ph/9810553 [hep-ph]];
C. Kolda and D. H. Lyth, _Quintessential difficulties_ , Phys. Lett. B 458
(1999) 197 [e-Print: hep-ph/9811375 [hep-ph]];
T. Chiba, _Quintessence, the gravitational constant, and gravity_ , Phys. Rev.
D 60 (1999) 083508 [doi: 10.1103/PhysRevD.60.083508];
Ph. Brax and J. Martin, _Quintessence and supergravity_ , Phys. Lett. B 468
(1999) 40 [e-Print:astro-ph/9905040];
A. Masiero, M. Pietroni, F. Rosati, _SUSY QCD and quintessence_ , Phys. Rev. D
61 (2000) 023504 [e-Print: hep-ph/9905346 [hep-ph]];
M. C. Bento and O. Bertolami, _Compactification, vacuum energy and
quintessence_ , Gen. Relativ. Gravit. 31 (1999) 1461 [e-Print: gr-qc/9905075
[gr-qc]];
F. Perrotta, C. Baccigalupi, S. Matarrase, _Extended quintessence_ , Phys.
Rev. D 61 (2000) 023507 [e-Print: astro-ph/9906066 [astro-ph]];
A. Arbey, J. Lesgourgues, and P. Salati, _Cosmological constraints on
quintessential halos_ , Phys. Rev. D 65 (2002) 083514 [e-Print: astro-
ph/0112324 [astro-ph]].

  * (14) J. E. Kim, H. P. Nilles, and M. Peloso, _Completing natural inflation_ , JCAP 01 (2005) 005 [e-Print: hep-ph/0409138 [hep-ph]].
  * (15) J. E. Kim, S. Nam, and B. Kyae, _The anomalous U(1) anom symmetry and flavors from an SU(5)$\times$SU(5)′ GUT in ${\bf Z}_{12-I}$ orbifold compactification_, Euro. Phys. J. C 77 (2017) 847 [e-Print: 1703.05345 [hep-ph]].
  * (16) See, for example, J. P. Conlon and S. Krippendorf, _Axion decay constants away from the lamppost_ , JHEP 04 (2016) 085 [e-print: 1601.00647 [hep-th]].
  * (17) Y. Nomura, T. Watari, and T. Yanagida, _Quintessence axion potential induced by electroweak instanton effects_ , Phys. Lett. B 484 (2020) 103 [e-Print:hep-ph/0004182].
  * (18) J. E. Kim, _Gaugino condensation scale of one family hidden SU(5) ′, dilaton stabilization and gravitino mass_, Phys. Lett. B 678 (2009) 197 [e-Print: 0904.1015 [hep-ph]].
  * (19) S. M. Barr and D. Seckel, _Planck-scale corrections to axion models_ , Phys. Rev. D 46 (1992) 539 [doi: 10.1103/PhysRevD.46.539];
M. Kamionkowski and J. March-Russell, _Are textures natural ?_ , Phys. Rev.
Lett. 69 (1992) 1485 [e-print: hep-th/9201063];
R. Holman, S. D. H. Hsu, E. W. Kolb, R. Watkins, and L. M. Widrow,
_Cosmological texture is sensitive to Planck-scale physics_ , Phys. Rev. Lett.
69 (1992) 1489 [doi:10.1103/PhysRevLett.69.1489];
R. Holman, S. D. H. Hsu, T. Kephart, E. W. Kolb, R. Watkins, and L. M. Widrow,
__, Phys. Lett. B 282 (1992) 132 [e-print: hep-th/9202003];
M. Kamionkowski and J. March-Russell, __, Phys. Lett. B 282 (1992) 137 [doi:
];
S. M. Lusignoli and M. Roncadelli, __, Phys. Lett. B 283 (1992) 278 [doi: ].

  * (20) R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , Phys. Rev. Lett. 38 (1977) 1440 [doi: 10.1103/PhysRevLett.38.1440].
  * (21) A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Yu. S. Tyupkin, _Pseudoparticle solutions of the Yang-Mills equations_ , Phys. Lett. B 59 (1975) 85 [doi: 10.1016/0370-2693(75)90163-X];
C. G. Callan, R. F. Dashen, and D, J. Gross, _The structure of the gauge
theory vacuum_ , Phys. Lett. B 63 (1976) 334 [doi:
10.1016/0370-2693(76)90277-X];
R. Jackiw and C. Rebbi, _Vacuum periodicity in a Yang-Mills quantum theory_ ,
Phys. Rev. Lett. 37 (1976) 172 [doi: 10.1103/PhysRevLett.37.172].

  * (22) D. Bourilkov, _Strong coupling running, gauge coupling unification and the scale of new physics_ , JHEP 11 (2015) 117 [e-print: 1508.04176 [hep-ph]];
D. R. T. Jones, _Asymptotic behaviour of supersymmetric Yang-Mills theories in
the two-loop approximation_ , Nucl. Phys. B 87 (1975) 127 [doi:
10.1016/0550-3213(75)90256-4].

  * (23) V. Khachatryan et al. (CMS Collaboration), _Constraints on parton distribution functions and extraction of the strong coupling constant from the inclusive jet cross section in pp collisions at sqrt(s) = 7 TeV_ , Euro. Phys. J. C 75 (2015) 288 [e-print: 1410.6765 [hep-ex]].
  * (24) V. Khachatryan et al. (CMS Collaboration), _Measurement of the inclusive 3-jet production differential cross section in proton-proton collisions at 7 TeV and determination of the strong coupling constant in the TeV range_ , Euro. Phys. J. C 75 (2015) 186 [e-print: 1412.1633 [hep-ex]].
  * (25) D. R. T. Jones, _Two loop diagrams in Yang-Mills theory_ , Nucl. Phys. B 75 (1974) 531 [doi: 10.1016/0550-3213(74)90093-5].
  * (26) W. E. Caldwell, _Asymptotic behavior of non-abelian gauge theories to two-loop order_ , Phys. Rev. Lett. 33 (1974) 244 [doi: 10.1103/PhysRevLett.33.244].
  * (27) S. M. Barr, _A new symmetry breaking pattern for SO(10) and proton decay_ , Phys. Lett. B 112 (1982) 219 [doi: 10.1016/0370-2693(82)90966-2].
  * (28) J. P. Derendinger, J. E. Kim, and D. V. Nanopoulos, _Anti-SU(5)_ , Phys. Lett. B 139 (1984) 170 [doi:10.1016/0370-2693(84)91238-3 ].
  * (29) K.-S. Choi and J. E. Kim, Quarks and Leptons from Orbifolded Superstring, 2nd ed., Lecture Notes in Physics Vol. 954 (Springer-Verlag, 2020).
  * (30) For a recent one, see, J. E. Kim, _Grand unfication models from SO(32) heterotic string_ , Int. J. Mod. Phys. A 35 (2020) 2050198 [ e-print: 2008.00367 [hep-th]], and references therein.
  * (31) N. Sakai and T. Yanagida, _Proton decay in a class of supersymmetric grand unified models_ , Nucl. Phys. B 197 (1981) 533 [doi: 10.1016/0550-3213(82)90457-6].
  * (32) L. J. Hall, J. D. Lykken, and S. Weinberg, _Supergravity as the messenger of supersymmetry breaking_ , Phys. Rev. D 27 (1983) 2359 [doi: 10.1103/PhysRevD.27.2359].
  * (33) S. Weinberg, _Baryon and lepton nonconserving processes_ , Phys. Rev. Lett. 43 (1979) 1566 [doi: 10.1103/PhysRevLett.43.1566].
  * (34) J. E. Kim and H. P. Nilles, _The $\mu$ problem and the strong CP problem_, Phys. Lett. B 138 (1984) 150 [doi:10.1016/0370-2693(84)91890-2].
  * (35) I. Antoniadis, J. R. Ellis, J. S. Hagelin, and D. V. Nanopoulos, _The flipped SU(5) x U(1) string model revamped_ , Phys. Lett. B 231 (1989) 65 [doi:10.1016/0370-2693(89)90115-9].
  * (36) J. E. Kim and B. Kyae, _Flipped SU(5) from ${\bf Z}_{(12-I)}$ orbifold with Wilson line_, Nucl. Phys. B 770 (2007) 47 [arXiv:hep-th/0608086].
  * (37) J. A. Casas and C. Munoz, _A Natural solution to the $\mu$ problem_, Phys. Lett. B 306 (1993) 288 [e-print:hep-ph/9302227 [hep-ph]].
  * (38) J. E. Kim, _A common scale for the invisible axion, local SUSY GUTs and saxino decay_ , Phys. Lett. B 136 (1984) 378 [doi:10.1016/0370-2693(84)92023-9].
  * (39) See, for example, T. Kobayashi, H. P. Nilles, F. Ploger, S. Raby, and M. Ratz, _Stringy origin of non-Abelian discrete flavor symmetries_ , Nucl. Phys. B 768 (2007) 135 [e-Print: hep-ph/0611020 [hep-ph]].


# Optimization of pencil beam scanning pattern for FLASH proton therapy

Sylvain Deffet  Molecular Imaging, Radiotherapy and Oncology, Institut de
Recherche Expérimentale et Clinique (IREC), Université catholique de Louvain,
1200 Woluwe-Saint-Lambert, Belgium  Edmond Sterpin  Molecular Imaging,
Radiotherapy and Oncology, Institut de Recherche Expérimentale et Clinique
(IREC), Université catholique de Louvain, 1200 Woluwe-Saint-Lambert, Belgium
Particle Therapy Interuniversity Center Leuven - PARTICLE, Leuven, 3000,
Belgium

###### Abstract

Background: The FLASH effect, which reduces the radiosensitivity of healthy
tissue while maintaining tumor control at high dose rates, has shown potential
for improving radiation therapy. While the mechanisms behind the effect are
not fully understood, it has been extensively studied with MeV electron beams
and high-energy proton beams. However, to achieve FLASH proton therapy,
changes to equipment and delivery systems are needed. Conformal FLASH proton
therapy involves advanced beam-shaping technologies and specialized nozzle
designs to confine the dose to the target volume. Optimizing the spot delivery
pattern and range modulators can enhance the local dose rate, and genetic
algorithms have been used to optimize scan patterns for stereotactic FLASH
proton therapy of early-stage lung cancer and lung metastases.
Purpose: Maximize the dose rate within regions of interest through an
efficient approach grounded in graph theory.
Methods: We have created a graph-based algorithm to optimize the trajectory of
proton spots to maximize the 100th percentile dose rate. Since this problem is
NP-hard, we have employed an approximation algorithm that can solve this kind
of Traveling Salesman Problem (TSP) efficiently.
Results: When compared to a conventional serpentine pattern, the optimized
scanning trajectory led to a doubling of the median dose rate, but only a
minor increase in DR95. Our approach is more efficient and requires fewer
evaluations of the objective function and hyper-parameters compared to
existing genetic algorithms.
Conclusions: The optimized scanning trajectory led to a doubling of the median
dose rate, but only a minor increase in DR95. The extent to which the dose
rate can be increased depends on the size and shape of the region of interest.
Future research could explore integrating FLASH objectives into treatment
planning and incorporating the proposed method into plan optimization.

_Keywords_ FLASH, Flash Proton Therapy, Dose Rate, PBS

##  1 Introduction

The FLASH effect, which refers to a significant reduction of the
radiosensitivity of healthy tissue while maintaining tumor control at
ultrahigh dose rates, has garnered considerable interest in the radiation
therapy community since it was reported in 2014[4]. Although FLASH has a huge
potential to increasing the therapeutic bandwidth of radiation therapy, the
mechanisms behind the effect are still not fully understood, but possible
explanations focus on the role of oxygen[5], radiochemistry, and the immune
system[14, 18, 12, 9, 8, 5]. The effect has been most extensively studied in
radiobiological experiments with MeV electron beams, and a first patient has
been successfully treated with FLASH electron beam therapy. Recent research
has demonstrated that the FLASH effect is also present in high-energy proton
beams[1], which may be particularly suited for deep-seated targets.

One of the main challenges in developing FLASH proton therapy is the need to
significantly increase the dose rate. Preclinical experiments are typically
conducted using small fields covered by a passive scattering method, but the
maximum achievable field size for a given dose rate is directly limited by the
maximum current that the system can output. However, the pencil beam scanning
(PBS) technique has the potential to locally achieve high dose rates due to
the limited size of the spots that are delivered individually.

To achieve a FLASH dose rate in PBS, a number of changes to the equipment and
delivery system used in intensity modulated proton therapy (IMPT) are needed.
In particular, the proton beam must be delivered at a much higher energy than
is generally used in IMPT to ensure a high transmission efficiency[10].
Additionally, using a treatment plan that involves multiple energies incurs
delays required by the system to switch from one energy to the next,
negatively impacting the average dose rate.

One approach to delivering a FLASH treatment involves shooting at maximum
energy through the patient with so-called transmission beams, but this results
in a significant amount of dose being delivered after the tumor, compromising
the superior dosimetric potential of protons[15].

Conformal FLASH proton therapy, on the other hand, involves the use of
advanced beam-shaping technologies and specialized nozzle designs to confine
the dose to the target volume, similar to IMPT. A patient-specific range
modulator, located between the nozzle and the patient, is used to tailor the
range of the proton beam, enabling a conformal treatment plan with a single
high-energy layer. Several methods have been proposed to optimize range
modulators[16, 13, 20, 3].

In addition to utilizing high energy beams, optimizing the spot delivery
pattern can enhance the local dose rate[11]. This optimization is a
combinatorial problem that requires the use of approximation algorithms when
the dimensionality of the problem increases. For instance, José Santo et
al.[11] applied genetic algorithms to optimize scan patterns for stereotactic
FLASH proton therapy of early-stage lung cancer and lung metastases. These
algorithms can easily adapt to complex cost functions such as maximizing the
dose rate in specific organs at risk. However, there are costs associated with
these methods, such as the need to tune hyper-parameters and the execution
time required.

We propose a fast and effective method to optimize the dose rate in specific
regions of interest (ROIs) based on graph theory. In an in silico study, the
optimized pattern is then compared with the conventional serpentine pattern
for a head and neck case in terms of computed dose rates.

##  2 Materials and Methods

###  2.1 Dose rate definition

In order to optimize the dose rate, an essential first step is to establish a
clear definition. In PBS proton therapy, local variations in dose rate occur
as each voxel receives dose contributions from nearby PBS spots. Therefore, it
is essential to consider the irradiation time of the spots and the time taken
for the pencil beam to move from one position to the next. These factors have
been incorporated into the definition of the PBS dose rate which has been
established by Folkerts et al.[7, 6]. This definition was later extended to
explicitely include a dose threshold expressed as a percentage of the dose
delivered to the voxel, resulting in the following percentile dose rate:[2]:

$\accentset{\mbox{.}}{D}^{P}_{i}=\frac{pD_{i}}{t_{1,i}-t_{0,i}}$ (1)

where $t_{1,i}=t_{i}(\frac{(1-p)}{2}2D_{i})$,
$t_{0,i}=t_{i}(\frac{(1+p)}{2}D_{i})$, $p$ is an arbitrary percentage, and
$D_{i}$ it the total accumulated dose received in voxel $i$.

The maximum percentile dose rate proposed by Deffet et al.[2] extends the
concept of percentile dose rate by considering all time windows in which the
accumulated dose is at least $pD_{i}$:

$\displaystyle\accentset{\mbox{.}}{D}^{p}_{i}$ $\displaystyle=$
$\displaystyle\max_{t_{0},t_{1}}\frac{\int_{t_{0}}^{t_{1}}d_{i}(t)dt}{t_{1}-t_{0}}$
$\displaystyle s.t.$ $\displaystyle\int_{t_{0}}^{t_{1}}d_{i}(t)dt\geq pD_{i}$
$\displaystyle t_{1}>t_{0}$

where $d_{i}(t)$ is the dose received in voxel $i$ at time $t$.

In the present study, we focus on optimizing the 100-percentile dose rate for
each voxel, which is the dose delivered to the voxel divided by the
corresponding time interval over which it is delivered. Our selection of the
100-percentile dose rate as the optimization objective is the result of the
graph representation that will be introduced later in the paper, and which
will be thoroughly discussed.

According to Eq. 1, the 100-percentile dose rate is:

$\accentset{\mbox{.}}{D}^{100}_{i}(\mathbf{S})=\frac{D_{i}}{t_{1}(\mathbf{S})-t_{0}(\mathbf{S})}=\frac{D_{i}}{T_{i}^{100}(\mathbf{S})}\\\
$ (3)

where $\mathbf{S}$ is the spot sequence, $t_{1}(\mathbf{S})=t_{i}(D_{i}^{-})$
and $t_{0}(\mathbf{S})=t_{i}(0^{+})$ and where we introduce
$T_{i}^{100}(\mathbf{S})=t_{1}(\mathbf{S})-t_{0}(\mathbf{S})$ as the quantity
will be used many times later.

$\accentset{\mbox{.}}{D}^{100}_{i}(\mathbf{S})$ is equivalent to the PBS and
percentile dose rates with a dose threshold of $0^{+}~{}\mathrm{Gy}$. The use
of a threshold of $0^{+}~{}\mathrm{Gy}$ instead of $0~{}\mathrm{Gy}$ is the
mathematical expression that we do not want to count the time when no dose is
given before the first dose contribution to the voxel and also after the last
contribution to the voxel.

###  2.2 Optimization problem

The objective is to maximize the dose rate in some specific ROIs, ie. to
solve:

$\arg\max_{\mathbf{S}}\sum_{i\in
ROI}\accentset{\mbox{.}}{D}^{100}_{i}(\mathbf{S})$ (4)

We relax this objective and rather consider the minimization of the time
required to deliver the dose to each voxel:

$\arg\min_{\mathbf{S}}\sum_{i\in ROI}T_{i}^{100}(\mathbf{S})$ (5)

In other words, we are going to determine the spot sequence that minimizes the
averaged time required to deliver the dose to each voxel. One very important
point is that for each voxel, we do not count the time when no dose is given
before any dose contribution to the voxel and also after that the voxel has
received its full dose.

####  2.2.1 Optimization of the delivery time of the whole field

Our initial focus is on the fundamentals of graph theory as applied to the
optimization of treatment field delivery time. While this is not our primary
objective, this exercise is helpful in laying the groundwork for many of the
definitions and concepts that will be utilized later on.

All the possibles delivery timing can be represented by means of a fully
connected graph, named $G$. The corresponding adjacency matrix is called
$\mathbf{M}$. $M_{s_{1},s_{2}}$ is the time, also named $T_{s_{1}\rightarrow
s_{2}}$, required to move from spot $s_{1}$ to spot $s_{2}$. $M_{s_{1},s_{1}}$
is the time, also named $T_{s_{1}}$, is the irradiation time associated to
spot $s_{1}$:

$\mathbf{M}=\begin{bmatrix}T_{s_{1}}&T_{s_{1}\rightarrow
s_{2}}&T_{s_{1}\rightarrow s_{3}}&\dots\\\ T_{s_{2}\rightarrow
s_{1}}&T_{s_{2}}&T_{s_{2}\rightarrow s_{3}}&\dots\\\
\vdots&\vdots&\ddots&\vdots\\\ \dots&\dots&\dots&\dots\end{bmatrix}$

According to our simple delivery model, $M$ is symmetrical, i.e.
$T_{s_{1}\rightarrow s_{2}}=T_{s_{2}\rightarrow s_{1}}$ for all $s_{1},s_{2}$.

In the context of optimizing the time required to deliver a treatment field,
the eligible delivery sequences refer to all the possible routes that visit
each spot exactly once. If a delivery sequence is feasible, the delivery time
of the whole treatment field can be represented by the integrated distance
along the sequence. Thus, we can optimize the time required to deliver the
whole field by finding the shortest path in $G(\mathbf{M})$. This can be
viewed as an instance of the well-known Travelling Salesman Problem (TSP), for
which numerous algorithms exist in the literature, and which can be selected
depending on the scale of the problem at hand.

An equivalent way to build such a graph is not to place the delivey times
$T_{s}$ on the vertices but to add them on the edges. The corresponding
adjancency matrix is thus:

$\mathbf{M}=\left[\begin{array}[]{c|ccc}0&T_{s_{1}}&T_{s_{2}}&\dots\\\
\hline\cr 0&0&T_{s_{1}\rightarrow s_{2}}+T_{s_{2}}&\dots\\\
0&T_{s_{2}\rightarrow s_{1}}+T_{s_{1}}&0\dots\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&\dots&\dots&\dots\end{array}\right]$ (6)

This matrix has zeros everywhere on the diagonal and is not symmetrical
anymore. This matrix has one more row and colums with respect to the previous
one. The first row ensures that the beam-on time of the first spot to be
delivered will be accounted for. The first column set to zero everywhere
ensures that the first row can only be reached once.

In this representation, a notable advantage is that considering the shortest
path on the corresponding graph results in the consideration of all possible
starting spots simultaneously, due to the addition of an extra, fictional spot
$s_{0}$. However, if a specific starting point is desired, the first row and
column can be removed, and the matrix can be reorganized such that the
starting spot occupies the first row and column.

####  2.2.2 Optimization of the local dose rate

In light of the above method for optimizing the delivery timing of the entire
treatment field, we can now focus on the specific task of optimizing the sum
of $T_{i}^{100}$ across all voxels within a designated region of interest
(ROI).

Considering a sequence of spots $\mathbf{S}$ which is a path $P$ in the graph,
the sum of $T_{i}^{100}$ can be obtained by computing, for each voxel, the
minimum spanning subgraph in $G$ which spans over the spots which contribute
to $i$:

$\displaystyle\sum_{i}T_{i}^{100}(P,\mathbf{M})$ $\displaystyle=$
$\displaystyle\sum_{i}P^{i}(\mathbf{M})$ (7)

where $P^{i}$ is the minimal subpath in $P$ that delivers the full dose to
voxel $i$, and $P^{i}(\mathbf{M})$ is the length of $P^{i}$ computed on
$G(\mathbf{M})$.

To determine the optimal sequence we could use brute force:

  1. 1.

Compute every feasible paths $P$;

  2. 2.

For each feasible path:

    1. (a)

For each voxel $i$, contribute the sets of spots which contribute to the dose
of the voxel;

    2. (b)

Find the minimum spanning supbath $P^{sub}\in P$ which spans over all the
spots which contribute to $i$ and compute its cumulative length $T_{i}$

    3. (c)

Add the length of this subpath to the cumulative sum $\sum_{i}T_{i,p}^{100}$

The optimal path is the one that minimizes $\sum_{i}T_{i,p}^{100}$ over all
possible paths $p$.

Given the infeasibility of using brute force to optimize problems with a
significant amount of delivery spots, we propose a modified approach to graph
optimization. To prioritize routes that contribute to the largest number of
voxels which should decrease $\sum_{i}T_{i}^{100}$, we adjust the edge weights
of the original graph used for optimizing the entire treatment field.
Specifically, we construct an adjacency matrix $\mathbf{E}$, where $E_{s,s}=0$
for all $s$, and $E_{s_{1},s_{2}}$ represents the ratio of (1) the time
required to deliver a dose from spot $s_{2}$ when starting from spot $s_{1}$,
and (2) the number of voxels in the ROI that receive a dose contribution from
spot $s_{2}$.

For the sake of conciseness, we define the number of voxels that receive
contributions from spots $s$ as follows.

$N_{s}=\sum_{i}\delta_{D_{i,s}>0}$ (9)

We now build the adjacency matrix $\mathbf{E}$ similarly to $\mathbf{M}$ but
with

$\displaystyle E_{s_{1},s_{2}}$ $\displaystyle=$
$\displaystyle\frac{T_{s_{1}\rightarrow s_{2}}+T_{s_{2}}}{N_{s_{2}}}$ (10)
$\displaystyle=$ $\displaystyle\frac{M_{s_{1},s_{2}}}{N_{s_{2}}}$ (11)

The adjacency matrix $\mathbf{E}$ is thus:

$\mathbf{E}=\mathbf{M}\left[\begin{array}[]{c|cccc}1&0&0&0&\dots\\\ \hline\cr
0&\frac{1}{N_{s_{1}}}&0&0&\dots\\\ 0&0&\frac{1}{N_{s_{2}}}&0&\dots\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&\dots&\dots&\dots&\dots\end{array}\right]$ (13)

The proposed approach aims to favor transitions based on the amount of voxels
that receive dose at the end of the transition.

The resulting algorithm for finding the optimal spot sequence is:

  1. 1.

Select the subsets $S$ of spots which contribute to the ROI

  2. 2.

Compute $T_{s_{i}}\forall s_{i}\in S$

  3. 3.

Compute $T_{s_{i}\rightarrow s_{j}}\forall s_{i},s_{j}\in S$

  4. 4.

Compute dose influence matrix $D_{i,s}\forall i\in ROI,s\in S$

  5. 5.

From $D_{i,s}$, compute $N_{s_{j}}\forall s_{j}\in S$. The utilization of a
threshold may be implemented to eliminate any contributions to the dose that
are deemed not significant enough to warrant inclusion in the calculation.

  6. 6.

Compute $\mathbf{E}$ from $T_{s_{i}},T_{s_{i}\rightarrow s_{j}},N_{s_{j}}$

  7. 7.

Solve optimal_ordering = TSP($\mathbf{E}$)

###  2.3 In silico assessment

We applied spot pattern optimization on conformal FLASH treatment plans
calculated using the methodology presented in our prior publication[3]. In
this paper, a treatment plan was optimized on a head and neck case which is
reused in the present publication. The PTV considered in the present study had
a prescription of 54 Gy. However, dose rates must be computed per fraction. As
FLASH treatments will most likely be hypofractionated[19], we considered that
the dose per fraction would be around 8 Gy.

The treatment plan and its associated range modulator were computed for an
energy of 226 MeV, a spot size of
$(\sigma_{x}=4.5~{}\mathrm{mm},\sigma_{y}=5~{}\mathrm{mm})$. As the spot
spacing has a direct impact on the scanning time and thus on the average dose
rate, we considered two spot spacings: 4 mm which is close to that
conventionnaly used, and 15 mm which is expected to yield higher dose rates.
It is to be noted that because of the higher scattering introduced by the
passive degradation of the beam, larger spot spacing may be used than in IMPT
without introducing significant degradation of the dose[2].

To facilitate the comparison of the different dose rate formulas, we used a
simple model where:

  1. 1.

the nozzle output current was considered constant;

  2. 2.

the time between 2 spots was proportional to the distance between the spots.

In other words,

  1. 1.

the irradiation time of a spot was the ratio between the charge of the spot
and the current;

  2. 2.

the time separating 2 spots was the ratio between the distance of the spots
and an assumed scanning speed.

The current was 500 nA (averaged on a pulse period) at the output of the
nozzle and the scanning speed was 8000 mm/s. We considered that 1 MU
corresponds to 152,880,000 protons.

A dose influence matrix was calculated using MCsquare[17] to determine the
contributions of each spot to the dose of each voxel.

In order to showcase the potential of the proposed approach, we applied it to
optimize the spot pattern on both regions of interest (ROIs) corresponding to
organs at risk as well as an extention of the PTV. The resulting dose rate
maps were then compared to that obtained with the conventional serpentine
pattern.

##  3 Results

In our study, the submandibular (brown contour) and parotid (red contour)
glands in the head and neck region were used as examples of regions of
interest in which we could optimize the dose rate. The submandibular gland,
being proximal to the target volume, receives a significant dose, while the
smaller parotid gland is mostly within the target volume. First, we optimized
the dose rate exclusively in the submandibular gland using our algorithm.
Then, we optimized the dose rate exclusively in the parotid. Finally, we
sought to maximize the dose rate in both the parotid and submandibular glands
using the algorithm. For each spot scanning pattern, we computed the
95-percentile dose rates and the corresponding DRVHs. The results are
presented in Fig. 2 and Fig. 1 for a spot spacing of 15 mm. The dashed DRVHs
represent the serpentine pattern and the solid DRVHs correspond to the
optimized spot delivery pattern. The yellow DRVH in the figure corresponds to
the PTV.

Figure 1: (a) Unoptimized spot pattern and (b) spot pattern optimized to
maximize the dose rate in both the parotid and submadibular glands with a spot
spacing of 15 mm. Figure 2: For a spot spacing of 15 mm, (a) 95-percentile
dose rate for the unoptimized (serpentine) spot pattern and (b) for the spot
pattern optimized to maximize the dose rate in both the parotid and
submadibular glands. (c) 95-percentile dose rate volume histograms for
unoptimized and optimized spot patterns.

The results show a much higher DR50 with the optimized scanning pattern than
with the conventional serpentine pattern. However, the DR95 is only slighty
improved. In addtion, we see that the improvement is much higher when only one
organ at risk is considered. Overall, this means that the improvement of the
dose rate values is limited by the size and the complexity of the volume in
which we want to maximize it, which is not suprising.

We conducted the same optimization procedure on a treatment plan with a spot
spacing of 4 mm. Fig. 4 shows that the resulting dose rate was lower than with
a spot spacing of 15 mm which is due to the accumulated scanning time being
larger as one can figure out by comparing spot maps in Fig. 1 and Fig. 3.
After optimizing the spot pattern, we found that a DR50 greater than 40 Gy/s
was achievable in either the parotid or the submandibular gland, but not both
simultaneously. These results demonstrate the crucial role that spot spacing
plays in determining the maximum achievable dose rate.

Figure 3: (a) Unoptimized spot pattern and (b) spot pattern optimized to
maximize the dose rate in both the parotid and submadibular glands with a spot
spacing of 4 mm. Figure 4: For a spot spacing of 4 mm, (a) 95-percentile dose
rate for the unoptimized (serpentine) spot pattern and (b) for the spot
pattern optimized to maximize the dose rate in both the parotid and
submadibular glands. (c) 95-percentile dose rate volume histograms for
unoptimized and optimized spot patterns.

##  4 Discussion

In PBS proton therapy, local variations in dose rate occur as each voxel
receives dose contributions from nearby PBS spots. In order to accurately
calculate the dose rate, it is essential to consider the irradiation time of a
spot and the time taken for the pencil beam to move from one position to the
next. These factors have been incorporated into the definition of the PBS dose
rate, which we refer to as the percentile dose rate when the dose threshold is
expressed as a percentage of the dose delivered to the voxel.

In our study, we have created a graph-based algorithm that optimizes the spot
trajectory with the goal of maximizing the 100-percentile dose rate. It is
widely acknowledged that solving such a graph-based problem is more efficient
than utilizing genetic algorithms and involves fewer hyper-parameters to be
fine-tuned. However, owing to the NP-hard characteristic of the TSP problem,
approximation algorithms are typically necessary. In addition, it should be
noted that genetic algorithms possess certain advantages too, such as
modularity and the capability to optimize various definitions of the dose
rate. For instance, the genetic algorithm proposed by José Santos et al. could
optimize both the 95- and 100-percentile dose rates. In contrast, it is not
feasible to directly optimize the 95-percentile dose rate using our proposed
method. Nonetheless, the 100-percentile dose rate can be considered the upper
limit of the percentile dose rate, and optimizing it is expected to lead to an
improvement in the 95-percentile dose rate too, as was observed in our in
silico study.

In comparison to the serpentine pattern, our optimized scanning trajectory
achieves a two-fold increase in the median dose rate. However, two important
caveats should be noted. Firstly, the increase in DR95 is significantly lower
than that of the median dose rate. Secondly, the degree to which the dose rate
can be increased depends on the size and shape of the ROI being considered. It
is unrealistic to expect that the simple optimization of the scanning
trajectory can double the dose rate accross the entire CTV and its extension.
Nevertheless, our results are promising with respect to organs at risk, which
we prioritize to benefit from the FLASH effect. Integrating FLASH objectives
into treatment planning and incorporating our proposed method into the
optimization of the plan are potential avenues for future research.

##  5 Conclusions

A graph-based algorithm has been developed to optimize spot trajectory to
maximize the dose rate in ROIs. The optimized scanning trajectory achieved a
two-fold increase in the median dose rate but only a limited increase in DR95.
The degree to which the dose rate can be increased depends on the size and
shape of the region of interest, and integrating FLASH objectives into
treatment planning and incorporating the proposed method into plan
optimization are potential future research avenues.

##  6 Acknowledgments

This work was supported by the Walloon Region of Belgium through technology
innovation partnership no. 8341 (EPT-1 – Emerging Proton Therapies Phase 1)
co-led by MecaTech and BioWin clusters.

##  7 Conflict of Interest Statement

We have no conflicts of interest to disclose.

## References

  * [1] Shannon Cunningham, Shelby McCauley, Kanimozhi Vairamani, Joseph Speth, Swati Girdhani, Eric Abel, Ricky A. Sharma, John P. Perentesis, Susanne I. Wells, Anthony Mascia, and Mathieu Sertorio.  Flash proton pencil beam scanning irradiation minimizes radiation-induced leg contracture and skin toxicity in mice.  Cancers, 13(5), 2021.
  * [2] Sylvain Deffet, Valentin Hamaide, and Edmond Sterpin.  Definition of flash dose rate in pencil beam scanning: a comparative study.  Medical Physics, submitted.
  * [3] Sylvain Deffet, Kevin Souris, and Edmond Sterpin.  Optimization of patient-specific range modulators for conformal flash proton therapy.  Physica Medica, submitted.
  * [4] Vincent Favaudon, Laura Caplier, Virginie Monceau, Frédéric Pouzoulet, Mano Sayarath, Charles Fouillade, Marie-France Poupon, Isabel Brito, Philippe Hupé, Jean Bourhis, Janet Hall, Jean-Jacques Fontaine, and Marie-Catherine Vozenin.  Ultrahigh dose-rate flash irradiation increases the differential response between normal and tumor tissue in mice.  Science Translational Medicine, 6(245):245ra93–245ra93, 2014.
  * [5] Vincent Favaudon, Rudi Labarbe, and Charles L. Limoli.  Model studies of the role of oxygen in the flash effect.  Medical Physics, 49(3):2068–2081, 2022.
  * [6] Michael Folkerts, Eric Abel, Simon Busold, Jessica PEREZ, Vidhya Krishnamurthi, and C. Clifton LING.  Defining dose rate for pencil beam scanning.  US patent, (US20210393981A1), 2020.
  * [7] Michael M. Folkerts, Eric Abel, Simon Busold, Jessica Rika Perez, Vidhya Krishnamurthi, and C. Clifton Ling.  A framework for defining flash dose rate for pencil beam scanning.  Medical Physics, 47(12):6396–6404, 2020.
  * [8] Anna A. Friedl, Kevin M. Prise, Karl T. Butterworth, Pierre Montay-Gruel, and Vincent Favaudon.  Radiobiology of the flash effect.  Medical Physics, 49(3):1993–2013, 2022.
  * [9] Jian-Yue Jin, Anxin Gu, Weili Wang, Nancy L. Oleinick, Mitchell Machtay, and Feng-Ming (Spring) Kong.  Ultra-high dose rate effect on circulating immune cells: A potential mechanism for flash effect?  Radiotherapy and Oncology, 149:55–62, 2020.
  * [10] Simon Jolly, Hywel Owen, Marco Schippers, and Carsten Welsch.  Technical challenges for flash proton therapy.  Physica Medica, 78:71–82, 2020.
  * [11] Rodrigo José Santo, Steven J.M. Habraken, Sebastiaan Breedveld, and Mischa S. Hoogeman.  Pencil-beam delivery pattern optimization increases dose rate for stereotactic flash proton therapy.  International Journal of Radiation Oncology*Biology*Physics, 2022\.
  * [12] Rudi Labarbe, Lucian Hotoiu, Julie Barbier, and Vincent Favaudon.  A physicochemical model of reaction kinetics supports peroxyl radical recombination as the main determinant of the flash effect.  Radiotherapy and Oncology, 153:303–310, 2020.  Physics Special Issue: ESTRO Physics Research Workshops on Science in Development.
  * [13] Ruirui Liu, Serdar Charyyev, Niklas Wahl, Wei Liu, Minglei Kang, Jun Zhou, Xiaofeng Yang, Filipa Baltazar, Martina Palkowitsch, Kristin Higgins, William Dynan, Jeffrey Bradley, and Liyong Lin.  An integrated biological optimization framework for proton sbrt flash treatment planning allows dose, dose rate, and let optimization using patient-specific ridge filters, 2022.
  * [14] Guillem Pratx and Daniel S Kapp.  A computational model of radiolytic oxygen depletion during FLASH irradiation and its effect on the oxygen enhancement ratio.  Physics in Medicine & Biology, 64(18):185005, sep 2019.
  * [15] Bethany Rothwell, Matthew Lowe, Erik Traneus, Miriam Krieger, and Jan Schuemann.  Treatment planning considerations for the development of flash proton therapy.  Radiotherapy and Oncology, 175:222–230, Oct 2022.
  * [16] Yuri Simeonov, Uli Weber, Christoph Schuy, Rita Engenhart-Cabillic, Petar Penchev, Veronika Flatten, and Klemens Zink.  Development, monte carlo simulations and experimental evaluation of a 3d range-modulator for a complex target in scanned proton therapy.  Biomedical Physics & Engineering Express, 8(3):035006, mar 2022\.
  * [17] Kevin Souris, J Lee, and E Sterpin.  Fast multipurpose monte carlo simulation for proton therapy using multi- and many-core cpu architectures.  Medical Physics, 43(4):1700–1712, 2016.
  * [18] Douglas R. Spitz, Garry R. Buettner, Michael S. Petronek, Joël J. St-Aubin, Ryan T. Flynn, Timothy J. Waldron, and Charles L. Limoli.  An integrated physico-chemical approach for explaining the differential impact of flash versus conventional dose rate irradiation on cancer and normal tissue responses.  Radiotherapy and Oncology, 139:23–27, Oct 2019.
  * [19] Steven van de Water, Sairos Safai, Jacobus M Schippers, Damien C Weber, and Antony J Lomax.  Towards FLASH proton therapy: the impact of treatment planning and machine characteristics on achievable dose rates.  Acta Oncol, 58(10):1463–1469, June 2019.
  * [20] Guoliang Zhang, Wenchao Gao, and Hao Peng.  Design of static and dynamic ridge filters for flash–impt: A simulation study.  Medical Physics, 49(8):5387–5399, 2022.