# Beamlet-free optimization for Monte Carlo based treatment planning in proton
therapy

Danah Pross  Université catholique de Louvain, Institut de Recherche
Expérimentale et Clinique (IREC), Center of Molecular Imaging, Radiotherapy
and Oncology, Louvain-La-Neuve, 1348, Belgium  Sophie Wuyckens  Université
catholique de Louvain, Institut de Recherche Expérimentale et Clinique (IREC),
Center of Molecular Imaging, Radiotherapy and Oncology, Louvain-La-Neuve,
1348, Belgium  Sylvain Deffet  Université catholique de Louvain, Institut de
Recherche Expérimentale et Clinique (IREC), Center of Molecular Imaging,
Radiotherapy and Oncology, Louvain-La-Neuve, 1348, Belgium  Edmond Sterpin
Université catholique de Louvain, Institut de Recherche Expérimentale et
Clinique (IREC), Center of Molecular Imaging, Radiotherapy and Oncology,
Louvain-La-Neuve, 1348, Belgium  KU Leuven, Department of Oncology, Laboratory
of Experimental Radiotherapy, Leuven, 3000, Belgium  Particle Therapy
Interuniversity Center Leuven - PARTICLE, Leuven, 3000, Belgium  John A. Lee
Université catholique de Louvain, Institut de Recherche Expérimentale et
Clinique (IREC), Center of Molecular Imaging, Radiotherapy and Oncology,
Louvain-La-Neuve, 1348, Belgium  Kevin Souris  Université catholique de
Louvain, Institut de Recherche Expérimentale et Clinique (IREC), Center of
Molecular Imaging, Radiotherapy and Oncology, Louvain-La-Neuve, 1348, Belgium
Ion Beam Applications SA, Louvain-La-Neuve, 1348, Belgium

###### Abstract

Background: Dose calculation and optimization algorithms in proton therapy
treatment planning often have high computational requirements regarding time
and memory. This can hinder the implementation of efficient workflows in
clinics and prevent the use of new, elaborate treatment techniques aiming to
improve clinical outcomes like robust optimization, arc and adaptive proton
therapy.
Purpose: A new method, namely, the beamlet-free algorithm, is presented to
address the aforementioned issue by combining Monte Carlo dose calculation and
optimization into a single algorithm and omitting the calculation of the time-
consuming and costly dose influence matrix.
Methods: The beamlet-free algorithm simulates the dose in proton batches of
randomly chosen spots and evaluates their relative impact on the objective
function at each iteration. Based on the approximated gradient, the spot
weights are then updated and used to generate a new spot probability
distribution. The beamlet-free method is compared against a conventional,
beamlet-based treatment planning algorithm on a brain case.
Results: The beamlet-free algorithm maintained a comparable plan quality while
reducing the computation time by 70% and the peak memory usage by 95%.
Conclusion: The implementation of a beamlet-free treatment planning algorithm
for proton therapy is feasible and capable of achieving a considerable
reduction of time and memory requirements.

Proton therapy, optimization, treatment planning, beamlet-free, stochastic
gradient descent, micro-optimization

##  I Introduction

In proton therapy, the state-of-the-art is to treat patients with intensity
modulated proton therapyKooy_2015  (IMPT) using pencil beam scanning.
Treatment planning typically involves the precalculation of individual beamlet
dose distributions and subsequent optimization of their respective intensity.
The total dose to voxel $i$ is then given as the linear combination of these
elementary dose distributions

$d_{i}=\sum_{j}D_{ij}x_{j}$ (1)

with $D_{ij}$ as the individual beamlet dose contribution of spot $j$ to voxel
$i$ and $x_{j}$ its corresponding intensity (‘spot weight’). The dose
influence matrix $D$ is the concatenation of the beamlet dose distributions
and usually reaches considerable size, due to the large number of spots in
IMPT.
While dose calculations can be performed using fast analytical dose
calculation algorithms Hong_1996 , their accuracy falls short of the more
elaborate Monte Carlo methods Paganetti_2012 ; Taylor2017 . The drawback of MC
dose engines is their high demand of computational resources that, especially
for large numbers of spots, can prolong treatment planning time considerably.
It is therefore an ongoing effort to reduce computation time for MC methods
and several fast MC dose engines have been introduced, often based on CPU
(Central Processing Unit) or GPU (Graphical Processing Unit) parallelization
Fippel2004 ; Tourovsky_2005 ; Kohno_2011 ; Jabbari2014 ; Souris2016 ;
Schiavi_2017 .
The storage of the beamlet dose distributions in the dose influence matrix,
necessary for optimization, also puts considerable strain on the the computers
memory capacities Ziegenhein_2008 .
Besides the time expended on the repeated calls to the dose calculation
engine, the optimization algorithm also places high demands on the
computational resources and prolongs computation time, due to the very costly
multiplications with the large dose influence matrix involved in the gradient
computation.
These high computational demands, both with respect to computation time and
memory usage, are hindrances in clinical routines that require to be addressed
to enable efficient workflows when a MC dose engine is involved. These
problems become additionally exacerbated in study cases that require a large
number of beamlets, such as robust optimization, adaptive proton therapy or
arc proton therapy. These increasingly elaborate methods seek to improve
treatment quality by reducing dose to organs at risk (OAR) and ensuring target
coverage, but are limited by the current computational demands.
The aim of this work is to present a new method of treatment plan optimization
that substantially reduces the required computation time and memory, while
maintaining the plan quality. The approach combines MC dose calculation and
spot weight optimization into a single efficient algorithm. This algorithm is
referred to as beamlet-free algorithm since it does not require the
computation of individual beamlet dose distributions. The method is applied
here to proton therapy with pencil beam scanning, where the number of beamlets
is typically higher than in photon therapy and where the superior accuracy of
MC methods over analytical dose calculation algorithms is especially
beneficial due to the steep dose gradients involved.

##  II Methods

###  II.A Cost function

Following a classical approach, the objectives penalize deviations of the dose
from a reference value in a region of interest (ROI)Oelfke2001 . For example,
to enforce a prescribed minimum dose $b_{\min}$, the following objective term
may be minimized:

$f_{\min}(\bm{d})=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\min\left(0,d_{i}(\bm{x})-b_{\min}\right)^{2}\enspace,$
(2)

where $N_{r}=|ROI|$ is the number of voxel in the ROI. A similar objective may
be used to enforce a maximum dose $b_{max}$:

$f_{\max}(\bm{d})=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\max\left(0,d_{i}(\bm{x})-b_{\max}\right)^{2}\enspace.$
(3)

Several objectives may be defined. Hence, the general treatment plan
optimization problem takes the following form:

$\displaystyle\min$ $\displaystyle F(\bm{x})=\sum_{k}f_{k}(\bm{d})$
$\displaystyle\mathrm{s.t.}$ $\displaystyle\bm{x}\geq 0\enspace,$

where the index $k$ iterates over the objectives and $f_{k}(\bm{d})$ are
either $f_{\min}(\bm{d})$ or $f_{\max}(\bm{d})$ depending on the objective
type.

###  II.B Beamlet-free algorithm

The idea of beamlet-free optimization is based on continually and
incrementally evaluating the objective function while the dose is being
simulated. Instead of computing each individual beamlet dose distribution, the
dose is continually scored during the MC simulation of the proton batches and
directly added to a total dose distribution. Consequently, spot weights in
their original sense as modulation factors cannot exist anymore. The
equivalent of the original spot weights with respect to the final dose output
would be the number of protons simulated towards each spot. To optimize this
quantity during the ongoing dose simulation, a set of unnormalized
probabilities $p$ is introduced, that are used to generate a spot probability
distribution $P$.
In each iteration of the proposed algorithm, a spot is randomly sampled from
this probability distribution $P$. A batch of protons is simulated towards the
chosen spot and the impact on the objective function is calculated to estimate
the gradient of the objective function with respect to the probability of the
chosen spot.
After all protons of the batch have been simulated, the estimated gradient,
with respect to the sampled spot, is used to update the spot probability,
similarly to Stochastic Gradient Descent (SGD) Robbins1951 .
This yields Algorithm 1, where $p_{j}$ is the probability for a spot $j$, $P$
is the corresponding probability distribution, $h$ is a step size, $c$ is the
batch cost calculated as

$c=\sum_{i}\sum_{k}F(d_{i,k-1})-F(d_{i,k})$ (6)

with $i$ an index iterating over all voxel, $k$ an index over all proton
contributions in a batch, $F(\cdot)$ the objective function as described in
equation II.A and $d_{i,k-1}$ and $d_{i,k}$ as the normalized dose to voxel
$i$ prior to and after the contribution of particle $k$. The particle cost
calculated in 6 corresponds to the finite difference approximation of the
objective function gradient.
The normalization of the dose during the cost computation is necessary to
evaluate the relative impact of simulating particles to a given spot $j$. The
total, unnormalized dose is multiplied with a normalization coefficient
$n(\bm{d})=\frac{\hat{d}}{b}$ where $\hat{d}$ is the mean target dose and $b$
is the prescribed target dose.

Algorithm 1 Beamlet-free algorithm

1:Initialize spot probability vector $\bm{p}$ to unit values

2:Generate probability distribution $P(\bm{p})$

3:Set $n_{\text{simulated}}\leftarrow 0$

4:while $n_{\text{simulated}}\leq n_{\text{total}}$ do

5: Sample spot $j$ from $P$

6: Initiate a batch of $n_{\text{batch}}$ particles for spot $j$

7: Initialize cost $c\leftarrow 0$

8: for $k\leq n_{\text{batch}}$ do

9: Simulate particle $k$

10: Update cost $c\leftarrow c+\sum_{i,k}F(d_{i,k-1})-F(d_{i,k})$

11: $n_{\text{simulated}}\leftarrow n_{\text{simulated}}+1$

12: end for

13: Update spot probability: $p_{j}\leftarrow p_{j}-hc$

14:end while

The quantity to be optimized in this procedure is the spot probability
$\bm{p}$, which is stochastically linked to the number of protons simulated
toward the corresponding spots, the spot weight $\bm{x}$. The objective
function is evaluated over the total dose distribution $\bm{d}$, which depends
on $\bm{x}$, and only indirectly on $\bm{p}$. The finite difference
approximation in 6 concerns the gradient of the objective function with
respect to $\bm{x}$.
There are therefore two separate stochastic processes taking place, the MC
simulation generating the dose distribution and the spot sampling at the
beginning of each iteration which causes the distribution of protons per spot,
$\bm{x}$, and the probability distribution $P$ to converge towards the same
final distribution.
The beamlet-free algorithm was implemented in C around the existing, fast
Monte Carlo code MCsquare Souris2016 .

###  II.C Phantom data

The beamlet-free method is expected to improve both computation time and
memory usage when large numbers of beamlets are involved. To highlight its
benefits and limitations, a complexity analysis was performed. Simulations are
performed on a water phantom containing a target region of various sizes. The
phantom is represented in figure 1. The multiple target sizes result in
geometries with varying numbers of spots. Spot and layer spacing are set to 4
mm and a target margin of 6 mm is applied. The influence of the number of
spots on computation time and memory usage are investigated.
A single beam with couch and gantry angles of $0\degree$ is used for the
analysis. The water phantom has dimensions of 150x150x150 mm3 and the target
region is cubic, with the center located at a depth of 110 mm in beam
direction and side lengths ranging between 20 mm and 60 mm. A representation
of the water phantom can be found in figure 1.
Simulations are repeated 5 times to estimate the impact of statistical
fluctuations.

Figure 1: Water phantom used for the scaling analysis. An examplary target
region with side length of 20 mm is represented in red.

###  II.D Patient data

To demonstrate the high potential of the beamlet-free algorithm, we chose to
benchmark it against a conventional, beamlet-based gradient descent
optimization algorithm. A brain case (base-of-skull tumor) was chosen for the
comparison. Local ethics committee approval was obtained for the use of the
patient data. A treatment plan was generated using a three-beam combination at
a fixed gantry angle G $=90^{\circ}$ and rotating the couch angle C at
$0^{\circ}$, $270^{\circ}$ and $315^{\circ}$. A total of 9812 spots compose
the plan. The ROIs upon which optimization objectives are placed are the
clinical target volume (CTV), optical chiasm (OC), brain stem (BS), right
optical nerve (RON) and left optical nerve (LON). The optimization objectives
are listed in table 1.

Table 1: Optimization objectives ROI | Objective Type | Dose limit | Weight
---|---|---|---
PTV | $D_{\text{min}}$ | 54.0 | 5
PTV | $D_{\text{max}}$ | 54.0 | 5
OC | $D_{\text{max}}$ | 60.0 | 1
BS | $D_{\text{max}}$ | 55.0 | 1
RON | $D_{\text{max}}$ | 55.0 | 1
LON | $D_{\text{max}}$ | 55.0 | 1

The treatment plan optimization is performed using the in-house developed,
research treatment planning system OpenTPSopenTPS2023  which can incorporate
both the beamlet-free method as well as established optimization algorithms
such as second order methods like the BFGS algorithm Broyden1970 ; Shanno1970
; Goldfarb1970 ; Fletcher1970 .

###  II.E Optimization parameters

The total number of particles needed in the Monte Carlo simulation, depend on
the number of spots. For the water phantom, a number of 40 million primary
protons was chosen, which was shown to result in good statistical accuracy for
all cases. For the patient case, which was both larger and more complex, a
higher number of 200 million particles was chosen. All other parameters were
chosen identically for the different simulations.
The beamlet-free algorithm simulates the primary protons in batches of 16
corresponding to $12.5\cdot 10^{6}$ iterations of the first order algorithm
for the patient case. For the comparison, a beamlet matrix with $5\cdot
10^{4}$ protons per beamlet is simulated and subsequently optimized using the
scipy scipy  FORTRAN implementation (wrapped in Python) of the second-order
optimization algorithm L-BFGS Liu1989  (Limited-memory Broyden-Fletcher-
Goldfarb-Shanno algorithm). The optimizer continues until it meets the
tolerance criterion which is defined as the absolute difference in ojective
function values between iterations and set to $10^{-4}$. A final dose
computation with 40 million protons is performed to match the statistical
accuracy of the beamlet-free method. Both the beamlet calculation and the
final dose computation are performed using the MCsquare dose engine that is
also utilized in the beamlet-free algorithm. The number of primaries for the
MC dose calculation in the beamlet-free algorithm was intentionally picked so
that the optimized plan quality would reach the one obtained with the beamlet-
based method and consequently ease the comparison.
Simulations were performed on a computation server with a total of 80 CPU
cores and 526 GB RAM.

##  III Results

###  III.A Scaling analysis

Computation time and peak memory usage are plotted against the number of spots
for both methods and compared in Figure 2. The values are calculated as the
mean over 5 runs and the errorbar represents the standard deviation. The
difference in target coverage was less than 0.20 Gy in all cases.

Figure 2: Memory usage (a) and computation time (b) required by the beamlet-
based and beamlet-free methods versus the number of spots for a water phantom
with different target sizes.

###  III.B Patient data

Dose distributions obtained by both methods are depicted in figure 3. In both
cases all OAR objectives could be met and the target D95/D5 is within $5\%$ of
the prescription, as it was reported in table 2.
For the OAR, $D_{\text{max}}$ objectives are considered as met, if the $D_{5}$
value of the ROI is below the dose limit. Similarly, $D_{\text{mean}}$
objectives are considered as met, if the $D_{\text{mean}}$ value of the ROI is
below the dose limit. The target objectives are considered as met if the
$D_{95}$/$D_{5}$ value is within $\pm 5\%$ of the prescription dose.

Figure 3: Comparison of the optimized dose distribution obtained by the beamlet-based (a) and beamlet-free method (b) and the computed dose volume histogram (c) for beamlet-based (dashed) and beamlet-free (solid). Table 2: Clinical goals evaluation for the beamlet-based and beamlet-free method. ROI | Criterium |  | Dose
---
limit
(Gy)
|  $D_{\text{Criterium}}$ (Gy)
---
Beamlet-
based
|  $D_{\text{Criterium}}$ (Gy)
---
Beamlet-
free
PTV | $D_{95}$ | 51.3 | 52.73 | 52.28
PTV | $D_{5}$ | 56.7 | 55.24 | 55.36
PTV | $D_{98}$ | 52.98 | 52.35 | 51.53
PTV | $D_{2}$ | 55.08 | 55.58 | 55.87
OC | $D_{5}$ | 60.0 | 55.08 | 55.09
OC | $D_{2}$ | 60.0 | 55.43 | 55.36
OC | $D_{mean}$ | 50.0 | 53.85 | 53.82
BS | $D_{5}$ | 55.0 | 53.11 | 52.55
BS | $D_{2}$ | 55.0 | 54.03 | 53.82
BS | $D_{mean}$ | 8.50 | 24.79 | 24.16
RON | $D_{5}$ | 55.0 | 51.74 | 49.99
RON | $D_{2}$ | 55.0 | 53.10 | 51.79
LON | $D_{5}$ | 55.0 | 54.57 | 54.58
LON | $D_{2}$ | 55.0 | 54.96 | 55.09

Due to the stochastic nature of both methods, small fluctuations in the final
dose distributions are to be expected. The particle number and the total
number of iterations have to be chosen sufficiently high to ensure an adequate
level of confidence. To this end, both methods were run 5 times. The
fluctuations on the target coverage,$D_{5}-D_{95}$, were less than $4\%$ for
both methods.
The objective function convergence over the course of the optimization is
shown in figure 4. The objective function values are plotted against the
number of iteration of the respective optimization algorithm. The relative
behaviour of the objective function can be compared, and a similar behaviour
of rapidly decreasing objective function values can be observed.

Figure 4: (a) Convergence plot of the beamlet-based algorithm (scipy L-BFGS)
and (b) convergence plot of the beamlet-free algorithm.

The peak memory required during the dose calculation and optimization, could
be reduced from $29$ GB for the beamlet-based method to about $1.561$ GB for
the beamlet-free method. It is worth noticing that the peak memory for the
beamlet-based method was recorded during the import of the individual beamlet
dose distributions, compressed as sparce matrices, into the dose influence
matrix (stacked as columns). Other steps in the treatment planning process,
however, required higher peak memory usage than the beamlet-free method needed
during dose calculation and optimization. This peak memory was recorded as
about $3.569$ GB. It depends on the used treatment planning system and might
be improved independently from the presented algorithm. A memory comparison is
drawn in Table 3. Besides the peak RAM usage during simulation, it also lists
the required storage disk space of the files needed for the dose calculation
and optimization. The difference between the two methods with respect to
utilized disk space can be explained through the storage of the large dose
influence matrix that requires 4.4 GB alone.

Table 3: Comparison of memory usage for the beamlet-based and beamlet-free method. | Beamlet-based | Beamlet-free
---|---|---
Peak RAM usage [GB] | $29$ | $1.561$
Storage [GB] | $6.4$ | $2.2$

The total computation times averaged over 5 runs are reported in table 4.

Table 4: Comparison of computation time for the beamlet-based and beamlet-free method averaged over 5 runs. Process |  $t$ [s]
---|---
BL-based | Plan Creation | $299\pm 2$
Beamlet Comp. | $3356\pm 24$
Optimization | $12601\pm 1714$
Final Dose Comp. | $301\pm 1$
Total | $\mathbf{16558\pm 1714}$
BL-free | Plan Creation | $299\pm 2$
Optimization | $4697\pm 129$
Total | $\mathbf{4996\pm 129}$

##  IV Discussion

A novel treatment planning method, not limited to proton therapy, has been
introduced. It merges dose calculation and optimization into a single process,
resulting in improved computational efficiency and accuracy. The so-called
beamlet-free algorithm was able to achieve comparable plan quality to a
conventional beamlet-based algorithm.
The comparability of the two methods is limited by the fact that different
parameters have to be chosen for each of them. These parameters have to be
chosen in a way that allows a balance between keeping computation time low
enough to enable clinical workflows while achieving a high statistical
confidence. The statistical uncertainty inevitably present in MC-based
treatment planning stems from the fact that MC methods, and in the case of the
beamlet-free algorithm also the optimization, are based on stochastic
processes. The two methods have different parameters that may be tuned to this
end. For this comparison, the same number of primary protons was therefore
used in the final dose calculation to ensure a comparable level of uncertainty
for the MC simulation. The iteration number of the optimization algorithms is
not comparable, since the direction of search is only partially estimated at
each step of the beamlet-free algorithm which consequently needs a lot more
iterations to converge. In both cases, the number of iterations utilized was
considered as realistic to produce adequate plan quality in a reasonable time-
frame. Another parameter that impacts computation time and quality of the
results is the number of protons per beamlet in the beamlet-based method and
the batch size in the beamlet-free. While they are hardly equivalent to one
another, they play a similar role during the optimization. They are used as
input for the optimization algorithm. Both need to be sufficiently high to
prevent optimization based on inaccurate dose estimations, but still
considerate of the total computation time. In beamlet-free treatment planning
the three following parameters, total number of primary protons, number of
iterations and batch size, depend on each other. However, they can be chosen
independently in the beamlet-based case. This places some restrictions on the
user, but might also be beneficial, since it reduces the risk of errors.
Both methods are able to fulfill all clinical objectives. Target coverage
worsened slightly for the beamlet-free case, while OAR sparing improved for
most organs. Altogether, the difference between the two methods is small, but
statistically significant. It might stem from a number of sources. For
example, a possible explanation could be misguided optimization steps based on
insufficiently approximated gradients.
Regarding computational performances, the peak memory required during
treatment planning could be reduced by $95\%$ through the use of the beamlet-
free algorithm. The total computation time could be reduced by $70\%$ through
use of the beamlet-free algorithm. Notice that the two methods profit
unequally from the number of available computation threads. The MC code,
MCsquare, utilizes CPU parallelization to improve computation time and is
slowed down considerably, if the number of CPU threads is reduced. The
beamlet-free algorithm also utilizes MCsquare for the dose simulation, but the
optimization slows down the MC simulation itself. The beamlet-free algorithm
still profits from having multiple CPU threads at its disposal, but especially
very high numbers of threads cannot achieve comparable simulation times to the
pure MC simulation. The update procedure of the mean target dose is here a
bottleneck in the parallelization, that prevents the beamlet-free algorithm to
draw the full benefit from very large numbers of CPU threads, probably due to
conflicting accesses to the same chunks of memory. It should be noted that the
presented data was obtained from simulations performed on a computation server
with 80 CPU cores and that the benefit of the beamlet-free algorithm is
considerably higher, for computations performed with a limited number of CPUs.
For a more thorough comparison, both optimizers should be implemented within
the same framework, to ensure identical conditions. Small differences can be
expected by using an existing external code for the beamlet-based part (i.e.
Scipy optimizer in our case) of the comparison.
The complexity analysis revealed that, in beamlet-based treatment planning,
computation time and peak memory usage increase almost linearly with the
number of spots. On the other hand, the beamlet-free method can achieve
comparable plan quality without being affected by the number of spots,
resulting in independent computation time and memory usage. The peak memory
usage from the beamlet-free method was recorded below that of the beamlet-
based simulation in all cases. With respect to computation time, for low
numbers of spots the beamlet-based method outperformed the beamlet-free
method, but for increasing numbers of spots, the beamlet-based method gets
slowed down considerably. It is possible to determine roughly from which
number of spots it would be beneficial to utilize the beamlet-free instead of
a beamlet-based method. It should be kept in mind, however, that factors like
the machine computations are performed on, the exact beamlet-based optimizer
and implementation chosen for the comparison and the depth location of the
target and number of simulated particles, can impact the two methods
differently and render it impossible to fix any specific number.
It can be observed that for very low numbers of spots, the beamlet-free method
can even bring disadvantages, but already for numbers of spots regularly
encountered in clinical routines, the beamlet-free method starts to outperform
conventional beamlet-based treatment planning. The benefits of the beamlet-
free method are particularly highlighted in methods that require a large
number of beamlets, such as robust optimization. In the robustness approach,
traditional strategies make use of a set of worst-case scenarios
Unkelbach_2007 ; Pflugfelder_2008 , compiling the treatments uncertainties. If
a beamlet-based approach is used, the very costly influence dose matrix must
be calculated for each scenario and this further brings a computational burden
during the optimization. Additionally, a more recent development sees the
emergence of arc proton therapy treatment modality DING20161107 ; gu2020novel
; Engwall_2022 ; WUYCKENS2022105609  that delivers the dose across a very
large set of irradiation angles, shaping the arc. However, increasing the
number of beams translates into increasing the number of spots, resulting in
producing a larger beamlet matrix than the “classical” one obtained with an
IMPT treatment. Finally, the beamlet-free method could also play an
interesting role for adaptive proton therapy where a new treatment plan might
be needed in good time for an anatomy that has changed.

##  V Conclusion

The proposed beamlet-free algorithm could successfully generate treatment
plans of comparable quality to conventional treatment planning methods. Memory
requirements could be reduced by $95\%$ and the computation time by $70\%$ on
an illustrative patient case. The scaling analysis revealed that the benefit
to be drawn from the beamlet-free method increases for higher numbers of
spots, which makes the method promising for new emerging treatment techniques
such as robust optimization, adaptive and arc proton therapy.

##  VI Aknowledgments

D.P. is funded through a Télévie grant by the F.R.S.-FNRS (Fonds de la
Recherche scientifique). S.W. is funded by the Walloon Region as part of the
Arc Proton Therapy convention (Pôles Mecatech et Biowin). J.A.L. is a Research
Director with the F.R.S.-FNRS (Fonds de la Recherche scientifique).

##  VII Conflict of Interest

Kevin Souris is an employee at Ion Beam Applications SA.

## References

  * (1) H. M. Kooy and C. Grassberger,  Intensity modulated proton therapy,  The British Journal of Radiology 88, 20150195 (2015),  PMID: 26084352.
  * (2) L. Hong, M. Goitein, M. Bucciolini, R. Comiskey, B. Gottschalk, S. Rosenthal, C. Serago, and M. Urie,  A pencil beam algorithm for proton dose calculations,  Physics in Medicine & Biology 41, 1305 (1996).
  * (3) H. Paganetti,  Range uncertainties in proton therapy and the role of Monte Carlo simulations,  Physics in Medicine & Biology 57, R99 (2012).
  * (4) K. S. F. . F. D. S. Taylor, P. A.,  Pencil Beam Algorithms Are Unsuitable for Proton Dose Calculations in Lung,  International journal of radiation oncology, biology, physics 99, 750–756 (2017).
  * (5) M. Fippel and M. Soukup,  A Monte Carlo dose calculation algorithm for proton therapy,  Medical Physics 31, 2263–2273 (2004).
  * (6) A. Tourovsky, A. J. Lomax, U. Schneider, and E. Pedroni,  Monte Carlo dose calculations for spot scanned proton therapy,  Physics in Medicine & Biology 50, 971 (2005).
  * (7) R. Kohno, K. Hotta, S. Nishioka, K. Matsubara, R. Tansho, and T. Suzuki,  Clinical implementation of a GPU-based simplified Monte Carlo method for a treatment planning system of proton beam therapy,  Physics in Medicine & Biology 56, N287 (2011).
  * (8) K. Jabbari and J. Seuntjens,  A fast Monte Carlo code for proton transport in radiation therapy based on MCNPX,  Journal of Medical Physics 39, 156–163 (2014).
  * (9) K. Souris, J. A. Lee, and E. Sterpin,  Fast multipurpose Monte Carlo simulation for proton therapy using multi- and many-core CPU architectures,  Medical Physics 43, 1700–1712 (2016).
  * (10) A. Schiavi, M. Senzacqua, S. Pioli, A. Mairani, G. Magro, S. Molinelli, M. Ciocca, G. Battistoni, and V. Patera,  Fred: a GPU-accelerated fast-Monte Carlo code for rapid treatment plan recalculation in ion beam therapy,  Physics in Medicine & Biology 62, 7482 (2017).
  * (11) P. Ziegenhein, J. J. Wilkens, S. Nill, T. Ludwig, and U. Oelfke,  Speed optimized influence matrix processing in inverse treatment planning tools,  Physics in Medicine & Biology 53, N157 (2008).
  * (12) U. Oelfke and T. Bortfeld,  Inverse planning for photon and proton beams,  Medical Dosimetry 26, 113–124 (2001).
  * (13) H. E. Robbins,  A Stochastic Approximation Method,  Annals of Mathematical Statistics 22, 400–407 (1951).
  * (14) S. Wuyckens, D. Dasnoy, G. Janssens, V. Hamaide, M. Huet, E. Loÿen, G. R. de Hertaing, B. Macq, E. Sterpin, J. A. Lee, K. Souris, and S. Deffet,  OpenTPS – Open-source treatment planning system for research in proton therapy, 2023.
  * (15) C. G. BROYDEN,  The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations,  IMA Journal of Applied Mathematics 6, 76–90 (1970).
  * (16) D. Shanno,  Conditioning of quasi-Newton methods for function minimization,  Mathematics of Computation. 24, 647–656 (1970).
  * (17) D. Goldfarb,  A family of variable-metric methods derived by variational means,  Mathematics of Computation. 24, 23–26 (1970).
  * (18) R. Fletcher,  A new approach to variable metric algorithms,  The Computer Journal 13, 317–322 (1970).
  * (19) P. Virtanen et al.,  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,  Nature Methods 17, 261–272 (2020).
  * (20) D. Liu and J. Nocedal,  On the limited memory BFGS method for large scale optimization,  Mathematical Programming 45, 503–528 (1989).
  * (21) J. Unkelbach, T. C. Y. Chan, and T. Bortfeld,  Accounting for range uncertainties in the optimization of intensity modulated proton therapy,  Physics in Medicine & Biology 52, 2755 (2007).
  * (22) D. Pflugfelder, J. J. Wilkens, and U. Oelfke,  Worst case optimization: a method to account for uncertainties in the optimization of intensity modulated proton therapy,  Physics in Medicine & Biology 53, 1689 (2008).
  * (23) X. Ding, X. Li, J. M. Zhang, P. Kabolizadeh, C. Stevens, and D. Yan,  Spot-Scanning Proton Arc (SPArc) Therapy: The First Robust and Delivery-Efficient Spot-Scanning Proton Arc Therapy,  International Journal of Radiation Oncology*Biology*Physics 96, 1107–1116 (2016).
  * (24) W. Gu, D. Ruan, Q. Lyu, W. Zou, L. Dong, and K. Sheng,  A novel energy layer optimization framework for spot-scanning proton arc therapy,  Medical physics 47, 2072–2084 (2020).
  * (25) E. Engwall, C. Battinelli, V. Wase, O. Marthin, L. Glimelius, R. Bokrantz, B. Andersson, and A. Fredriksson,  Fast robust optimization of proton PBS arc therapy plans using early energy layer selection and spot assignment,  Physics in Medicine & Biology 67, 065010 (2022).
  * (26) S. Wuyckens, M. Saint-Guillain, G. Janssens, L. Zhao, X. Li, X. Ding, E. Sterpin, J. A. Lee, and K. Souris,  Treatment planning in arc proton therapy: Comparison of several optimization problem statements and their corresponding solvers,  Computers in Biology and Medicine , 105609 (2022).


#  $f\left(R,\square R\right)$-gravity and equivalency with the modified GUP
Scalar field models

Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of Systems Science,
Durban University of Technology, Durban 4000, South Africa  Departamento de
Matemáticas, Universidad Católica del Norte, Avda. Angamos 0610, Casilla 1280
Antofagasta, Chile

###### Abstract

Inspired by the generalization of scalar field gravitational models with a
minimum length we study the equivalent theory in modified theories of gravity.
The quadratic Generalized Uncertainty Principle (GUP) gives rise to a deformed
Heisenberg algebra in the application, resulting in the emergence of
additional degrees of freedom described by higher-order derivatives. The new
degrees of freedom can be attributed to the introduction of a new scalar
field, trasforming the resulting theory into a two-scalar field theory. Thus,
in order to describe all the degrees of freedom we investigate special forms
of the sixth-order modify $f\left(R,\square R\right)-$theory of gravity, where
the gravitational Lagrangian has similar properties to that of the GUP scalar
field theory. Finally, the cosmological applications are discussed, and we
show that the de Sitter universe can be recovered without introducing a
cosmological constant.

Cosmology; $f\left(R,\square R\right)$-gravity; Generalized Uncertainty
Principle

##  I Introduction

Dark energy is an exotic matter source introduced in Einstein’s field
equations of General Relativity in order to explain the cosmic acceleration as
it is observed by the cosmological data rr1 ; Teg ; Kowal ; Komatsu ; suzuki11
. Currently, there are no observable phenomena directly linked to the nature
and characteristics of dark energy. As a result, the physical nature and the
origin of dark energy it is up for debate in the scientific community.
Introducing the cosmological constant in the Einstein-Hilbert Action Integral
represents one of the simplest approaches to address the issue of dark energy.
The $\Lambda$CDM cosmology is an analytic solution of the field equations of
General Relativity for a spatially flat Friedmann–Lemaître–Robertson–Walker
(FLRW) geometry with a nonzero cosmological constant term and a pressureless
dust fluid. $\Lambda$CDM has achieved considerable success in describing a
large span of astronomical and cosmological data. Nevertheless, the
cosmological constant can not explain the complete evolution of the cosmic
history l01 ; l02 , we refer the reader to the recent discussion l03 .

In order to overcome the problems of the cosmological constant and to explain
the late-time cosmic acceleration, in recent years, cosmologists have proposed
various models which can categorized to two large families of theories, the
dark energy theories and the modified theories of gravity. For the dark energy
models, an energy-momentum tensor that attributes the new degrees of freedom
is introduced in field equations of General Relativity. The dynamics driven by
the newly introduced degrees of freedom provide an explanation for various
cosmological phenomena. Some of the most common dark energy models are the
quintessence scalar field ratra ; q5 ; qq1 ; qq4 , phantom scalar field q22 ;
q23 , chameleon mechanism qq01 , scalar-tensor models st1 ; st2 ; st3 ,
Galileons st4 ; st5 ; multi-scalar field models st6 ; st7 ; st8 , Chaplygin
gas-like fluids q20 ; q21 ; q21a ; q21b ; q21c , k-essence q26 ; q26a ,
tachyons q25 ; q25a ; q25b . On the other hand, in modified theories of
gravity the Einstein-Hilbert Action Integral is modified with the introduction
of geometric invariants md1 ; md2 ; md3 . As a result, the gravitational field
equations are modified such that new geometrodynamical components to be
introduced and provide an effective geometric matter source to explain the
acceleration of the universe cp1 . The ”zoology” of modified theories of
gravity can be categorized based on the geometric invariant used to modify the
gravitational Action Integral and the order of derivatives involved. Within
the realm of modified theories of gravity, a particular family of interest
comprises the so-called $f\left(X\right)~{}$-theories, where the gravitational
Action Integral is a function $f$ of the geometric invariant $X$. The latter
can be the Ricci scalar leads to $f\left(R\right)~{}$-gravity Buda , the
torsion scalar $T$ of the Weitzenböck connection in teleparallelism Ferraro ,
the non-metricity scalar $Q$ in symmetric teleparallel theory f6 , the Gauss-
Bonnet term bb1 , for other proposed theories see for instance bb2 ; bb3 ; bb4
; bb5 ; bb6 ; bb7  and references therein.

$f\left(R\right)~{}$-gravity Sotiriou ; odin1  has been widely studied in the
literature with many interesting results in cosmological studies cf1 ; cf2 ;
cf3  as also important results in the description of astrophysical objects cf4
; cf5 ; cf6 . The quadratic theory of gravity $f\left(R\right)~{}=R+qR^{2}$
qua1 ; qua2 ; qua3 , has been used successfully planck2015  for the
description of another inflationary epoch of our universe. Specifically, the
quadratic $f\left(R\right)$-gravity has been used as a mechanism for
acceleration in the early stages of the universe star . The quadratic term
$R^{2}$ follows from the analytic expressions for the quantum-gravitational
effects in the one-loop approximation. Indeed the origin of the $R^{2}$ term
is the vacuum polarization of the physical space df1 ; df2 .

The existence of a minimum length, i.e. maximum energy in nature, is supported
by various approaches to quantum gravity, like string theory, doubly special
relativity and the black hole physics. The incorporation of a minimum length
necessitates a modification of Heisenberg’s Uncertainty Principle, giving rise
to the Generalized Uncertainty Principle (GUP) Maggiore . By modifying the
Uncertainty Principle governing quantum observables, we arrive at a revised
definition of the Heisenberg algebra, which in turn leads to adjustments in
the Poisson brackets in the classical limits. cas1 ; cas2 ; cas3 . The
modified Poisson brackets revise the equations of motion such that new degrees
of freedom to be introduced. There is a plethora of studies of GUP in
gravitational physics, see for instance mc1 ; mc2 ; mc3 . As far as the
cosmological constant is concerned, it has been found that it is related to
the GUP and the minimum length cc1 ; cc2 . Another extension of the
Uncertainty Principle in the Extended Uncertainty Principle (EUP) ep1 ; ep2 ,
however in this work we are focus in GUP. We refer the reader to the recent
review gup11

The quintessence model, employed to describe dark energy in angup1
incorporates modifications to the scalar field Lagrangian through the
application of the quadratic GUP, that is, in the equation of motion for the
scalar field, i.e. the Klein-Gordon equation, new higher-order derivatives
have been introduced, as a result of the deformed Heisenberg algebra. These
new terms are related to the existence of the minimum length. The higher-order
derivatives can be described by a second-scalar field and this modified
quintessence model is equivalent to a multi-scalar field cosmological model
with interaction between the two scalar fields. It was found that the GUP
components affect the cosmological evolution of a FLRW geometry not only in
the early stages of the universe but also in the late-universe, while the
effects for the existence of the deformed Heisenberg algebra are observable in
the cosmological perturbations angup2 . The case where a matter source is
included in the field equations coupled to the scalar field was the subject of
study in angup3 . Furthermore, in angup4  the deformed Heisenberg algebra for
the quadratic GUP was considered to study the effects of the minimum length in
the case of a gravitational theory, which satisfies Mach’s Principle.
Specifically, the case of the Bran-Dicke scalar field was considered modified
by the quadratic GUP. It was found that the modified theory is a multi-scalar
field theory with different dynamics and evolution from the unmodified theory.
An important result is that the nature of the asymptotic solutions does not
depend on the value of the Brans-Dicke parameter, which is different from the
case of the unmodified theory. Last but not least, the effects of the GUP
modification in the Brans-Dicke theory are observable not only in the early
stages of the universe but also in the late-time.

There exists a unique connection between some modified theories of gravity and
scalar field models. Indeed, the new degrees of freedom provided by the
geometric scalars in modified theories of gravity can be attributed to scalar
fields. For instance $f\left(R\right)~{}$-gravity Sotiriou  is equivalent with
the so-called O’Hanlon theory ohan  which belongs to the family of scalar-
tensor theories far1  and specifically to the Brans-Dicke gravity brans .
Recall that $f\left(R\right)$-gravity is fourth-order theory for a nonlinear
function $f$, while it reduces to a second-order theory in the limit of
General Relativity when $f$ is a linear function. Thus, with the introduction
of a Lagrange multiplier the field equations can be written in the equivalent
order of second-order derivatives but in the same time the number of the
dependent variables increases. In general, Lagrange multipliers are applied
for the introduction of constraints for the cosmological model. There are
various studies in the literature for the application of Lagrange multipliers
in gravitational theories, see for instance lan1 ; lan2 ; lan3 .

In this study, we investigate if there exists a modified theory of gravity
where we can recover the quadratic GUP corrections that follow from the
deformation of the Heisenberg algebra in scalar field theories. With such
analysis, we will be able to find the geometric equivalent of GUP in modified
theories of gravity. The GUP scalar field models studied before, in the
Einstein and Jordan frames angup1 ; angup4 , are second-order multi-scalar
field models, specifically the matter source is attributed to two-scalar
fields. Hence, in order to be able to recover a two scalar field theory in
modified theories of gravity we shall consider a sixth-order theory.
$F\left(R,\square R\right)~{}$-theory has been introduced before sx1 ; sx2  as
a sixth-order gravity and extension of $f\left(R\right)~{}$ theory, where
$\square$ is the Laplace operator. The introduction of higher-order
derivatives in the gravitational Action Integral is in agreement with quantum
gravity sx1 . The effects of the $\square R$ terms in $F\left(R,\square
R\right)~{}-$theory has been widely studied before in the description of
inflation sx3 ; sx4 ; sx4a . A detailed analysis of the cosmological dynamics
in $F\left(R,\square R\right)$ performed recently in sx5  where it was found
that higher-order terms can dominate the evolution of the universe. For more
applications of $F\left(R,\square R\right)~{}$-theory in gravitational physics
we refer the reader to sx6 ; sx7 ; sx8 ; sx9 ; sx10  and references therein.
The structure of the paper is as follows.

In Section II we present the basic properties and definitions of GUP and we
focus on the case of quadratic GUP. We define the deformed Heisenberg algebra,
and we derive the modified Klein-Gordon equation for a spin-0 particle. The
latter modified Lagrangian is used in Section III in order to introduce the
effects of the deformed Heisenberg algebra in scalar field cosmological
models. The modified quintessence and modified Brans-Dicke models are
presented. In Section IV we consider the $F\left(R,\square R\right)$-gravity
which is a theory of gravity of sixth-order. We introduce Lagrange multipliers
in order to increase the number of the dependent variables with the
introduction of scalar fields, and at the same time reduce the theory into a
second-order gravitational model. We found that a separable function
$F\left(R,\square R\right)~{}$ is equivalent to two-scalar field theory with
similar properties to that of GUP scalar field models. In Section V we focus
in the case of $F\left(R,\square R\right)~{}=R+K\left(~{}\square R\right)~{}$
gravity, where the term $K\left(~{}\square R\right)~{}$ introduces similar
corrections terms in the field equations as that of the minimum length of GUP.
Indeed, the $K\left(~{}\square R\right)~{}$ we can say that follows from the
deformation algebra of the Einstein-Hilbert Action Integral. For a spatially
flat FLRW background geometry in Section V it was found that the de Sitter
universe is a unique attractor for the cosmological solution without necessary
introduce a cosmological constant term. Thus, the existence of a minim length
in the early universe leads to an accelerated universe. Finally, in Section VI
we draw our conclusions.

##  II Generalized Uncertainty Principle

The existence of a minimum length leads to the modification of the
Heisenberg’s uncertainty principle as

$\Delta X_{i}\Delta P_{j}\geqslant\frac{\hbar}{2}[\delta_{ij}(1+\beta
P^{2})+2\beta P_{i}P_{j}],$ (1)

where parameter $\beta$ is the deformed parameter defined as
$\beta={\beta_{0}}/{M_{Pl}^{2}c^{2},}$ where $M_{Pl}$ is the Planck mass and
$M_{Pl}c^{2}$ is the Planck energy, or equivalently
$\beta={\beta_{0}\ell_{Pl}^{2}}/{\hbar^{2}}$ where $\ell_{Pl}$ ($\approx
10^{-35}~{}m)$ is the Planck length. Parameter $\beta_{0}$ usually is selected
to be positive and equal to one, however in order to have observable quantum
effects, the parameter $\beta_{0}$ can have different values Vagenas .

The usual choice of the parameter $\beta_{0}$ is $\beta_{0}=1$, however, such
a choice could lead do not observable quantum effects; however in Vagenas , it
has been shown that the dimensionless parameter $\beta_{0}$ could has upper
bound such as $\beta_{0}>>1$. However, there are studies where shown that the
deformation parameter $\beta_{0}$ can be negative neg1 ; neg2 ; neg5 ; neg6 .

Therefore, the modified uncertainty principle (1) leads to the deformed
Heisenberg algebra Kemph1 ; Kemph2

$[X_{i},P_{j}]=i\hbar[\delta_{\alpha\beta}(1+\beta_{0}\frac{\ell_{Pl}^{2}}{2\hbar^{2}}P^{2})+\beta_{0}\frac{\ell_{Pl}^{2}}{\hbar^{2}}P_{\alpha}P_{\beta}].$
(2)

Hence, the coordinate representation of the momentum operator, which satisfies
the commutation relation (2), can be defined as $P_{i}=p_{i}(1+\beta p^{2}),$
where we have selected to keep underformed the position underformed, that is,
$X_{i}=x_{i}$. Representation $(x,p)$ is the canonical representation
satisfying $[x_{i},p_{j}]=i\hbar\delta_{ij}$.

In the relativistic limit the commutation relation (2) reads Moayedi

$[X_{\mu},P_{\nu}]=-i\hbar[(1+\beta_{0}\frac{\ell_{Pl}^{2}}{2\hbar^{2}}(\eta^{\mu\nu}P_{\mu}P_{\nu}))\eta_{\mu\nu}+\beta_{0}\frac{\ell_{Pl}^{2}}{\hbar^{2}}P_{\mu}P_{\nu}],$
(3)

where $\eta_{\mu\nu}$ is the flat metric. Therefore, the corresponding
deformed operators for (3) are
$P_{\mu}=p_{\mu}(1-\beta(\eta^{\alpha\gamma}p_{\alpha}p_{\gamma}))~{},~{}~{}X_{\nu}=x_{\nu}.$

In the relativistic limit the equation of motion for a spin-0 particle with
rest mass zero is

$\left[\eta^{\mu\nu}P_{\mu}P_{\nu}-\left(mc\right)^{2}\right]\Psi=0,$ (4)

that is,

$\square\Psi-2\beta\hbar^{2}\square\left(\square\Psi\right)+\left(\frac{mc}{\hbar}\right)^{2}\Psi+O\left(\beta^{2}\right)=0.$
(5)

where the Laplace operator for the metric $\eta_{\mu\nu}$ is marked with the
symbol $\square$ . For a generic metric tensor $g_{\mu\nu}$ the Laplace
operator is defined as $\Delta=\square$, where
$\Delta=\frac{1}{\sqrt{-g}}\partial_{\mu}\left(g^{\mu\nu}\sqrt{-g}\partial_{\mu}\right)$.
In equation (5) term $2\beta\hbar^{2}\square\left(\square\Psi\right)$ is the
quantum correction term which follows from the existence of the minimum
length. Because the quantum correction terms introduce the fourth-order
derivatives in the Klein-Gordon equation, the $\beta^{2}\rightarrow 0$,
equation (5) is a singular perturbative system which means that there is a
solution where the term $2\beta\hbar^{2}\square\left(\square\Psi\right)$
dominates and drives the dynamics. For the discussion of inner and outer
solutions of singular pertubative differential equations we refer the reader
in Tikhonov .

The modified Klein-Gordon equation (5) is a fourth-order partial differential
equation. It can be written in the equivalent form of two second-order
differential equations by introducing the new scalar field $\Phi=\square\Psi$.
Hence, equation (5) in the limit $\beta^{2}\rightarrow 0$ becomes,

$\square\Psi-\Phi=0,$ (6)
$2\beta\hbar^{2}\square\Phi+\left(mc\right)^{2}\Psi+\Phi=0.$ (7)

Let us now derive the Lagrangian for the two scalar field model with equations
of motion (6) and (7).

We consider the Action Integral of the equation (5)

$S_{KG}^{\operatorname{mod}}=\int
dx^{4}\left(\frac{1}{2}\eta^{\mu\nu}\mathcal{D}_{\mu}\Psi\mathcal{D}_{\nu}\Psi-\frac{1}{2}\left(\frac{mc}{\hbar}\right)^{2}\Psi^{2}\right),$
(8)

where now $\mathcal{D}_{\mu}$ is the deformed operator defined as

$\mathcal{D}_{\mu}=\nabla_{\mu}+\beta\hbar^{2}\nabla_{\mu}\square.$ (9)

Hence, expression (8) becomes

$S_{KG}^{\operatorname{mod}}=\int
dx^{4}\left(\frac{1}{2}\eta^{\mu\nu}\left(\nabla_{\mu}\Psi\nabla_{\nu}\Psi+2\beta\hbar^{2}\left(\nabla_{\nu}\nabla_{\mu}\square\right)\Psi\right)-\frac{1}{2}\left(\frac{mc}{\hbar}\Psi\right)^{2}\right),$
(10)

We include the Lagrange multiplier $\lambda$ into the dynamical system and the
second scalar field $\Phi$ with constraint $\Phi=\square\Psi$, therefore the
latter Action Integral after integration by parts is written as follows

$S_{KG}^{\operatorname{mod}}=\int
dx^{4}\sqrt{-g}\left(\frac{1}{2}\eta^{\mu\nu}\nabla_{\mu}\Psi\nabla_{\nu}\Psi+2\beta\hbar^{2}\eta^{\mu\nu}\nabla_{\mu}\Psi\nabla_{\nu}\Phi+\beta\hbar^{2}\Phi^{2}-\frac{1}{2}\left(\frac{mc}{\hbar}\right)^{2}\Psi^{2}\right).$
(11)

##  III Scalar field theories modified by GUP

We review previous results on the application of the quadratic GUP in scalar
field theories.

Consider the Scalar-tensor Action Integral far1

$S=\int
dx^{4}\sqrt{-g}\left[\frac{1}{2}F\left(\phi\right)^{2}R-L_{\phi}\right],$ (12)

where $F\left(\phi\right)$ is the coupling function and $L_{\phi}$ is the
Lagrangian for the scalar field, that is

$L_{\phi}\left(x^{\mu},\phi,\nabla_{\mu}\phi\right)=\frac{\omega}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi+V\left(\phi\right)$
(13)

where $\omega$ is a parameter which defines the nature of the scalar field.
Lagrangian function (13) is that of the underformed Heisenberg algebra. Hence,
after the application of the quadratic GUP, we end with Lagrangian function
(11) which reads

$L_{\phi}^{GUP}\left(x^{\mu},\phi,\psi,\nabla_{\mu}\phi,\nabla_{\mu}\psi\right)=\frac{\omega}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi+V\left(\phi\right)+\beta\left(2g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\psi+\psi^{2}\right)$
(14)

where $\psi$ is the second scalar field which follows from the higher-order
derivatives.

We replace in (14) , (12) and we end with the Action Integral

$S=\int
dx^{4}\sqrt{-g}\left[\frac{1}{2}F\left(\phi\right)R-\frac{\omega}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-V\left(\phi\right)-\beta\left(2g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\psi+\psi^{2}\right)\right].$
(15)

###  III.1 Quintessence and phantom cosmologies

In the case of a spatially flat FLRW geometry

$ds^{2}=-dt^{2}+a^{2}\left(t\right)\left(dx^{2}+dy^{2}+dz^{2}\right),$ (16)

and for a minimally coupled scalar field, that is $F\left(\phi\right)=1$, the
cosmological field equations are

$3H^{2}-\frac{1}{2}\dot{\phi}^{2}-2\beta\dot{\phi}\dot{\psi}-V\left(\phi\right)+\beta\psi^{2}=0,$
(17)
$2\dot{H}+3H^{2}+\frac{1}{2}\dot{\phi}^{2}+2\beta\dot{\phi}\dot{\psi}-\left(V\left(\phi\right)-\beta\psi^{2}\right)=0,$
(18) $\displaystyle\ddot{\phi}+3H\dot{\phi}-\psi$ $\displaystyle=0,$ (19)
$\displaystyle\beta\left(\ddot{\psi}+3H\dot{\psi}\right)+\frac{1}{2}\left(\psi+V_{,\phi}\right)$
$\displaystyle=0.$ (20)

where $H=\frac{\dot{a}}{a}$ is the Hubble function.

For $\omega=1$, the scalar field is a quintessence while for $\omega=-1$
scalar field $\psi$ is a phantom field. However, as it was found before in
angup1 , the equation of state parameter for the effective parameter can cross
the phantom divide line and for the case of a quintessence field, that is,
because the second scalar field $\psi$ can dominates such that $w_{eff}<-1$.
Last but not least, from the analysis of the dynamics in angup2  it was found
that the de Sitter universe is always a late-time attractor independent from
the nature of the scalar field potential.

Furthermore, the field equations (17)-(20) can derive from the variation of
the point-like Lagrangian

$L_{Q}=3a\dot{a}^{2}-\frac{1}{2}a^{3}\dot{\phi}^{2}-2\beta
a^{3}\dot{\phi}\dot{\psi}+a^{3}V\left(\phi\right)-\beta a^{3}\psi^{2}.$ (21)

###  III.2 Brans-Dicke cosmology

Brans-Dicke model is recovered when $F\left(\phi\right)=\phi^{2}$, where
$\omega=\bar{\omega}_{BD}$ is now the Brans-Dicke parameter. In this case the
cosmological field equations in the case of a FLRW spacetime

$6H^{2}+12H\left(\frac{\dot{\phi}}{\phi}\right)+\frac{\bar{\omega}_{BD}}{2}\left(\dot{\phi}^{2}+2\beta\hbar^{2}\left(\frac{\dot{\phi}}{\phi}\right)\left(\frac{\dot{\psi}}{\phi}\right)-\frac{\psi^{2}}{\phi^{2}}\right)+\frac{V\left(\phi\right)}{\phi^{2}}=0,$
(22)
$2\dot{H}+3H^{2}+4\left(\frac{\dot{\phi}}{\phi}\right)H+\left(\frac{\bar{\omega}_{BD}}{4}-2\right)\left(\frac{\dot{\phi}}{\phi}\right)^{2}+\frac{V\left(\phi\right)}{2\phi^{2}}+\beta\frac{\bar{\omega}_{BD}}{4\phi^{2}}\left(\psi^{2}-2\dot{\psi}\dot{\phi}\right)=0,$
(23)
$\bar{\omega}_{BD}\left(\beta\ddot{\psi}+\ddot{\phi}\right)+3\bar{\omega}_{BD}H\left(\beta\dot{\psi}+\dot{\phi}\right)+V_{,\phi}+12\phi\left(\dot{H}+2H^{2}\right)=0.$
(24) $\ddot{\phi}+3H\dot{\phi}+\psi=0,$ (25)

Similarly as before the Lagrangian function which reproduces the field
equations reads

$L\left(a,\dot{a},\phi,\dot{\phi},\psi,\dot{\psi}\right)=-6a\phi^{2}\dot{a}^{2}-12\phi
a^{2}\dot{a}\dot{\phi}-\frac{\bar{\omega}_{BD}}{2}a^{3}\left(\dot{\phi}^{2}+2\beta\hbar^{2}\dot{\phi}\dot{\psi}-\psi^{2}\right)+a^{3}V\left(\phi\right).$
(26)

In angup4  it was found that the modified field equations provide a different
cosmological history. The GUP introduces significant changes to the dynamics
and asymptotic solutions of the field equations, resulting in newfound degrees
of freedom. These modifications diverge from those of the unmodified model.
Moreover, the physical characteristics of the asymptotic solutions are
contingent upon the exponent of the potential function rather than the value
of the Brans-Dicke parameter, which is the case in the unmodified model.

##  IV $F\left(R,\square R\right)$-gravity

In this Section we incoporate Lagrange multipliers in order to study the
scalar-tensor description of $F\left(R,\square R\right)$-gravity. In our
study, we focus into the scenario where scalar field models with GUP are
derived from $F\left(R,\square R\right)$-gravity.

###  IV.1 $f\left(R\right)$-gravity

The Ricciscalar $R$ is defined by second-order derivatives of the metric
tensor; hence a Lagrangian which is an arbitrary function of $R$ leads to
fourth-order equations of motions. Indeed, $f\left(R\right)$-gravity is a
fourth-order theory and equivalent with a special case of the scalar-tensor
theory.

Consider a Riemannian manifold with metric $g_{\mu\nu}$ and Ricciscalar $R$.
The Action Integral of the $f\left(R\right)$-gravity is as follows

$S=\int dx^{4}\sqrt{-g}f\left(R\right)$ (27)

when $f\left(R\right)=R-2\Lambda$ the limit of General Relativity is
recovered.

We make use of the Lagrange multiplier $\lambda$, thus the Action Integral
(27) reads

$S=\int
dx^{4}\sqrt{-g}\left(f\left(\chi\right)+\lambda\left(\chi-R\right)\right).$
(28)

The equation of motion for $\lambda$ is $\frac{\partial S}{\partial\chi}=0$,
that is, $\lambda=-f_{,\chi}$.

By replacing $\lambda$ in (28) we find

$S=\int dx^{4}\left[f_{,\chi}R+\left(f-\chi f_{,\chi}\right)\right],$ (29)

or equivalently

$S=\int dx^{4}\left[\varphi
R+V\left(\varphi\right)\right],~{}\varphi=f_{,\chi},$ (30)

and potential function $V\left(\varphi\right)=\left(f-\chi
f_{,\chi}\right)\,$.

###  IV.2 $F\left(R,\square R\right)$-gravity

We consider the sixth-order theory of gravity with Action Integral

$S=\int dx^{4}\sqrt{-g}F\left(R,\square R\right)$ (31)

where $\square R=g^{\mu\nu}R_{;\mu\nu}$.

We consider the new variables $\chi,\zeta$ and we introduce the Lagrange
multipliers $\lambda_{1},~{}\lambda_{2}$. Therefore, the Action Integral (31)
is expressed as follows

$S=\int
dx^{4}\sqrt{-g}\left[F\left(\chi,\zeta\right)+\lambda_{1}\left(\chi-R\right)+\lambda_{2}\left(\zeta-\square
R\right)\right].$ (32)

Variation with respect to the Lagrange multipliers provide the equations of
motions $\lambda_{1}+F_{,\chi}=0$ and $\lambda_{2}+F_{,\zeta}=0$.

We replace in the gravitational Action Integral (32) and we find

$S=\int
dx^{4}\sqrt{-g}\left[\left(F_{,\chi}+F_{,\zeta}\square\right)R+\left(F-\chi
F_{,\chi}-\zeta F_{,\zeta}\right)\right]$ (33)

Integration by parts of the second term of (33) gives

$\int dx^{4}\sqrt{-g}\left(F_{,\zeta}\square R\right)=-\int
dx^{4}\sqrt{-g}\left(F_{,\zeta\zeta}g^{\mu\nu}\chi_{;\mu}\zeta_{;\nu}\right).$

Consequently we end with the Gravitational Action Integral

$S=\int
dx^{4}\sqrt{-g}\left[F_{,\chi}R-F_{,\zeta\zeta}g^{\mu\nu}\nabla_{\mu}\chi\nabla_{\nu}\zeta+\left(F-\chi
F_{,\chi}-\zeta F_{,\zeta}\right)\right]$ (34)

As we expected we have two scalar fields, the $\left\\{\chi,\zeta\right\\}$,
because $F\left(R,\square R\right)$-gravity is a sixth-order gravity. The
field $\chi$ is a non-minimally coupled field and the field $\zeta$ is an
extra field which can be seen as perturbation effects. What is of special
interest is the non-diagonal term
$F_{,\zeta\zeta}g^{\mu\nu}\chi_{;\mu}\zeta_{;\nu}$ which has similarities with
the non-diagonal term in (15).

####  IV.2.1 $F\left(R,\square R\right)=f\left(R\right)+K\left(\square
R\right)$

Assume now the case where $F\left(R,\square R\right)$ is a separable function,
i.e. $F_{,\chi\zeta}=0$. Hence, by replacing $F\left(R,\square
R\right)=f\left(R\right)+K\left(\square R\right)$ in (34) we derive

$S=\int
dx^{4}\sqrt{-g}\left[M\left(\chi\right)R-g^{\mu\nu}\chi_{;\mu}\psi_{;\nu}+V\left(\chi\right)+\hat{V}\left(\psi\right)\right]$
(35)

where $\psi=K_{,\zeta}$, $M\left(\chi\right)=f_{,\chi}\left(\chi\right),$
$V\left(\chi\right)=f\left(\chi\right)-\chi f_{,\chi}\left(\chi\right)$ and
$\hat{V}\left(\psi\right)=K\left(\zeta\right)-\zeta
K_{,\zeta}\left(\zeta\right)$.

We observe that the Action Integral (35) is of the form of the GUP scalar-
tensor theory (15) for zero value of parameter $\omega,~{}$that is,
$\omega=0$. Moreover, for $K\left(\square
R\right)=\frac{1}{2\hat{\beta}}\left(\square R\right)^{2}$, we find
$\hat{V}\left(\psi\right)=-\frac{\hat{\beta}}{2}\psi^{2}$, which means that
that the separable $F\left(R,\square
R\right)=f\left(R\right)+\frac{1}{2\hat{\beta}}\left(\square R\right)^{2}$
introduce new geometrodynamical terms in the field equations which can be
attributed to the corrections provided by the quadratic GUP in scalar-tensor
theory. In our case, the quantum corrections of $f\left(R\right)$-gravity are
on the limit of O’Hanlon theory where $\omega=0$. Parameter $\hat{\beta}$ is
the analogue of the deformed parameter $\beta$ of GUP.

In the special case where $f\left(R\right)$ is a linear function and for the
$F\left(R,\square R\right)=R-2\Lambda+\frac{1}{2\alpha}\left(\square
R\right)^{2}$ theory, from (35) it follows

$S=\int
dx^{4}\sqrt{-g}\left[R-2\Lambda-\hat{\beta}g^{\mu\nu}\chi_{;\mu}\psi_{;\nu}-\frac{\hat{\beta}}{2}\psi^{2}\right].$
(36)

##  V Cosmological solutions

We proceed with the analysis of cosmological evolution for the physical
parameters for the gravitational Action Integral (35) with background geometry
the FLRW spacetime (16).

For the spatially flat FLRW geometry, the field equations can be derived from
the variation of the point-like Lagrangian function

$L\left(a,\dot{a},\chi,\dot{\chi},\psi,\dot{\psi}\right)=6aM\left(\chi\right)\dot{a}^{2}+6a^{2}M_{,\phi}\left(\chi\right)\dot{a}\dot{\chi}+a^{3}\dot{\chi}\dot{\psi}+a^{3}\left(V\left(\chi\right)+\hat{V}\left(\psi\right)\right),$
(37)

in which the first modified Friedmann equation is the constraint equation

$6M\left(\chi\right)H^{2}+6M_{,\chi}\left(\chi\right)H\dot{\chi}+\dot{\chi}\dot{\psi}-\left(V\left(\chi\right)+\hat{V}\left(\psi\right)\right)=0,$
(38)

which can be seen as the conservation laws of energy for the three-dimensional
dynamical system described by the point-like Lagrangian (37). The rest of the
field equations are the Euler-Lagrange equations
$\left(\frac{d}{dt}\frac{\partial}{\partial\mathbf{\dot{y}}}-\frac{\partial}{\partial\mathbf{y}}\right)L\left(\mathbf{y,\dot{y}}\right)=0$,
where now $\mathbf{y=}\left(a,\chi,\psi\right)$.

Consequently, the second-order differential equations are

$2M\left(2\dot{H}+3H^{2}\right)+2M_{,\chi}\left(2H\dot{\chi}+\ddot{\chi}\right)+2\dot{\chi}^{2}M_{,\chi\chi}-\dot{\chi}\dot{\psi}-\left(V\left(\chi\right)+\hat{V}\left(\psi\right)\right)=0,$
(39) $\ddot{\chi}+3H\dot{\chi}-\hat{V}\left(\psi\right)_{,\psi}=0,$ (40)
$\ddot{\psi}+3H\dot{\psi}+6M_{,\chi}\left(\dot{H}+2H^{2}\right)-V_{,\chi}=0.$
(41)

Let us focus in the case where $M\left(\chi\right)=1$ and
$V\left(\chi\right)=0$. Then the latter gravitational field equations read

$3H^{2}+\frac{1}{2}\left(\dot{\chi}\dot{\psi}-\hat{V}\left(\psi\right)\right)=0,$
(42)
$2\dot{H}+3H^{2}-\frac{1}{2}\left(\dot{\chi}\dot{\psi}+\hat{V}\left(\psi\right)\right)=0,$
(43) $\ddot{\chi}+3H\dot{\chi}-\hat{V}\left(\psi\right)_{,\psi}=0,$ (44)
$\ddot{\psi}+3H\dot{\psi}+6\left(\dot{H}+2H^{2}\right)=0.$ (45)

Therefore, the effective energy density and pressure components are defined as

$\rho_{eff}=\frac{1}{2}\left(-\dot{\chi}\dot{\psi}+\hat{V}\left(\psi\right)\right),$
(46)
$p_{eff}=-\frac{1}{2}\left(\dot{\chi}\dot{\psi}+\hat{V}\left(\psi\right)\right),$
(47)

where the effective equation of state parameter is written

$w_{eff}=\frac{\left(\dot{\chi}\dot{\psi}+\hat{V}\left(\psi\right)\right)}{\left(\dot{\chi}\dot{\psi}-\hat{V}\left(\psi\right)\right)}.$
(48)

Indeed, when the potential function $\hat{V}\left(\psi\right)$ dominates the
asymptotic solution is that of the de Sitter universe with $w_{eff}=-1$.

If we define the new scalar fields
$\chi=\phi+\varpi,~{}\dot{\psi}=\phi-\varpi$, then equation (42) reads

$3H^{2}+\frac{1}{2}\left(\left(\dot{\phi}^{2}-\dot{\varpi}^{2}\right)-\hat{V}\left(\phi-\varpi\right)\right)=0,$
(49)

which means that the model is equivalent to the quintom cosmological model
with a mixed potential term, which is different from the analysis presented in
quin11 . For a complete analysis of the dynamics in quintom cosmology we refer
the reader in quin22 .

Below we present the phase-space analysis for the field equations (42)-(45) by
using dimensionless variables in the $H$-normalization approach q5 .

###  V.1 Phase-space analysis

We define the new dimensionless variables

$x=\frac{\dot{\chi}}{H}~{},~{}z=\frac{\dot{\psi}}{6H}\text{
},~{}w=\frac{\hat{V}\left(\psi\right)}{6H^{2}}~{},~{}\lambda=\frac{\hat{V}_{,\psi}}{\hat{V}},~{}\tau=\ln
a$ (50)

and we rewrite the field equations (42)-(45) in the form of the subsequent
algebraic-differential system

$\displaystyle\frac{dx}{d\tau}$
$\displaystyle=\frac{3}{2}\left(w\left(4\lambda-x\right)-x\left(1+x+z\right)\right),$
(51) $\displaystyle\frac{dz}{d\tau}$
$\displaystyle=-\frac{3}{2}z\left(1+w+xz\right),$ (52)
$\displaystyle\frac{dw}{d\tau}$
$\displaystyle=3w\left(1-w-\left(x-2\lambda\right)z\right),$ (53)
$\displaystyle\frac{d\lambda}{d\tau}$
$\displaystyle=6\lambda^{2}z\left(\Gamma\left(\lambda\right)-1\right)~{},~{}\Gamma\left(\lambda\right)=\frac{V_{\psi\psi}V}{\left(V_{,\psi}\right)^{2}}\text{\thinspace},$
(54)

with constraint

$1+xz-w=0.$ (55)

By utilizing the constraint equation, we can effectively reduce the dimension
of the aforementioned dynamical system by one.

Hence, we arrive at the three-dimensional system

$\displaystyle\frac{dx}{d\tau}$
$\displaystyle=3\left(2\lambda-x\right)\left(1+xz\right),$ (56)
$\displaystyle\frac{dz}{d\tau}$ $\displaystyle=-3z\left(1+xz\right),$ (57)
$\displaystyle\frac{d\lambda}{d\tau}$
$\displaystyle=6\lambda^{2}z\left(\Gamma\left(\lambda\right)-1\right)\text{.}$
(58)

At each stationary point
$P=\left(x\left(P\right),z\left(P\right),\lambda\left(P\right)\right)~{}$of
the dynamical system (56)-(58) describes an asymptotic solution where the
effective fluid source has the equation of state parameter

$w_{eff}=-1-2x\left(P\right)z\left(P\right)\text{.}$ (59)

###  V.2 Exponential potential

In order to reduce further the dimension of the dynamical system we consider
the simple case where $\hat{V}\left(\psi\right)$ is the exponential potential,
that is, $\hat{V}\left(\psi\right)=V_{0}e^{\lambda_{0}\psi}$. For this
potential we calculate $\Gamma\left(\lambda\right)=1$, and
$\lambda=\lambda_{0}$ is always a constant.

Hence, the stationary points of the dynamical system (56), (57) are

$P_{1}=\left(x_{1},-\frac{1}{x_{1}}\right)~{},~{}P_{2}=\left(2\lambda,0\right)\text{.}$

For the family of points $P_{1}$ we derive $w_{eff}\left(P_{1}\right)=1$, this
implies that the points represent a family of stiff fluid solutions.
Furthermore, the asymptotic solution at point $P_{2}$ describes the
accelerated de Sitter universe, because $w_{eff}\left(P_{2}\right)=-1$.

Let us now proceed with the investigation of the stability properties of the
points. For the linearization system (56), (57) around the stationary points
$P_{1}$ we determine the eigenvalues
$e_{1}\left(P_{1}\right)=6\left(1-\frac{\lambda}{x_{1}}\right)$,
$e_{2}\left(P_{1}\right)=0$. For $1-\frac{\lambda}{x_{1}}>0$ points $P_{1}$
are sources, however for $1-\frac{\lambda}{x_{1}}<0$ we will infer about the
stability of the points from the phase-space portraits. For point $P_{2}$ the
eigenvalues of the linearized system are $e_{1}\left(P_{2}\right)=-3$ and
$e_{2}\left(P_{2}\right)=-3$, which means that the de Sitter universe is a
future attractor for the dynamical system.

In Fig. 1 we present the phase-space portrait for the two-dimensional
dynamical system (56), (57). We observe that the family of points $P_{1}$
describe always unstable solutions and they are the boundaries where the
trajectories move to the infinity.

Figure 1: Phase-space portrait for the two dynamical system (56), (57) for
$\lambda=1$ (left fig.) and $\lambda=-1$ (right fig.). With red lines are the
family of points $P_{1}$ and with red dot is marked the de Sitter attractor
$P_{2}$. We observe that points $P_{1}$ describe always unstable solutions and
are the limits where the trajectories move to the infinity regime.

####  V.2.1 Analysis at infinity

We define the Poincare variables

$x=\frac{X}{\sqrt{1-X^{2}-Z^{2}}}~{},~{}z=\frac{Z}{\sqrt{1-X^{2}-Z^{2}}}$

and the new independent variable $dT=\sqrt{1-X^{2}-Z^{2}}d\tau$.

In the Poincare variables the field equations read

$\displaystyle\frac{dX}{dT}$
$\displaystyle=-3\left(X^{2}+Z^{2}-1-XZ\right)\left(2\lambda-X\left(2\lambda
X+\sqrt{1-X^{2}-Z^{2}}\right)\right),$ (60) $\displaystyle\frac{dZ}{dT}$
$\displaystyle=3Z\left(X^{2}+Z^{2}-1-XZ\right)\left(2\lambda
X+\sqrt{1-X^{2}-Z^{2}}\right),$ (61)

and the $w_{eff}\left(x,z\right)$ becomes

$w_{eff}\left(X,Z\right)=-1-\frac{2XZ}{1-X^{2}-Z^{2}}\,.$ (62)

Hence, in order the solutions at the infinity to be physical accepted it
should be $XZ\geq 0$. For $XZ>0$ at infinity the solutions describe Big Rip
singularities with $w_{eff}\rightarrow-\infty$, while in the limit $XZ=0$, the
de Sitter universe is recovered.

The stationary points at the infinity are calculated $Q_{1}^{\pm}=\left(\pm
1,0\right)$ and $Q_{2}^{\pm}=\left(0,\pm 1\right)$, which means that
$w_{eff}\left(Q_{1}^{\pm}\right)=-1$ and $w_{eff}\left(Q_{1}^{\pm}\right)=-1$.
The linearized system at points $Q_{1}^{\pm},Q_{2}^{\pm}$ are
$e_{1}\left(Q_{1}^{\pm}\right)=0$, $e_{2}\left(Q_{1}^{\pm}\right)=0$ and
$e_{1}\left(Q_{2}^{\pm}\right)=\pm 6\lambda_{2}$,
$e_{2}\left(Q_{2}^{\pm}\right)=0$ respectively. From the phase-space portrait
of Fig. 2, we remark that the stationary points at infinity describe always
unstable solutions and the unique attractor of the dynamical system is the de
Sitter solution described by the stationary point $P_{2}$. The regions outside
the red lines (points $P_{1}$) lead to unphysical solutions with
$w_{eff}\left(X,Z\right)>1,$ thus the accepted initial conditions of the
dynamical system are those inside the red lines.

Figure 2: Phase-space portrait for the two dynamical system (56), (57) for
$\lambda=1$ (left fig.) and $\lambda=-1$ (right fig.) in the Poincare
variables $\left(X,Z\right)$. With red lines are the family of points $P_{1}$
and with red dot are marked the de Sitter attractor $P_{2}$ and the stationary
points at the infinity $Q_{1}^{\pm}$ and $Q_{2}^{\pm}.$ The regions outside
the red lines (points $P_{1}$) lead to unphysical solutions with $w_{eff}>1,$
thus the accepted initial conditions of the dynamical system are those inside
the red lines.

##  VI Conclusions

Scalar field gravitational theories which are modified by the presence of a
minimum length in the Uncertainty Principle are equivalent with two-scalar
field theories of gravity. In this work we investigate the equivalent modified
theory of gravity which can be described by the latter GUP scalar field
theory. We assumed $F\left(R,\square R\right)$-gravitational model which is a
sixth-order theory. We found that the separable model $F\left(R,\square
R\right)=f\left(R\right)+K\left(\square R\right)$ can be in comparison with
the GUP scalar field theory. Specifically, for the $K\left(\square
R\right)=\frac{1}{2\hat{\beta}}\left(\square R\right)^{2}$ model, there are
introduced geometrodynamical terms in the field equations which are in the
same form with that introduced by the deformed Heisenberg algebra of the
quadratic GUP. Indeed, parameter $\hat{\beta}$ is linearly related with the
deformation parameter $\beta$ of GUP.

We focused in the case of $F\left(R,\square R\right)=R+K\left(\square
R\right)$ where now the $K\left(\square R\right)$-terms can attribute the
quantum corrections in the limit of General Relativity. For the latter theory
and for a spatially flat FLRW background geometry we investigated the
cosmological dynamics. From the analysis of the phase-space it was found that
there exist always a future attractor which corresponds to the de Sitter
solution without to introduce a cosmological constant term. We conclude that
the quantum corrections related by the $K\left(\square R\right)$ can drive the
dynamics such that inflation to occurred.

In a future work we plan to investigate further the effects of the
introduction of higher-order terms $K\left(\square R\right)$ in
$f\left(R\right)$f theory for the cosmological evolution and investigate the
case of EUP in gravitational models.

Data Availability Statements: Data sharing is not applicable to this article
as no datasets were generated or analyzed during the current study.

###### Acknowledgements.

The author thanks the support of Vicerrectoría de Investigación y Desarrollo
Tecnológico (Vridt) at Universidad Católica del Norte through Núcleo de
Investigación Geometría Diferencial y Aplicaciones, Resolución Vridt No -
098/2022.

## References

  * (1) A.G. Riess et al., Astron. J. 116, 1009 (1998)
  * (2) M. Tegmark et al., Astrophys. J. 606, 702 (2004)
  * (3) M. Kowalski et al., Astrophys. J. 686, 749 (2008)
  * (4) E. Komatsu et al., Astrophys. J. Suppl. Ser. 180, 330 (2009)
  * (5) N. Suzuki et al., Astrophys. J. 746, 85 (2012)
  * (6) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989)
  * (7) T. Padmanabhan, Phys. Rept. 380, 235 (2003)
  * (8) L. Perivolaropoulos and F. Skara, New Astronomy Reviews, 95, 101659 (2022)
  * (9) P. Ratra and L. Peebles, Phys. Rev. D 37, 3406 (1988)
  * (10) E.J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006)
  * (11) C. Rubano and J.D. Barrow, Phys. Rev. D 64, 127301 (2001)
  * (12) S. Basilakos and G. Lukes-Gerakopoulos, Phys. Rev. D 78, 083509 (2008)
  * (13) R.C. Caldwell, M. Kamionkowski and N.N. Weinberg, Phys. Rev. Lett. 91, 071301 (2003)
  * (14) F. Briscese, E. Elizalde, S. Nojiri and S.D. Odintsov, Phys. Lett. B 646, 105 (2007)
  * (15) J. Khoury and A. Wetlman, Phys. Rev. D 69, 044026 (2004)
  * (16) J.E.M. Aguilar, M. Montes and A. Bernal, Phys. Scr. 98, 035021 (2023)
  * (17) C.H.-T. Wang, J.A. Reid, A.St.J. Murphy, D. Rodrigues, M. Al Alawi, R. Bingham, J.T. Mendonca and T.B. Davies, Phys. Lett. A 380, 3761 (2016)
  * (18) A. Zucca, L. Pogosian, A. Silvestri, Y. Wang and G.-B. Zhao, Phys. Rev. D 101, 043518 (2020)
  * (19) G. Koutsoumbas, K. Ntrekis, E. Papantonopoulos and E.N. Saridakis, JCAP 18, 003 (2018)
  * (20) R. Gannouji and M. Sami, Phys. Rev. D 82, 024011 (2010)
  * (21) A. Paliathanasis, Class. Quantum Grav. 37, 195014 (2020)
  * (22) G. Leon, A. Paliathanasis and J.L. Morales-Martinez, Eur. Phys. J. C 78, 753 (2018)
  * (23) T. Pati, S. Panda, M. Sharma and Ruchika, Eur. Phys. J. C 83, 131 (2023)
  * (24) A. Al Mamon, A. Paliathanasis and S. Saha, l, Eur. Phys. J. C 82, 232 (2021)
  * (25) S. Basilakos, G. Lukes-Gerakopoulos, Phys. Rev. D 78, 083509 (2008)
  * (26) H. Li, W. Yang and L. Gai, A&A 623, A28 (2019)
  * (27) V. A. Popov, Phys. Lett. B 686, 211 (2010)
  * (28) R. von Marttens, D. Babosa and J. Alcaniz, JCAP 04, 052 (2023)
  * (29) C. Deffayet, X. Gao, D.A. Steer and G. Zahariade, Phys. Rev. D 84, 064039 (2011)
  * (30) B. Chen, Y. Wu, J. Chi, W. Liu and Y. Hu, Universe 8, 520 (2022)
  * (31) A. Kundu, S.D. Pathak and V.K. Ojha, Comm. Theor. Phys. 73, 025402 (2021)
  * (32) J.S.Bagla, H.K.Jassal and T.Padmanabhan, Phys. Rev. D 67, 063504 (2003)
  * (33) E. Nouri, H. Motavalli and A.R. Akbarieh, Int. J. Mod. Phys. D 31, 2250020 (2022)
  * (34) S. Nojiri, S.D. Odintsov and V.K. Oikonomou, Phys.Rept. 692, 1 (2017)
  * (35) A. Joyce, L. Lombriser and F. Schmidt, Annu. Rev. Nucl. Part. Sci. 66, 122 (2016)
  * (36) T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, Phys. Rept. 513, 1 (2012)
  * (37) S. Capozziello and V. Faraoni, Gen. Relativ. Grav. 40, 357 (2008)
  * (38) H.A. Buchdahl, Mon. Not. Roy. Astron. Soc. 150, 1 (1970)
  * (39) R. Ferraro and F. Fiorini, Phys. Rev. D 75, 084031 (2007)
  * (40) J. B. Jimenez, L. Heisenberg and T. Koivisto, Phys. Rev. D 98, 044048 (2018)
  * (41) G. Cognola, E. Elizalde, S. Nojiri, S.D. Odintsov and S. Zerbini, Phys. Rev. D 73, 084007 (2006)
  * (42) G. Kofinas and E.N. Saridakis, Phys. Rev. D 90, 084044 (2014)
  * (43) T. Harko, F.S.N. Lobo, G. Otalora and E.N. Saridakis, JCAP 12, 021 (2014)
  * (44) A. Paliathanasis and G. Leon, Eur. Phys. J. Plus 136, 1092 (2021)
  * (45) H.R. Kausar, R. Saleem and A. Ilyas, Phys. Dark Universe 26, 100401 (2019)
  * (46) G.J. Olmo, Int. J. Mod. Phys. D 20, 413 (2011)
  * (47) A. De Felice, J.-M. Gerard and T. Suyama, Phys. Rev. D 82, 063526 (2010)
  * (48) T.P. Sotiriou and V. Faraoni Rev. Mod. Phys. 82 451 (2010)
  * (49) S. Nojiri and S.D. Odintsov, Phys. Rep. 505 59 (2011)
  * (50) V. K. Oikonomou and I. Giannakoudi, Int. J. Mod. Phys. D 31, 2250075 (2022)
  * (51) S.P. Hatkar, P.S. Dudhe and S.D. Katore, Foundations in Physics 49, 1067 (2019)
  * (52) S. Nojiri and S. Odintsov, Int. J. Geom. Meth. Mod. Phys. 11, 1460006 (2014)
  * (53) A. de la Cruz-Dombriz, A. Dobado and A. L. Maroto, Phys. Rev. D 80, 124011 (2009)
  * (54) Kh Jafarzade, E Rezaei and S H Hendi, Progress of Theoretical and Experimental Physics 2023, 053E01 (2023)
  * (55) G. G. L. Nashed and S. Nojiri, Phys. Rev. D 102, 124022 (2020)
  * (56) H. Nariai and K. Tomita, Prog. Theor. Phys. 46, 776 (1971)
  * (57) G.V. Bicknell, J. Phys. A.: Math. Nucl. Gen. 7, 1061 (1974)
  * (58) J.D. Barrow, Nucl. Phys. B 296, 679 (1988)
  * (59) P.A.R. Ade et al. (Planck 2015 Collaboration), A.&A. 594, A20 (2016)
  * (60) A.A. Starobinsky, Phys. Lett. B 91, 99 (1980)
  * (61) H. Nariai, Prog. Theor. Phys. 46, 433 (1971)
  * (62) F.C. Michel, Ann. Phys. 76, 281 (1973)
  * (63) M. Maggiore, Phys. Lett. B 304, 65 (1993).
  * (64) R. Casadio and F. Scardigli, Phys. Lett. B 807, 135558 (2020)
  * (65) S. Benczik, L.N. Chang, D. Minic, N. Okamura, S. Rayyan and T. Takeuchi, Phys. Rev.D 66, 026003 (2002)
  * (66) A.F. Ali, Class. Quantum Grav. 28, 065013 (2011)
  * (67) E.C. Vagenas, A.F. Ali, M. Hemeda and H. Alshal, Eur. Phys. J. C 79, 398 (2019)
  * (68) A.F. Ali, S. Das and E.C. Vagenas, Phys. Lett. B 678, 497 (2009)
  * (69) Z. Silagadge, Phys. Lett. A 373, 2643 (2009)
  * (70) Y.-G. Miao and Y.-J. Zhao, Int. J. Phys. D. 23, 1450062 (2014)
  * (71) L.N. Chang, D. Minic, N. Okamura and T. Takeuchi, Phys. Rev D 65, 125028 (2002)
  * (72) C. Bambi and F.R. Urban, Class. Quantum Grav. 25, 095006 (2008)
  * (73) S. Aghababaei, H. Moradpour and E.C. Vagenas, Eur. Phys. J. Plus 136, 997 (2021)
  * (74) P. Bosso, G.G. Luciano, L. Petruzziello and F. Wagner, 30 years in: Quo vadis generalized uncertainty principle? (2023) [arXiv:2305.16193]
  * (75) A. Paliathanasis, S. Pan and S. Pramanik, Class. Quant. Grav. 32, no.24, 245006 (2015)
  * (76) A. Giacomini, G. Leon, A. Paliathanasis and S. Pan, Eur. Phys. J. C 80, 931 (2020)
  * (77) A. Paliathanasis, G. Leon, W. Khyllep, J. Dutta and S. Pan, Eur. Phys. J. C 81, 607 (2021)
  * (78) A. Paliathanasis and G. Leon, Gen. Relativ. Grav. 55, 12 (2023)
  * (79) J. O’Hanlon, Phys. Rev. Lett. 29, 137 (1972)
  * (80) V. Faraoni, Cosmology in Scalar-Tensor Gravity. Dordrecht, Kluwer Academic (2004)
  * (81) C.H Brans and R.H. Dicke, Phys. Rev. 124, 925 (1961)
  * (82) C. Gao, Y. Gan, X. Wang and X. Chen, Phys. Lett. B 702, 107 (2011)
  * (83) S. Nojiri, S.D. Odintsov, V.K. Oikonomou and T. Paul, Phys. Rev. D 100, 084056 (2019)
  * (84) D. F. Jimenez, L. N. Granda, and E. Elizalde, Int. J. Mod. Phys. D 28, 1950171 (2019)
  * (85) S. Gottlober, H.-J. Schmidt and A.A. Starobinsky, Class. Quamtum Grav. 7, 893 (1990)
  * (86) D. Wands, Class. Quantum Grav. 11, 269 (1994)
  * (87) A.L. Berkin and K.-i. Meida, Phys. Lett. B 245, 348 (1990)
  * (88) A.R.R. Castellanos, F. Sobreira, I.L. Shapiro and A.A. Starobinsky, JCAP 12, 007 (2018)
  * (89) L. Amendola, A.B. Mayer, S. Capozziello, F. Occhionero and S. Gottlober, Class. Quantum Grav. 10, L43 (1993)
  * (90) S. Carloni, J.L. Rosa and J.P.S. Lemos, Phys. Rev. D 99, 104001 (2019)
  * (91) I. Quandt and H.-J. Schmidt, Astron. Nachr. 312, 97 (1991)
  * (92) S. Nesseris and A. Mazumbar, Phys. Rev. D 79, 104006 (2009)
  * (93) M. Skugoreva, A. Toporensky and P. Tretyakov, Grav. Cosmol. 17, 110 (2011)
  * (94) M. Iihoski, JCAP 02, 022 (2011)
  * (95) S. Rani, M.B.A. Sulehri and A. Jawad, Phys. Dark Univ. 29, 100555 (2020)
  * (96) S. Das and E.C. Vagenas, Phys. Rev. Lett. 101, 221301 (2008)
  * (97) P. Jizba, H. Kleinert and F. Scardigli, Phys. Rev. D 81, 084030 (2010)
  * (98) L. Buoninfante, G.G. Luciano and L. Petruzziello, Eur. Phys. J. C 79, 663 (2019).
  * (99) L. Buoninfante, G. Lambiase, G.G. Luciano and L. Petruzziello, Eur. Phys. J. C 80, 853 (2020)
  * (100) S. Giardino and V. Salzano, Eur. Phys. J. C 81, 110 (2021)
  * (101) A. Kempf, J. Phys. A: Math. Gen. 30, 2093 (1997)
  * (102) H. Hinrichen and A. Kempf, J. Math. Phys. 37, 2121 (1996)
  * (103) S.K. Moayedi, M.R. Setare and H Moayeri, Int. J. Theor. Phys. 49, 2080 (2010)
  * (104) A.N. Tikhonov, Sbornic. Math. 31, 575 (1952)
  * (105) E.N. Saridakis, J.M. Weller, Phys. Rev. 81, 123523 (2010)
  * (106) G. Leon, A. Paliathanasis and J.L. Morales-Martinez, Eur. Phys. J. C 75, 753 (2018)