Consider some $P_{i}$. $P_{i}$ is a lexicographically shortest path in $G$. We
have $|P_{i}|\leq h$ with high probability, since $H\subseteq C$ hits each of
the $O(n^{2})$ lexicographically shortest paths that have more than $h$ hops
in every of $O(\operatorname{poly}{n})$ versions of $G$ independent of $H$. As
a result, we have

$\delta_{G}(u_{i},v_{i})=|P_{i}|=\delta_{G}^{h}(u_{i},v_{i}).$

We now prove that $w_{X}(u_{i}v_{i})=|P_{i}|$. This will imply that there
exists an $s\to t$ path in $X$ of weight
$\sum_{i}w_{X}(u_{i}v_{i})\leq\sum_{i}|P_{i}|=\delta_{G}(s,t)$. If
$(u_{i},v_{i})\in Y_{s}$, then $w_{X}(u_{i}v_{i})=\delta_{G}^{h}(u_{i},v_{i})$
by construction.

Suppose that $(u_{i},v_{i})\notin Y_{s}$. Then we simultaneously have
$u_{i}\in C\setminus D$, $v_{i}\in V\setminus D$, and either
$\pi_{u_{i},v_{i}}$ does not exist in $\Pi$ or
$\pi_{u_{i},v_{i}}\in\Pi\setminus\Pi_{D}$, i.e., $\pi_{u_{i},v_{i}}$ does not
pass through a vertex of $D$. Let $u_{i}=c_{j}\in C$, where $c_{j}$ and
$G_{j}\subseteq G_{0}$ are defined as in Theorem 5.1. Since $P_{i}$ is
internally disjoint with $C\cup D$, and $v_{i}\notin D$, $P_{i}$ contains no
edges inserted in the current phase. Hence, $P_{i}\subseteq G_{j}$. As a
result, ${\delta_{G_{j}}^{h}(u_{i},v_{i})\leq|P_{i}|<\infty}$. This proves
that $\pi_{u_{i},v_{i}}=\pi_{c_{j},v_{i}}$ exists in $\Pi$. But we also have
$\pi_{u_{i},v_{i}}\subseteq G_{j}-D\subseteq G_{0}-D\subseteq G-D$, so
$\delta_{G}(u_{i},v_{i})\leq|\pi_{u_{i},v_{i}}|$. Finally, by Theorem 5.1, we
obtain:

$\delta_{G}(u_{i},v_{i})\leq|\pi_{u_{i},v_{i}}|=\delta_{G_{j}}^{h}(u_{i},v_{i})\leq|P_{i}|=\delta_{G}(u_{i},v_{i}).$

This proves that $w_{G}(u_{i},v_{i})=\delta_{G}(u_{i},v_{i})=|P_{i}|$ also
when $(u_{i},v_{i})\notin Y_{s}$. ∎

We obtain the following combined worst-case update/query bound:

$\widetilde{O}\left(\frac{hn^{2+\rho}}{b}+\Delta
h^{2}n+\frac{n^{3+\rho}}{b\Delta}+\frac{n^{2}b}{h\Delta}+h^{2}n^{1+\rho}+\frac{n^{2}}{h}\right).$

Note that there is no point in setting $b<\Delta$, as increasing $b$ to
$\Delta$ would not make harm because of the $n^{2}/h$ term. So $b\geq\Delta$.
If we aim at equal update and query time, there is no point in putting
$b>\Delta$ either. For $b=\Delta$ we need to balance $\Delta h^{2}n$,
$\frac{n^{3+\rho}}{\Delta^{2}}$ and $n^{2}/h$. We get $\Delta=n/h^{3}$ and
finally we balance $n^{1+\rho}h^{6}$ and $n^{2}/h$. For $h=n^{(1-\rho)/7}$,
the update/query time is $O(n^{(13+\rho)/7})=O(n^{1.933})$. We have thus
proved.

######  Theorem 9.3.

Let $G=(V,E)$ be unweighted digraph. There exists a Monte Carlo randomized
data structure supporting fully dynamic single-edge updates to $G$ and
computing exact distances from an arbitrary query source $s\in V$ to all other
vertices in $G$ in $O(n^{1.933})$ worst-case time.

If $\omega=2$, the bound becomes $\widetilde{O}(n^{2-1/14})$.

## References

  * [ABF23] Amir Abboud, Karl Bringmann, and Nick Fischer.  Stronger 3-sum lower bounds for approximate distance oracles via additive combinatorics.  In Proceedings of the 55th Annual ACM Symposium on Theory of Computing, STOC 2023, Orlando, FL, USA, June 20-23, 2023, pages 391–404. ACM, 2023.
  * [ABKZ22] Amir Abboud, Karl Bringmann, Seri Khoury, and Or Zamir.  Hardness of approximation in p via short cycle removal: cycle detection, distance oracles, and beyond.  In STOC, pages 1487–1500. ACM, 2022.
  * [ACG12] Ittai Abraham, Shiri Chechik, and Cyril Gavoille.  Fully dynamic approximate distance oracles for planar graphs via forbidden-set distance labels.  In STOC, pages 1199–1218. ACM, 2012.
  * [ACK17] Ittai Abraham, Shiri Chechik, and Sebastian Krinninger.  Fully dynamic all-pairs shortest paths with worst-case update-time revisited.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2017, Barcelona, Spain, Hotel Porta Fira, January 16-19, pages 440–452. SIAM, 2017.
  * [ACT14] Ittai Abraham, Shiri Chechik, and Kunal Talwar.  Fully dynamic all-pairs shortest paths: Breaking the o(n) barrier.  In APPROX-RANDOM, volume 28 of LIPIcs, pages 1–16. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2014.
  * [AGMN92] Noga Alon, Zvi Galil, Oded Margalit, and Moni Naor.  Witnesses for boolean matrix multiplication and for shortest paths.  In 33rd Annual Symposium on Foundations of Computer Science, Pittsburgh, Pennsylvania, USA, 24-27 October 1992, pages 417–426. IEEE Computer Society, 1992.
  * [AHR+19] Bertie Ancona, Monika Henzinger, Liam Roditty, Virginia Vassilevska Williams, and Nicole Wein.  Algorithms and hardness for diameter in dynamic graphs.  In ICALP, volume 132 of LIPIcs, pages 13:1–13:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.
  * [AIMN91] Giorgio Ausiello, Giuseppe F. Italiano, Alberto Marchetti-Spaccamela, and Umberto Nanni.  Incremental algorithms for minimal length paths.  J. Algorithms, 12(4):615–638, 1991.
  * [AvdB23] Anastasiia Alokhina and Jan van den Brand.  Fully dynamic shortest path reporting against an adaptive adversary.  CoRR, abs/2304.07403, 2023.
  * [AW14] Amir Abboud and Virginia Vassilevska Williams.  Popular conjectures imply strong lower bounds for dynamic problems.  In FOCS, pages 434–443. IEEE Computer Society, 2014.
  * [AW21] Josh Alman and Virginia Vassilevska Williams.  A refined laser method and faster matrix multiplication.  In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 - 13, 2021, pages 522–539. SIAM, 2021.
  * [BC16] Aaron Bernstein and Shiri Chechik.  Deterministic decremental single source shortest paths: beyond the o(mn) bound.  In STOC, pages 389–397. ACM, 2016.
  * [BC17] Aaron Bernstein and Shiri Chechik.  Deterministic partially dynamic single source shortest paths for sparse graphs.  In SODA, pages 453–469. SIAM, 2017.
  * [Ber09] Aaron Bernstein.  Fully dynamic (2 + epsilon) approximate all-pairs shortest paths with fast query and close to linear update time.  In FOCS, pages 693–702. IEEE Computer Society, 2009.
  * [Ber16] Aaron Bernstein.  Maintaining shortest paths under deletions in weighted directed graphs.  SIAM J. Comput., 45(2):548–574, 2016.  Announced at STOC‘13.
  * [Ber17] Aaron Bernstein.  Deterministic partially dynamic single source shortest paths in weighted graphs.  In ICALP, volume 80 of LIPIcs, pages 44:1–44:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2017.
  * [BFN22] Jan van den Brand, Sebastian Forster, and Yasamin Nazari.  Fast deterministic fully dynamic distance approximation.  In 63rd IEEE Annual Symposium on Foundations of Computer Science, FOCS 2022, Denver, CO, USA, October 31 - November 3, 2022, pages 1011–1022. IEEE, 2022.
  * [BGS20] Aaron Bernstein, Maximilian Probst Gutenberg, and Thatchaphol Saranurak.  Deterministic decremental reachability, scc, and shortest paths via directed expanders and congestion balancing.  In FOCS, pages 1123–1134. IEEE, 2020.
  * [BGS21] Aaron Bernstein, Maximilian Probst Gutenberg, and Thatchaphol Saranurak.  Deterministic decremental SSSP and approximate min-cost flow in almost-linear time.  In FOCS, pages 1000–1008. IEEE, 2021.
  * [BGW20] Aaron Bernstein, Maximilian Probst Gutenberg, and Christian Wulff-Nilsen.  Near-optimal decremental SSSP in dense weighted digraphs.  In FOCS, pages 1112–1122. IEEE, 2020.
  * [BHG+21] Thiago Bergamaschi, Monika Henzinger, Maximilian Probst Gutenberg, Virginia Vassilevska Williams, and Nicole Wein.  New techniques and fine-grained hardness for dynamic near-additive spanners.  In SODA, pages 1836–1855. SIAM, 2021.
  * [BHS07] Surender Baswana, Ramesh Hariharan, and Sandeep Sen.  Improved decremental algorithms for maintaining transitive closure and all-pairs shortest paths.  J. Algorithms, 62(2):74–92, 2007.
  * [BN19] Jan van den Brand and Danupon Nanongkai.  Dynamic approximate shortest paths and beyond: Subquadratic and worst-case update time.  In 60th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2019, Baltimore, Maryland, USA, November 9-12, 2019, pages 436–455. IEEE Computer Society, 2019.
  * [BNS19] Jan van den Brand, Danupon Nanongkai, and Thatchaphol Saranurak.  Dynamic matrix inverse: Improved algorithms and matching conditional lower bounds.  In FOCS, pages 456–480. IEEE Computer Society, 2019.
  * [BR11] Aaron Bernstein and Liam Roditty.  Improved dynamic algorithms for maintaining approximate shortest paths under deletions.  In SODA, pages 1355–1365. SIAM, 2011.
  * [Bra] Jan van den Brand.  Complexity term balancer.  www.ocf.berkeley.edu/~vdbrand/complexity/.  Tool to balance complexity terms depending on fast matrix multiplication.
  * [BS19] Jan van den Brand and Thatchaphol Saranurak.  Sensitive distance and reachability oracles for large batch updates.  In FOCS, pages 424–435. IEEE Computer Society, 2019.
  * [Chu21] Julia Chuzhoy.  Decremental all-pairs shortest paths in deterministic near-linear time.  In STOC, pages 626–639. ACM, 2021.
  * [CK19] Julia Chuzhoy and Sanjeev Khanna.  A new algorithm for decremental single-source shortest paths with applications to vertex-capacitated flow and cut problems.  In STOC, pages 389–400. ACM, 2019.
  * [CS21] Julia Chuzhoy and Thatchaphol Saranurak.  Deterministic algorithms for decremental shortest paths via layered core decomposition.  In SODA, pages 2478–2496. SIAM, 2021.
  * [CZ21] Shiri Chechik and Tianyi Zhang.  Incremental single source shortest paths in sparse digraphs.  In SODA, pages 2463–2477. SIAM, 2021.
  * [CZ23] Shiri Chechik and Tianyi Zhang.  Faster deterministic worst-case fully dynamic all-pairs shortest paths via decremental hop-restricted shortest paths.  In Proceedings of the 2023 ACM-SIAM Symposium on Discrete Algorithms, SODA 2023, Florence, Italy, January 22-25, 2023, pages 87–99. SIAM, 2023.
  * [DFNV22] Michal Dory, Sebastian Forster, Yasamin Nazari, and Tijn de Vos.  New tradeoffs for decremental approximate all-pairs shortest paths.  CoRR, abs/2211.01152, 2022.
  * [DI01] Camil Demetrescu and Giuseppe F. Italiano.  Fully dynamic all pairs shortest paths with real edge weights.  In FOCS, pages 260–267. IEEE Computer Society, 2001.
  * [DI02] Camil Demetrescu and Giuseppe F. Italiano.  Improved bounds and new trade-offs for dynamic all pairs shortest paths.  In ICALP, volume 2380 of Lecture Notes in Computer Science, pages 633–643. Springer, 2002.
  * [DI04] Camil Demetrescu and Giuseppe F. Italiano.  A new approach to dynamic all pairs shortest paths.  J. ACM, 51(6):968–992, 2004.  Announced at STOC‘03.
  * [DI05] Camil Demetrescu and Giuseppe F. Italiano.  Trade-offs for fully dynamic transitive closure on dags: breaking through the o(n${}^{\mbox{2}}$ barrier.  J. ACM, 52(2):147–156, 2005.
  * [DI08] Camil Demetrescu and Giuseppe F. Italiano.  Mantaining dynamic matrices for fully dynamic transitive closure.  Algorithmica, 51(4):387–427, 2008.
  * [DWZ22] Ran Duan, Hongxun Wu, and Renfei Zhou.  Faster matrix multiplication via asymmetric hashing.  CoRR, abs/2210.10173, 2022.
  * [EFGW21] Jacob Evald, Viktor Fredslund-Hansen, Maximilian Probst Gutenberg, and Christian Wulff-Nilsen.  Decremental APSP in unweighted digraphs versus an adaptive adversary.  In ICALP, volume 198 of LIPIcs, pages 64:1–64:20. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.
  * [ES81] Shimon Even and Yossi Shiloach.  An on-line edge-deletion problem.  J. ACM, 28(1):1–4, 1981.
  * [FGH21] Sebastian Forster, Gramoz Goranci, and Monika Henzinger.  Dynamic maintenance of low-stretch probabilistic tree embeddings with applications.  In SODA, pages 1226–1245. SIAM, 2021.
  * [FGNS23] Sebastian Forster, Gramoz Goranci, Yasamin Nazari, and Antonis Skarlatos.  Bootstrapping dynamic distance oracles.  In 31st Annual European Symposium on Algorithms, ESA 2023, September 4-6, 2023, Amsterdam, The Netherlands, volume 274 of LIPIcs, pages 50:1–50:16. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023\.
  * [FMNZ01] Daniele Frigioni, Tobias Miller, Umberto Nanni, and Christos D. Zaroliagis.  An experimental study of dynamic algorithms for transitive closure.  ACM J. Exp. Algorithmics, 6:9, 2001.
  * [GR21] Yong Gu and Hanlin Ren.  Constructing a distance sensitivity oracle in o(n^2.5794 M) time.  In ICALP, volume 198 of LIPIcs, pages 76:1–76:20. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.
  * [GU18] Francois Le Gall and Florent Urrutia.  Improved rectangular matrix multiplication using powers of the coppersmith-winograd tensor.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2018, New Orleans, LA, USA, January 7-10, 2018, pages 1029–1046. SIAM, 2018.
  * [GW20a] Maximilian Probst Gutenberg and Christian Wulff-Nilsen.  Decremental SSSP in weighted digraphs: Faster and against an adaptive adversary.  In SODA, pages 2542–2561. SIAM, 2020.
  * [GW20b] Maximilian Probst Gutenberg and Christian Wulff-Nilsen.  Fully-dynamic all-pairs shortest paths: Improved worst-case time and space bounds.  In SODA, pages 2562–2574. SIAM, 2020.
  * [GWW20] Maximilian Probst Gutenberg, Virginia Vassilevska Williams, and Nicole Wein.  New algorithms and hardness for incremental single-source shortest paths in directed graphs.  In STOC, pages 153–166. ACM, 2020.
  * [HK95] Monika Rauch Henzinger and Valerie King.  Fully dynamic biconnectivity and transitive closure.  In FOCS, pages 664–672. IEEE Computer Society, 1995.
  * [HKN13] Monika Henzinger, Sebastian Krinninger, and Danupon Nanongkai.  Sublinear-time maintenance of breadth-first spanning tree in partially dynamic networks.  In ICALP (2), volume 7966 of Lecture Notes in Computer Science, pages 607–619. Springer, 2013.
  * [HKN14] Monika Henzinger, Sebastian Krinninger, and Danupon Nanongkai.  Sublinear-time decremental algorithms for single-source reachability and shortest paths on directed graphs.  In STOC, pages 674–683. ACM, 2014.
  * [HKN16] Monika Henzinger, Sebastian Krinninger, and Danupon Nanongkai.  Dynamic approximate all-pairs shortest paths: Breaking the O(mn) barrier and derandomization.  SIAM J. Comput., 45(3):947–1006, 2016.  Announced at FOCS‘13.
  * [HKNS15] Monika Henzinger, Sebastian Krinninger, Danupon Nanongkai, and Thatchaphol Saranurak.  Unifying and strengthening hardness for dynamic problems via the online matrix-vector multiplication conjecture.  In STOC, pages 21–30. ACM, 2015.
  * [HP98] Xiaohan Huang and Victor Y. Pan.  Fast rectangular matrix multiplication and applications.  J. Complex., 14(2):257–299, 1998.
  * [Ita86] Giuseppe F. Italiano.  Amortized efficiency of a path retrieval data structure.  Theor. Comput. Sci., 48(3):273–281, 1986.
  * [JX22] Ce Jin and Yinzhan Xu.  Tight dynamic problem lower bounds from generalized BMM and omv.  In STOC ‘22: 54th Annual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 - 24, 2022, pages 1515–1528. ACM, 2022.
  * [Kin99] Valerie King.  Fully dynamic algorithms for maintaining all-pairs shortest paths and transitive closure in digraphs.  In FOCS, pages 81–91. IEEE Computer Society, 1999.
  * [KL19] Adam Karczmarz and Jakub Lacki.  Reliable hubs for partially-dynamic all-pairs shortest paths in directed graphs.  In ESA, volume 144 of LIPIcs, pages 65:1–65:15. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.
  * [KL20] Adam Karczmarz and Jakub Lacki.  Simple label-correcting algorithms for partially dynamic approximate shortest paths in directed graphs.  In 3rd Symposium on Simplicity in Algorithms, SOSA 2020, Salt Lake City, UT, USA, January 6-7, 2020, pages 106–120. SIAM, 2020.
  * [KMG22] Rasmus Kyng, Simon Meierhans, and Maximilian Probst Gutenberg.  Incremental SSSP for sparse digraphs beyond the hopset barrier.  In SODA, pages 3452–3481. SIAM, 2022.
  * [KMS22] Adam Karczmarz, Anish Mukherjee, and Piotr Sankowski.  Subquadratic dynamic path reporting in directed graphs against an adaptive adversary.  In STOC ‘22: 54th Annual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 - 24, 2022, pages 1643–1656. ACM, 2022.
  * [KS02] Valerie King and Garry Sagert.  A fully dynamic algorithm for maintaining the transitive closure.  J. Comput. Syst. Sci., 65(1):150–167, 2002.  Announced at STOC‘99.
  * [KS23] Adam Karczmarz and Piotr Sankowski.  Sensitivity and dynamic distance oracles via generic matrices and Frobenius form.  CoRR, abs/2308.08870, 2023.
  * [Lac11] Jakub Lacki.  Improved deterministic algorithms for decremental transitive closure and strongly connected components.  In SODA, pages 1438–1445. SIAM, 2011.
  * [LO88] Christos Levcopoulos and Mark H. Overmars.  A balanced search tree with _O_ (1) worst-case update time.  Acta Informatica, 26(3):269–277, 1988.
  * [Mao23] Xiao Mao.  Fully-dynamic all-pairs shortest paths: Likely optimal worst-case update time.  CoRR, abs/2306.02662, 2023.
  * [Mun71] J. Ian Munro.  Efficient determination of the transitive closure of a directed graph.  Inf. Process. Lett., 1(2):56–58, 1971.
  * [Rod08] Liam Roditty.  A faster and simpler fully dynamic transitive closure.  ACM Trans. Algorithms, 4(1):6:1–6:16, 2008.
  * [RZ11] Liam Roditty and Uri Zwick.  On dynamic shortest paths problems.  Algorithmica, 61(2):389–401, 2011.
  * [RZ12] Liam Roditty and Uri Zwick.  Dynamic approximate all-pairs shortest paths in undirected graphs.  SIAM J. Comput., 41(3):670–683, 2012.  Announced at FOCS‘04.
  * [San04] Piotr Sankowski.  Dynamic transitive closure via dynamic matrix inverse (extended abstract).  In 45th Symposium on Foundations of Computer Science (FOCS 2004), 17-19 October 2004, Rome, Italy, Proceedings, pages 509–517. IEEE Computer Society, 2004.
  * [San05] Piotr Sankowski.  Subquadratic algorithm for dynamic shortest distances.  In Computing and Combinatorics, 11th Annual International Conference, COCOON 2005, Kunming, China, August 16-29, 2005, Proceedings, volume 3595 of Lecture Notes in Computer Science, pages 461–470. Springer, 2005.
  * [San07] Piotr Sankowski.  Faster dynamic matchings and vertex connectivity.  In SODA, pages 118–126. SIAM, 2007.
  * [SM10] Piotr Sankowski and Marcin Mucha.  Fast dynamic transitive closure with lookahead.  Algorithmica, 56(2):180–197, 2010.
  * [Tar72] Robert Endre Tarjan.  Depth-first search and linear graph algorithms.  SIAM J. Comput., 1(2):146–160, 1972.
  * [Tho04] Mikkel Thorup.  Fully-dynamic all-pairs shortest paths: Faster and allowing negative cycles.  In SWAT, volume 3111 of Lecture Notes in Computer Science, pages 384–396. Springer, 2004.
  * [Tho05] Mikkel Thorup.  Worst-case update times for fully-dynamic all-pairs shortest paths.  In Proceedings of the 37th Annual ACM Symposium on Theory of Computing, Baltimore, MD, USA, May 22-24, 2005, pages 112–119. ACM, 2005.
  * [UY91] Jeffrey D. Ullman and Mihalis Yannakakis.  High-probability parallel transitive-closure algorithms.  SIAM J. Comput., 20(1):100–125, 1991.
  * [WXXZ23] Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, and Renfei Zhou.  New bounds for matrix multiplication: from alpha to omega.  CoRR, abs/2307.07970, 2023.
  * [Zwi02] Uri Zwick.  All pairs shortest paths using bridging sets and rectangular matrix multiplication.  J. ACM, 49(3):289–317, 2002.

##  Appendix A Algebraic data structure (Theorem 3.1)

In this section we show the algebraic subroutine used by our dynamic graph
algorithms. The dynamic algorithm is based on the reduction from dynamic
distance to dynamic matrix inverse from [San05] and it‘s deterministic
equivalent from [BFN22]. The dynamic matrix inverse algorithm will imply the
following dynamic distance result.

See 3.1

Theorem 3.1 is obtained from the following dynamic matrix inverse data
structure. Note that unlike Theorem 3.1, Lemma A.1 measures the complexity in
the number of arithmetic operations. This is because depending on what field
or ring the matrix $\mathbf{M}$ is from, one iteration may take more than
$O(1)$ time. In particular when we later reduce dynamic distance to dynamic
matrix inverse in order to prove Theorem 3.1, one arithmetic operation will
take $\tilde{O}(h^{2})$ time.

######  Lemma A.1.

Let $\mathbf{M}$ be a dynamic $n\times n$ matrix receiving entry updates that
is promised to stay non-singular throughout all updates. Let
$Y\subseteq[n]\times[n]$ be a dynamic set. For any $\mu\in[0,1]$, there exists
a deterministic data structure maintaining the values
$(\mathbf{M}^{-1})_{s,t}$ for all $(s,t)\in Y$ and supporting:

  * •

single-entry updates to $\mathbf{M}$ in
$O(n^{\omega(1,\mu,1)-\mu}+n^{1+\mu}+|Y|)$ operations in worst-case,

  * •

rank $n^{\delta}$ updates to $\mathbf{M}$ for any $\delta\in[0,1]$ in
$O(n^{\omega(1,\delta,1)})$ operations in worst-case,

  * •

single-element insertions to $Y$ in $O(n^{\mu})$ operations,

  * •

removing elements from $Y$ in $O(1)$ time.

###### Proof.

The data structure internally maintains $n\times n$ matrices $\mathbf{N}$ and
$n\times k$ matrices $\mathbf{L},\mathbf{R}$ (with $k\leq n^{\mu}$) such that

$\displaystyle\mathbf{M}^{-1}=\mathbf{N}-\mathbf{L}\mathbf{R}^{\top}.$ (3)

Further, the data structure explicitly maintains the values
$\mathbf{M}^{-1}_{s,t}$ for all $s,t\in Y$.

#### Initialization.

During initialization, the data structure computes $\mathbf{M}^{-1}$ in
$O(n^{\omega})$ operations and stores this inverse as $\mathbf{N}$. The
matrices $\mathbf{L},\mathbf{R}$ are initialized as all-0 matrices. In
particular, since we compute $\mathbf{M}^{-1}$ explicitly, we know
$\mathbf{M}^{-1}_{s,t}$ for all $s,t\in Y$.

#### Update to set $Y$.

When an index pair $(s,t)$ is added to set $Y$, we must compute
$\mathbf{M}^{-1}_{s,t}$. By representation (3) and matrices
$\mathbf{L},\mathbf{R}$ consisting of at most $n^{\mu}$ columns, this takes
$O(n^{\mu})$ operations. When deleting an index pair from $Y$, we do not need
to compute anything so such an update takes $O(1)$ time.

#### Entry update to $\mathbf{M}$.

When matrix $\mathbf{M}$ changes, we must update the representation (3). By
the Sherman-Morrison identity we have for any two vectors $u,v$ that

$\displaystyle(\mathbf{M}+uv^{\top})^{-1}=\mathbf{M}^{-1}-\frac{\mathbf{M}^{-1}u~{}v^{\top}\mathbf{M}^{-1}}{1+v^{\top}\mathbf{M}^{-1}u}.$
(4)

When $\mathbf{M}$ changes in a single entry, then vectors $u$, $v$ have only
one non-zero entry so the products $\mathbf{M}^{-1}u$ and
$v^{\top}\mathbf{M}^{-1}$ consist of reading one row and column of
$\mathbf{M}^{-1}$. By the representation (3) and $\mathbf{R},\mathbf{L}$
having at most $n^{\mu}$ columns, this takes $O(n^{1+\mu})$ operations.

Let $u^{\prime}:=\mathbf{M}^{-1}u/(1+v^{\top}\mathbf{M}^{-1}u)$ and
$v^{\prime}=(v^{\top}\mathbf{M}^{-1})^{\top}$, then we can set
$\mathbf{L}\leftarrow[\mathbf{L}|u^{\prime}]$,
$\mathbf{R}\leftarrow[\mathbf{R}|v^{\prime}]$ (i.e. append these vectors to
matrices $\mathbf{L},\mathbf{R}$) and thus have

$\displaystyle(\mathbf{M}+uv^{\top})^{-1}=$
$\displaystyle~{}\mathbf{M}^{-1}-\frac{\mathbf{M}^{-1}u~{}v^{\top}\mathbf{M}^{-1}}{1+v^{\top}\mathbf{M}^{-1}u}$
$\displaystyle=$ $\displaystyle~{}\mathbf{M}^{-1}-u^{\prime}v^{\prime\top}$
$\displaystyle=$
$\displaystyle~{}\mathbf{N}-\mathbf{L}\mathbf{R}^{\top}-u^{\prime}v^{\prime\top}$
$\displaystyle=$
$\displaystyle~{}\mathbf{N}-[\mathbf{L}|u^{\prime}][\mathbf{R}|v^{\prime}]^{\top}$

so we again satisfy (3).

Note that every $n^{\mu}$ updates, the matrices $\mathbf{L},\mathbf{R}$ have
$n^{\mu}$ columns. To guarantee these matrices never have more columns than
that, we set $\mathbf{N}\leftarrow\mathbf{N}+\mathbf{L}\mathbf{R}^{\top}$
using $O(n^{\omega(1,1,\mu)})$ operations every $n^{\mu}$ updates, which is
$O(n^{\omega(1,1,\mu)-\mu})$ amortized time per update.

At last, note that we must maintain $\mathbf{M}^{-1}_{s,t}$ for all $s,t\in
Y$. This can be done via

$\displaystyle(\mathbf{M}+uv^{\top})^{-1}_{s,t}=$
$\displaystyle~{}\mathbf{M}^{-1}_{s,t}-u^{\prime}_{s}v^{\prime}_{t}.$

Which takes $O(|Y|)$ time per update.

#### Low rank update to $\mathbf{M}$.

Consider a rank $n^{\delta}$ update to $\mathbf{M}$, i.e., we are given two
matrices $\mathbf{U},\mathbf{V}$ of size $n\times n^{\delta}$ and want to
update $\mathbf{M}\leftarrow\mathbf{M}+\mathbf{U}\mathbf{V}^{\top}$. For these
updates we update matrix $\mathbf{N}$ in (3) but not $\mathbf{L}$ or
$\mathbf{R}$. The Woodbury identity states that

$\displaystyle(\mathbf{M}+\mathbf{U}\mathbf{V}^{\top})^{-1}=\mathbf{M}^{-1}-\mathbf{M}^{-1}\mathbf{U}(\mathbf{I}+\mathbf{V}^{\top}\mathbf{M}^{-1}\mathbf{U})^{-1}\mathbf{V}^{\top}\mathbf{M}^{-1}.$
(5)

Here we can compute
$\mathbf{M}^{-1}\mathbf{U}=(\mathbf{N}-\mathbf{L}\mathbf{R}^{\top})\mathbf{U}$
and
$\mathbf{V}^{\top}\mathbf{M}^{-1}=\mathbf{V}^{\top}(\mathbf{N}-\mathbf{L}\mathbf{R}^{\top})$
in $O(n^{\omega(1,1,\delta)})$ operations. In total, computing
$\mathbf{M}^{-1}\mathbf{U}(\mathbf{I}+\mathbf{V}^{\top}\mathbf{M}^{-1}\mathbf{U})^{-1}\mathbf{V}^{\top}\mathbf{M}^{-1}$
takes
$O(n^{\omega(1,1,\delta)}+n^{\omega(1,\delta,\delta)}+n^{\delta\omega})=O(n^{\omega(1,1,\delta)})$
operations. We subtract the result from $\mathbf{N}$ in $O(n^{2})$ operations,
which is subsumed by the previous complexity terms. After this update, we
still satisfy (3) because

$\displaystyle(\mathbf{M}+\mathbf{U}\mathbf{V}^{\top})^{-1}=$
$\displaystyle~{}\mathbf{M}^{-1}-\mathbf{M}^{-1}\mathbf{U}(\mathbf{I}+\mathbf{V}^{\top}\mathbf{M}^{-1}\mathbf{U})^{-1}\mathbf{V}^{\top}\mathbf{M}^{-1}$
$\displaystyle=$
$\displaystyle~{}\underbrace{\mathbf{N}-\mathbf{M}^{-1}\mathbf{U}(\mathbf{I}+\mathbf{V}^{\top}\mathbf{M}^{-1}\mathbf{U})^{-1}\mathbf{V}^{\top}\mathbf{M}^{-1}}_{\text{new
}\mathbf{N}}-\mathbf{L}\mathbf{R}^{\top}.$

#### Worst-case update time

Note that the rank $n^{\delta}$ update does not affect the size of
$\mathbf{L}$ or $\mathbf{R}$, so the reset happens after exactly every
$n^{\mu}$ element updates. So the amortized
$O(n^{\omega(1,1,\mu)-\mu}+n^{1+\mu}+|Y|)$ update time can be made worst-case
via standard technique, see e.g. [BNS19, Appendix B]. ∎

Lemma A.1 can be extended to dynamic matrix inverse for polynomial matrices.
Note that for a polynomial matrix, the inverse might not exist since it‘s a
matrix over a ring instead of a field. The following Corollary A.2 verifies
that all the required inverses are well-defined.

######  Corollary A.2.

Let $\mathbf{A}$ be a dynamic $n\times n$ matrix receiving entry updates,
$h\in[1,n]$, and
$\mathbf{M}:=(\mathbf{I}-X\mathbf{A})\in(\mathbb{F}[X]/\langle
X^{h+1}\rangle)^{n\times n}$. Let $Y\subseteq[n]\times[n]$ be a dynamic set.
For any $\mu\in[0,1]$, there exists a deterministic data structure maintaining
the values $(\mathbf{M}^{-1})_{s,t}$ for all $(s,t)\in Y$ and supporting:

  * •

single-entry updates to $\mathbf{A}$ in
$\widetilde{O}(h\cdot(n^{\omega(1,\mu,1)-\mu}+n^{1+\mu}+|Y|))$ operations in
worst-case,

  * •

rank $n^{\delta}$ updates to $\mathbf{A}$ for any $\delta\in[0,1]$ in
$\widetilde{O}(h\cdot n^{\omega(1,\delta,1)})$ operations in worst-case,

  * •

single-element insertions to $Y$ in $\widetilde{O}(h\cdot n^{\mu})$
operations,

  * •

removing elements from $Y$ in $O(1)$ time.

###### Proof.

The proof follows the common technique of extending dynamic matrix inverse to
polynomial matrices [San04, BNS19, BFN22]. First, note that the inverse of
$\mathbf{M}$ exists and is given by

$(\mathbf{I}-X\mathbf{A})^{-1}=\sum_{k=0}^{h}X^{k}\mathbf{A}^{k}.$

This is because

$(\mathbf{I}-X\mathbf{A})\left(\sum_{k=0}^{h}X^{k}\mathbf{A}^{k}\right)=\sum_{k=0}^{h}X^{k}\mathbf{A}^{k}-\sum_{k=1}^{h+1}X^{k}\mathbf{A}^{k}=\mathbf{I}-X^{h+1}\mathbf{A}^{h+1}=\mathbf{I},$

where the last step uses that the entries of the matrix come from
$\mathbb{F}[X]/\langle X^{h+1}\rangle$, i.e. all arithmetic of the matrix
entries is performed modulo $X^{h+1}$, so $X^{h+1}\mathbf{A}^{h+1}$ is an all-
zero matrix.

We obtain Corollary A.2 by running Lemma A.1 on the matrix
$\mathbf{M}:=(\mathbf{I}-X\mathbf{A})$. Note that the complexity increases by
a factor of $\widetilde{O}(h)$ because the arithmetic operations are now
performed on polynomials of degree at most $h$.

The algorithm from Lemma A.1 works on polynomial matrices, because the
division in (4) and inversion of (5) are well defined. To see this, note that
an update to $\mathbf{M}=(\mathbf{I}-X\mathbf{A})$ is given by an update to
$\mathbf{A}$. So we can assume that vector $u$ in (4) (and likewise matrix
$\mathbf{U}$ in (5)) is a multiple of $X$. Thus, the denominator in (4) is a
polynomial of form $(1+X\cdot P(x))$ (and the to be inverted matrix in (5) is
of form $\mathbf{I}+X\mathbf{B}$, where $\mathbf{B}$ is a polynomial matrix)
whose inverse is given by $\sum_{k=0}^{h}X^{k}(-P(X))^{k}$ (and
$\sum_{k=0}^{h}X^{k}(-\mathbf{B})^{k}$). ∎

We can now prove Theorem 3.1 via the reduction from dynamic distance to
dynamic matrix inverse by [San05].

###### Proof of Theorem 3.1.

We maintain the distances in a dynamic graph via a reduction by Sankowski
[San05] (and its deterministic variant from [BFN22]). Note that

$(\mathbf{I}-X\mathbf{A})^{-1}=\sum_{k=0}^{h}X^{k}\mathbf{A}^{k}\in(\mathbb{F}[X]/\langle
X^{h+1}\rangle)^{n\times n}$

and when $\mathbf{A}$ is an adjacency matrix, we have that
$(\mathbf{A}^{k})_{s,t}$ is the number of (not necessarily simple) paths from
$s$ to $t$ of length $k$. In particular, the smallest $k$ such that
$\mathbf{A}^{k}_{s,t}\neq 0$ is the $st$-distance. So by maintaining the
inverse of $(\mathbf{I}-X\mathbf{A})$ for adjacency matrix $\mathbf{A}$ via
data structure Lemma A.1 we can maintain $h$-bounded distances by just
checking the entries of $(\mathbf{I}-X\mathbf{A})^{-1}$.

The maximum number of paths is bounded by $n^{h}$, so we can pick as field
$\mathbb{F}=\mathbb{Z}_{p}$ for some prime of size $n^{h}\leq p\leq 2n^{h}$.
Here every arithmetic operation in $\mathbb{Z}_{p}$ takes $O(h)$ time in the
Word-RAM model with word size $O(\log n)$. Thus the time complexity increases
by a factor $O(h)$ compared to Corollary A.2. ∎

##  Appendix B Omitted proofs

See 4.1

###### Proof.

Suppose we identify $V$ with $\\{1,\ldots,n\\}$. Let us partition $V$ into
$q=O(n/b)$ blocks $V_{1},\ldots,V_{q}$ of $\leq b$ consecutive integers.
Moreover, let $V^{\prime}$ be a copy of $V$, $V\cap V^{\prime}=\emptyset$,
such that each $v\in V$ has a unique corresponding vertex $v^{\prime}\in
V^{\prime}$. Set $V_{\leq i}:=\bigcup_{j=1}^{i}V_{j}$ for $i\geq 0$. For each
$i=1,\ldots,q$, define:

$G_{i}:=\left(V\cup V^{\prime},E(G)\cup\left\\{uv^{\prime}:uv\in E(G)\text{
and }u\in V_{\leq i}\right\\}\right).$

Note that with graphs $G_{i}$ defined as above, a single-edge update to $G$
leads to at most two single-edge updates in each $G_{i}$. Similarly, if vertex
updates are supported, then a vertex update to $G$ can be reflected in each
$G_{i}$ using $O(1)$ vertex updates therein. Moreover, we have $G_{q}=G$.

Observe that for any $u,v\in V$, $u\neq v$, there is a 1-1 correspondence
between $u\to v$ paths $P=P^{\prime}\cdot xv$ in $G$ with $x\in V_{\leq i}$
and paths $u\to v^{\prime}=P^{\prime}\cdot xv^{\prime}$ in $G_{i}$. As a
result, we have $\delta_{G}^{h}(u,v)=\delta_{G_{i}}^{h}(u,v^{\prime})$ iff
some shortest $h$-bounded path in $G$ has the penultimate vertex in $V_{\leq
i}$. Moreover, every path in $G$ is preserved in $G_{i}[V]$, so
$\delta_{G}^{h}(u,v)=\delta_{G_{i}}^{h}(u,v)$ for all $i$.

We construct and maintain data structures
$\mathcal{D}(G_{1}),\ldots,\mathcal{D}(G_{q})$. Since the $q=O(n/b)$ graphs
$G_{i}$ are at most twice as large as $G$, the initialization time is
$O\left(\frac{n}{b}\cdot I_{\mathrm{source}}(n,h)\right)$ and the update time
is $O\left(\frac{n}{b}\cdot U_{\mathrm{source}}(n,h)\right)$.

Upon query $(S,Y)$, for each $s\in S$ we issue a single-source query to each
$\mathcal{D}(G_{i})$, so that the values $\delta_{G_{i}}^{h}(s,t)$ for all
$(s,t)\in S\times V$ and $i$ are computed. This takes $O((n/b)\cdot|S|\cdot
Q_{\mathrm{source}}(n,h))$ time. Within the same bound, for each $(s,t)\in
S\times V$, we (naively) locate the smallest index $i_{s,t}$ such that
$\delta_{G_{i_{s,t}}}^{h}(s,t^{\prime})=\delta_{G_{q}}^{h}(s,t^{\prime})$.

For each $(s,t)\in Y$, we find the predecessor on some shortest $h$-bounded
$s\to t$ path in $G$ as follows. Find the first vertex $p_{s,t}\in
V_{i_{s,t}}$ such that $p_{s,t}t\in E(G)$ and
$\delta_{G_{q}}^{h}(s,t^{\prime})=\delta_{G_{q}}^{h}(s,p_{s,t})+1$, and report
$p_{s,t}$ as the sought predecessor. Since $|V_{i_{s,t}}|\leq b$, through all
$(s,t)$ this takes $O(|Y|b)$ time.

We now argue that the above is correct. By the definition of $i_{s,t}$, we
have either $i_{s,t}=1$ or
$\delta_{G_{i_{s,t}-1}}^{h}(s,t)>\delta_{G_{q}}^{h}(s,t)$. As a result, (1) no
shortest $h$-bounded $s\to t$ path in $G=G_{q}$ has the penultimate vertex in
$V_{\leq(i_{s,t}-1)}$, (2) at least one such path has the penultimate vertex
in $V_{\leq i_{s,t}}$. But $V_{\leq i_{s,t}}\setminus
V_{\leq(i_{s,t}-1)}=V_{i_{s,t}}$ and thus the penultimate vertex $z$ on each
such path has to satisfy $z\in V_{i_{s,t}}$, $zt\in E(G)$ and
$\delta_{G}^{h}(s,z)+1=\delta_{G}^{h}(s,t)$. On the other hand, any $z$
satisfying this condition is a valid predecessor vertex. By
$\delta_{G_{q}}^{h}(s,t^{\prime})=\delta_{G}^{h}(s,t)$ and
$\delta_{G_{q}}^{h}(s,z)=\delta_{G}^{h}(s,z)$ we obtain that the found vertex
$p_{s,t}$ is a valid predecessor. ∎

See 4.2

###### Proof.

We will only describe the changes needed to make to the construction of Lemma
4.1. We use only one data structure $\mathcal{D}(G)$; for the $q=O(n/b)$
auxiliary graphs $G_{i}$, we use data structures
$\mathcal{D}^{\prime}(G_{1}),\ldots,\mathcal{D}^{\prime}(G_{q})$. The update
and initialization bounds follow easily.

Recall that in the proof of Lemma 4.1, once we have computed, for each
$(s,t)\in Y$, the smallest index $i_{s,t}$ such that
$\delta^{h}_{G_{i_{s,t}}}(s,t^{\prime})=\delta_{G_{q}}^{h}(s,t^{\prime})=\delta_{G}^{h}(s,t)$,
locating the desired predecessors required $O(|Y|b)$ additional time. The
values $\delta_{G}^{h}(s,t)$ for all $(s,t)\in Y$ can be clearly computed by
issuing $|S|$ single-source queries to $\mathcal{D}(G)$, i.e., in $O(|S|\cdot
Q_{\mathrm{source}}(n,h))$ time. By $G_{1}\subseteq
G_{2}\subseteq\ldots\subseteq G_{q}$, we have
$\delta^{h}_{G_{1}}(s,t^{\prime})\geq\delta^{h}_{G_{2}}(s,t^{\prime})\geq\ldots\geq\delta^{h}_{G_{q}}(s,t^{\prime})=\delta^{h}_{G}(s,t)$.
As a result, the desired smallest $i_{s,t}$ such that
$\delta_{G_{i_{s,t}}}^{h}(s,t^{\prime})=\delta_{G}^{h}(s,t)$ can be located
via binary search using $O(\log{n})$ queries about
$\delta_{G_{i}}^{h}(s,t^{\prime})$ issued to $\mathcal{D}^{\prime}({G_{i}})$.
Each of these queries costs $Q_{\mathrm{pair}}(n,h)$ time. ∎

See 4.4

###### Proof.

Let a graph $G_{W}=(V^{\prime},E^{\prime})$ be obtained as follows. For each
$v\in V$, there are two copies $v^{+}$ and $v^{-}$ of $v$ in $V^{\prime}$, and
an auxiliary edge $e_{v}=v^{+}v^{-}\in E^{\prime}$ iff $v\in W$. Moreover,
each edge $uv\in E$ gives rise to the edge $u^{-}v^{+}$ in $E^{\prime}$. This
way, $|V^{\prime}|=2|V|$ and $|E^{\prime}|=|E|+|W|$.

Clearly, an edge update to $G$ can be simulated using a single edge update to
$G_{W}$. Similarly, a vertex update to $G$ can be simulated using $O(1)$
vertex updates to $G_{W}$. Finally, an element update to $W$ translates to a
single auxiliary edge update to $G_{W}$. We set
$\mathcal{D}^{\prime}(G):=\mathcal{D}(G_{W})$.

Observe a 1-1 correspondence of $k$-edge paths $P=s\to t$ in $G[W]$ and
$2k$-edge paths $s^{+}\to t^{+}$ in $G_{W}$. Moreover, no odd-length $s^{+}\to
t^{+}$ path exists in $G_{W}$. As a result, we have
$\delta_{G}^{h}(s,t)=\frac{1}{2}\cdot\delta_{G}^{2h}(s^{+},t^{+})$. The
submatrix predecessor query $(S,Y)$ in $G[W]$ is translated to the query
$\left(\\{s^{+}:s\in S\\},\\{(s^{+},t^{+}):(s,t)\in Y\\}\right)$ issued to
$G_{W}$. By the path correspondence and the fact that all incoming edges of
$t^{+}$ in $G_{W}$ correspond to original incoming edges of $t$ in $G$, we can
indeed answer a submatrix predecessor query on $G[W]$ in $O(Q(|S|,|Y|))$ time.
∎


††thanks: Corresponding Author††thanks: Accepted for publication Phys. Rev. E
102 06210 (2020). DOI: 10.1103/PhysRevE.102.062102

# Susceptibility of quasiclassical Brownian motion in harmonic nonlinear
potentials

Pedro J. Colmenares<EMAIL_ADDRESS>Departamento de Química. Universidad
de Los Andes. Mérida 5101, Venezuela

###### Abstract

This work sets the exact equations for the quasiclassical response function
and susceptibility of a Brownian particle immersed in a bath of quantum
harmonic oscillators driving by nonlinear harmonic potentials. A delta force
perturbation gives rise to a response whose susceptibility is the combination
of a linear term, own of the harmonic oscillator, plus a nonlinear one
involving an integral equation. It is provided a recursion method to find its
solutions based on functional equations in the Banach space. The ODE for the
response function is a highly nonlinear damped non-autonomous Duffing equation
for which the aforementioned method is used to get its solution.

Quantum statistical mechanics, Brownian motion, Stochastic processes, Fokker
Planck equation.

###### pacs:

05.30.Ch; 05.30.?d; 05.40.Jc; 02.50.Ey

##  I Introduction

The response of a system to an external disturbance depends on the
characteristics of the dynamics and of the applied force. Linear response
theory allows one to determine the resonant frequencies of the system with
those of the disturbance. Therefore, it is expected that for nonlinear
dynamics, the resonant frequencies different from those of the external field
will appear as nonlinear terms depending on the magnitude of the force. Other
properties, such as the correlation function of the physical variables Grabert
_et al._ (1984), as well as transport coefficients Colmenares (2019), can be
determined once the susceptibility has been calculated.

The statistical analysis of fluctuations in the Markovian classical Langevin
equation with non-linear harmonic potentials is a well-studied problem. The
Fokker-Planck equation associated with the dynamics of the system is
determined with stochastic techniques through the Ito formula Gardiner (1985).
Once the eigenvalues of the FPE operator are determined as it is done in
Risken’s book Risken (1984), other properties such as the Kramers escapes rate
and diffusion coefficients, to cite a few, are estimated. This has been done
for the normal and metastable (inverted) double-well potential and the
asymmetric one, too. A review by Landa et al Landa and MacClintock (2000)
provides a pedagogical presentation of Kramer’s formula for the mean first
passage time of the Brownian particle leaving any of the bistable wells as a
function of its initial position. None of these works attack the calculation
of the susceptibility induced by nonlinear potentials. In view of this, it
will amenable to extend these well known works to determine their associated
nonlinear susceptibilities. It is important that the theory to be developed
should reduce to the known non trivial result of the harmonic oscillator (HO)
already estimated by Grabert et al Grabert  _et al._ (1984).

There is an extensive literature about the theory of nonlinear response. In
general, they are suited to the problem at hand. For instance, Diezemann
Diezemann (2018) searched the susceptibility in the Markov approximation of
the dielectric relaxation in the reorientational motion of molecules by
expanding the probability density function (PDF) in terms of the external
field. For quantum systems, Grond et. al. Grond  _et al._ (2013) applied a
general linear-response many-body theory to trapped interacting bosonic
systems and Kryvohuz et. al. for chemical reactions Kryvohuz and Mukamel
(2014). Corberi et al. considered the two-time correlation functions by
studying their relation with the second-order susceptibility in ferromagnetic
substances Corberi  _et al._ (2010) while Bouchaud et al. made it for glassy
ones Bouchaud and Biroli (2005). Marini et al. Marconi  _et al._ (2008)
addressed an extensive review of response theory in stochastic and Fokker-
Planck models (see the references there in).

The objective of this work is to formally characterize with stochastic
calculus and standard statistical mechanics tools, the solution, response
function and susceptibility associated to the quasiclassical Langevin equation
(QCLE) Colmenares (2018) with a general nonlinear harmonic potential bias in
the configuration space. The equilibrium system is disturbed with an external
delta force such that a linear response term accompanies a nonlinear one. The
nonlinearity shows up as an integral term depending on the property to be
determined and the dissipation throughout the quantum noise correlation
function. It is in some respect, a partial complement to the analysis gained
in the dynamical description of the damped HO by Graber et al. Grabert  _et
al._ (1984) to the reign of nonlinear harmonic oscillators. The generalized
potential used in this article predicts the results of the HO as a particular
case. Therefore, any extra terms are due to the non-linearity of the
potential.

The paper is structured as follows. The relevant equations are presented in
Sec. II divided in three subsections. Section II.1 covers the system dynamical
equation and its solution in particular, the probability distribution function
(pdf) and associated moments and standard deviation of a Brownian particle
submerged in a quantum harmonic oscillator bath under the effect of a general
external nonlinear potential. Subsequently in Sec. II.2, is shown the
derivation of the response function differential equation and susceptibility
using the solution of the QCLE. Finally, Sec. III shows that the response ODE
is in fact a special kind of the Duffing equation. It is provided its solution
with the method used in Sec II.2. The paper ends in Sec. IV with some final
remarks.

##  II Relevant equations

###  II.1 Solution of the QCLE

The dynamics of a Brownian particle immersed in a reservoir of quantum
harmonic oscillators at the temperature $T$ with a frequency dependent damping
spectra in a position-time dependent potential is of a great interest for its
application in quantum systems. This has been done through the path integral
approach. Two approaches have been developed: either Caldeira et al.’s
Caldeira and Leggett (1983) or Grabert et al.’s Grabert  _et al._ (1988).
Their difference lies in the way the initial total density matrix is handled.
In the former, the total density is factorized such that there is not
correlation between the initial states of the particle and the bath while in
the latter is incorporated. This leads that Caldeira’s equation of motion is
not homogeneous leaving an extra term due to the preparation procedure.
Grabert formalism is chosen in this research. In it, is supposed that the
field only affects the particle Hamiltonian leaving unaltered the bath and
that in addition, the time dependent friction term obeys Kubo’s second
fluctuation-dissipation (SFD) theorem. It is an approximated theory in view of
the molecular dynamics results of Daldrop et al Daldrop  _et al._ (2017). In
it, a solute-water molecules system in a harmonic oscillator field shows that
the friction coefficient depends on the strength of the field. This physical
fact requires one to modify the bath Hamiltonian to include it. Recently, Lisý
and Tóthova Lisý and Tóthová (2019) did that and showed that in the framework
of the generalized Langevin equation (GLE), this extra term in the Hamiltonian
modifies the SFD and therefore, the time dependent friction will depend upon
this strength. A similar result for a time dependent field was previously
found by Olivares and Colmenares Olivares and Colmenares (2016) . The
incorporation of the field strength will require the reformulation of the
functional integral approach that is out of scope in this research.

Using the premises of Grabert et al Grabert  _et al._ (1988), the resulting
equation of motion is the quasiclassical GLE with a noise correlation function
provided by the theory itself and defined later; that is

$\ddot{q}(t)=-\int_{0}^{t}dy\,\Gamma(t-y)\dot{q}(y)-\frac{\partial(V(q,t)}{\partial
t}+\xi(t),$ (1)

where $q$ is the position, $V(q,t)$ the external potential and the friction
kernel $\Gamma(t)$ is the time dependent damping.

Assuming a frequency independent friction $\Gamma(t)=2\gamma\,\delta(t)$ with
$\gamma$ being the static friction coefficient and a time independent general
nonlinear potential $V(q)=\eta\,q^{2}/2+\alpha\,q^{4}/4+\epsilon\,q$,
encompassing parabolic ($\eta=1;\,\alpha=\epsilon=0$), bistable
($\eta\,=\,1;\,\,\alpha\,>\,0;\,\,\epsilon=0$) and asymmetric bistable
($\eta\,=-1\,\,\alpha\,>\,0;\,\,\epsilon>0$), the resulting dynamics is
controlled by the following Ohmic quasiclassical Langevin equation (QCLE).

$\ddot{q}(t)=-\gamma\,\dot{q}(t)-\eta\,q(t)-\alpha\,q(t)^{3}-\epsilon+\xi(t).\\\
$ (2)

The quantum noise $\xi(t)$ depends on the initial position of the Brownian
particle and the initial momentum and position of the bath oscillators Ingold
(2002). It is Gaussian with zero mean and with a two–time correlation function
given by

$\displaystyle\left<\xi(t)\,\xi(s)\right>=$ $\displaystyle-$
$\displaystyle\left(\frac{\gamma\,T}{2}\right)\nu\sinh^{-2}\left[\frac{1}{2}\nu(t-s)\right]$
(3) $\displaystyle+$
$\displaystyle\,i\frac{2\,\pi\,\gamma}{\nu}\dot{\delta}(t-s),$

where $\nu$ is the reduced Matsubara frequency Colmenares (2018).

Unlike the classical Markovian Langevin equation, the noise has a colored
spectrum and is also correlated with the initial position
$q_{{}_{\mathbf{0}}}$ because of the preparation procedure. The latter is
Schramm  _et al._ (1985):

$\displaystyle\left<\xi(t)q_{{}_{\mathbf{0}}}\right>=-2\,\gamma\,T\sum_{n=1}^{\infty}\frac{\nu_{{}_{n}}}{(\nu_{{}_{n}}\\!+\\!\lambda_{{}_{\mathbf{1}}})(\nu_{{}_{n}}\\!+\\!\lambda_{{}_{\mathbf{2}}})}\mbox{e}^{-\nu_{{}_{n}}t},$
(4)

where $\nu_{{}_{n}}=n\nu$ and frequencies
$\lambda_{{}_{\mathbf{1,2}}}=(\gamma\pm(\gamma^{2}-4)^{1/2}\,)/2$.

Equation (2) is obtained from the physical QCLE by scaling energy, position
and time by the parameters $\kappa\,q_{m}^{2}$, $q_{m}$ and
$(M/\kappa)^{1/2}$, respectively, where $M$ is the mass of the particle,
$\kappa$ is the stiffness of the harmonic term of the potential and $q_{m}$ is
chosen arbitrarily. For the bistable and asymmetric bistable potentials could
be the minimum located at the right, while for the HO is of free choice except
$q_{m}=0$. $T=k_{{}_{\mathbf{B}}}\widetilde{T}(\kappa\,q_{m}^{2})^{-1}$ where
$\widetilde{T}$ is the actual temperature in kelvin and $k_{{}_{\mathbf{B}}}$
is the Boltzmann constant.

The QCLE differs from the Markovian classical LE in which the noise
correlation is still Gaussian but delta correlated. The noise average of the
HO QCLE is indistinguishable from the classical. Thus, any calculations made
on it will give the classical result. The QCLE extends to the nonlinear reign
of the HO Markovian classical Langevin equation.

The statistical characterization of the dynamics driven by Eq. (2) is based on
the adaptation of the method developed in Colmenares (2018) for the HO to find
the PDF and its standard deviation.

The Laplace transformation of the QCLE is:

$\displaystyle\widehat{q}_{{}_{\mathrm{L}}}(s)$ $\displaystyle=$
$\displaystyle\widehat{\chi}_{{}_{\mathbf{q}}}\\!(s)\,q_{{}_{\mathbf{0}}}+\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s)\,v_{{}_{\mathbf{0}}}-\frac{\epsilon}{s}\,\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s)$
(5) $\displaystyle+$
$\displaystyle\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s)\bigg{(}\widehat{\xi}(s)-\alpha\,\widehat{H}_{{}_{\mathrm{L}}}(s)\bigg{)},$

where $q_{{}_{\mathbf{0}}}$ and $v_{{}_{\mathbf{0}}}$ are the initial position
and velocity,

$\displaystyle\widehat{\chi}_{{}_{\mathbf{q}}}\\!(s)$ $\displaystyle=$
$\displaystyle(s+\gamma)\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s),$ (6)
$\displaystyle\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s)$ $\displaystyle=$
$\displaystyle\frac{1}{s^{2}+\gamma\,s+\eta},$ (7)

and function
$\widehat{H}_{{}_{\mathrm{L}}}(s)=\mathcal{L}\big{\\{}q^{3}(t)\big{\\}}$ obeys
the identity

$\displaystyle\mathcal{L}\big{\\{}q^{3}(t)\big{\\}}$ $\displaystyle=$
$\displaystyle\Bigg{(}\mathcal{L}\bigg{\\{}\big{(}q^{3}\star q^{3}\star
q^{3}\big{)}(t)\bigg{\\}}\Bigg{)}^{1/3},$ (8)

with the triple convolution given by

$\big{(}q^{3}\star q^{3}\star
q^{3}\big{)}(t)=\int_{0}^{t}\\!\\!dy\,q^{3}(t-y)\int_{0}^{y}\\!\\!dx\,q^{3}(y-x)\,q^{3}(y).$
(9)

After applying the inverse transformation it is found that

$\displaystyle q(t)$ $\displaystyle=$
$\displaystyle\bar{q}(t)+\varphi_{{}_{\\!\mathbf{q}}}(t).$ (10) $\displaystyle
v(t)$ $\displaystyle=$
$\displaystyle\bar{v}(t)+\varphi_{{}_{\\!\mathbf{v}}}(t),$ (11)
$\displaystyle\bar{q}(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi_{{}_{\mathbf{v}}}\\!(t)\,v_{{}_{\mathbf{0}}}+\chi_{{}_{\mathbf{q}}}\\!(t)\,q_{{}_{\mathbf{0}}}}$
(12)
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{t}dy\bigg{[}\epsilon\chi_{{}_{\mathbf{v}}}\\!(y)-\alpha\chi_{{}_{\mathbf{v}}}\\!(t-y)\,H(y)\bigg{]}},$
$\displaystyle\bar{v}(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\dot{\chi}_{{}_{\mathbf{v}}}\\!(t)\,v_{{}_{\mathbf{0}}}+\dot{\chi}_{{}_{\mathbf{q}}}\\!(t)\,q_{{}_{\mathbf{0}}}+\epsilon\,\chi_{{}_{\mathbf{v}}}\\!(t)}$
(13) $\displaystyle-$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{0}^{t}dy\,\alpha\,\dot{\chi}_{{}_{\mathbf{v}}}\\!(t-y)\,H(y)},$

where the noises $\varphi_{{}_{\\!\mathbf{q}}}(t)$ and
$\varphi_{{}_{\\!\mathbf{v}}}(t)$ are defined by

$\displaystyle\varphi_{{}_{\\!\mathbf{q}}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dy\,\chi_{{}_{\mathbf{v}}}\\!(t-y)\,\xi(y),$ (14)
$\displaystyle\varphi_{{}_{\\!\mathbf{v}}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dy\,\dot{\chi}_{{}_{\mathbf{v}}}\\!(t-y)\,\xi(y).$
(15)

The over bar on $q(t)$ and $v(t)$ denotes their average values and the
susceptibilities $\chi_{{}_{\mathbf{q}}}\\!(t)$ and
$\chi_{{}_{\mathbf{v}}}\\!(t)$ read as

$\displaystyle\chi_{{}_{\mathbf{q}}}\\!(t)$ $\displaystyle=$
$\displaystyle\mbox{e}^{-\gamma\,t/2}\Bigg{[}\cosh\bigg{[}\frac{\omega_{0}}{2}t\bigg{]}+\frac{\gamma}{\omega_{0}}\sinh\bigg{[}\frac{\omega_{0}}{2}t\bigg{]}\Bigg{]},$
(16) $\displaystyle\chi_{{}_{\mathbf{v}}}\\!(t)$ $\displaystyle=$
$\displaystyle\frac{2}{\omega_{0}}\mbox{e}^{-\gamma\,t/2}\sinh\bigg{[}\frac{\omega_{0}}{2}t\bigg{]},$
(17)

with $\omega_{0}=\sqrt{\gamma^{2}-4\,\eta}$. The presence of the nonlinear
term in the potential shows up through $\widehat{H}_{{}_{\mathrm{L}}}(s)$.
Therefore, the integral equation in the Laplace space given by Eq. (5) has to
be first numerically solved in order to get all dynamical properties related
to $q(t)$. This partial result is not a simple numerical problem of computing.
Apart from that, it is own of the mathematical formulation.

Since the equations for $q(t)$ and $v(t)$ have the same structure as those
used in Ref. Colmenares (2018), the conditional distribution
$p(q,t|q_{{}_{\mathbf{0}}})$ is Gaussian with mean $G(t)$ and standard
deviation $\sigma^{2}(t)$ given by

$\displaystyle G(t)$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi_{{}_{\mathbf{q}}}\\!(t)\,q_{{}_{\mathbf{0}}}\\!\\!+\\!\\!\int_{0}^{t}\\!\\!\\!dy\bigg{[}\epsilon\chi_{{}_{\mathbf{v}}}\\!(y)-\alpha\chi_{{}_{\mathbf{v}}}\\!(t\\!-\\!y)\,H(y)\bigg{]}},$
(18)
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}(t)}$
$\displaystyle=$ $\displaystyle
T\,\chi_{{}_{\mathbf{v}}}\\!^{2}(t)+2\int_{0}^{t}dy\Bigg{[}\int_{0}^{y}dz\,\big{<}\varphi_{{}_{\\!\mathbf{v}}}(z)\varphi_{{}_{\\!\mathbf{v}}}(y)\big{>}$
(19) $\displaystyle+$
$\displaystyle\chi_{{}_{\mathbf{q}}}\\!(y)\,\big{<}\varphi_{{}_{\\!\mathbf{v}}}(y)q_{{}_{\mathbf{0}}}\big{>}\Bigg{]},$

respectively. The $\eta$ dependence on $G(t)$ is controlled by the
susceptibilities $\chi_{{}_{\mathbf{q}}}\\!(t)$ and
$\chi_{{}_{\mathbf{v}}}\\!(t)$. Similarly, the amplitude of the non-linear
term affects the standard deviation thorough
$\chi_{{}_{\mathbf{v}}}\\!^{2}(t)$.

It is worth mentioning that the preceding procedure can be applied to the
generalized QCLE. Important to keep in mind is that
$\widehat{\chi}_{{}_{\mathbf{v}}}\\!(s)$, Eq. (7), will not depend on $\gamma$
but on the Laplace transform of the kernel, that is, $\widehat{\Gamma}(s)$. In
order to get a closed equation, is imperative to define first the functional
form of the time-dependent kernel Olivares and Colmenares (2016).

###  II.2 Response Function ODE and Susceptibility

According to Ehrenfest’s theorem, Eq. (2) reads:

$\displaystyle F(t)$ $\displaystyle=$
$\displaystyle\Big{<}\ddot{q}(t)\Big{>}+\gamma\,\Big{<}\dot{q}(t)\Big{>}+\Big{(}3\,\alpha\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}(t)}+\eta\Big{)}\,\Big{<}q(t)\Big{>}$
(20) $\displaystyle+$ $\displaystyle\alpha\,\Big{<}q(t)\Big{>}^{3}+\epsilon,$

where the external driving force $F(t)$ was added . The third moment of $q(t)$
was written as
$\big{<}q^{3}(t)\big{>}=\big{(}\big{<}q(t)\big{>}^{3}+3\,\big{<}q(t)\big{>}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}(t)}\big{)}$
with $\sigma^{2}(t)$ given by Eq. (19).

The response function $R(t)$ of the system as the result of the perturbation
driving by $F(t)$ is defined by
$\big{<}q(t)\big{>}=\int_{-\infty}^{\infty}dy\,R(t-y)\,F(y)$. Then for a
driving force $F(t)=F_{0}\,\delta(t)$, we obtain the following ODE for $R(t)$

$\displaystyle\delta(t)$ $\displaystyle=$
$\displaystyle\ddot{R}(t)+\gamma\,\dot{R}(t)+\Big{(}3\,\alpha\,\sigma^{2}(t)+\eta\Big{)}R(t)$
(21) $\displaystyle+$
$\displaystyle\alpha\,F_{0}^{2}\,R^{3}(t)+\frac{\epsilon}{F_{0}},$

with initial conditions $R(0)=0$ and $\dot{R}(0)=1.$

This equation for finite values of $\alpha$ and $\epsilon$ gives rise to a non
linear response due to the dependence on $F_{0}$. For the parabolic potential,
the perturbation excites a response of the same frequency Reichl (1998).

The susceptibility $\chi(\omega)$ is defined as the Fourier transform (FT) of
the response function, that is,
$\chi(\omega)=\int_{0}^{\infty}dt\,R(t)\mbox{e}^{\imath\,\omega\,t}$.
Recalling the convolution identities

$\displaystyle\mathcal{F}\Big{\\{}R(t)\,\sigma(t)\Big{\\}}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dy\,\chi(\omega-y)\,\widehat{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}}_{{}_{\\!\mathrm{F}}}(y),$
(22) $\displaystyle\mathcal{F}\Big{\\{}R^{3}(t)\Big{\\}}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dy\,\chi(\omega-y)$ (23)
$\displaystyle\times$
$\displaystyle\int_{-\infty}^{\infty}dz\,\chi(z)\,\chi(y-z),$

then, the FT of the response function’s ODE, Eq. (21), gives the
susceptibility,

$\displaystyle\chi(\omega)$ $\displaystyle=$
$\displaystyle=\phi(\omega)+\psi\big{(}\omega\big{)},$ (24)
$\displaystyle\phi(\omega)$ $\displaystyle=$
$\displaystyle\widetilde{\chi}(\omega)\bigg{(}1{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}\frac{\epsilon}{F_{0}}\delta(\omega)\Bigg{)},$
(25) $\displaystyle\psi\big{(}\omega,[\chi]\big{)}$ $\displaystyle=$
$\displaystyle\widetilde{\chi}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\,\chi(\omega-\omega^{\prime})\Bigg{[}3\,\widehat{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}}_{{}_{\\!\mathrm{F}}}(\omega^{\prime})$
(26)
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\,F_{0}^{2}\int_{-\infty}^{\infty}d\omega^{\prime\prime}\,\chi(\omega^{\prime\prime})\,\chi(\omega^{\prime}-\omega^{\prime\prime})\Bigg{]},$

where
$\widetilde{\chi}(\omega)=(\eta-\omega^{2}-\imath\,\gamma\,\omega)^{-1}$. It
is a rather complicated equation for non-linear potentials that involves
solving a Fredholm integral equation of the second kind. However, for the HO
($\eta=1\,;\,\alpha=\epsilon=0$) $\chi(\omega)=\widetilde{\chi}(\omega)$
agrees with that of the reported in the literature Grabert  _et al._ (1984);
Reichl (1998).

The proposed solution of this integral equation uses a method developed by
Daftardar-Gejji et al. Daftardar-Gejji and Jafari (2006). Its application to a
variety of integral equations including ODEs, are of a high order of accuracy
Daftardar-Gejji and Jafari (2006); Al-Jawary and Abd-Al-Razaq (2016). It has
been chosen because of its simplicity. In order to show the scope of the
algorithm, it will be fully described next.

Considering $\psi\big{(}\omega,[\cdot]\big{)}$ as the Banach operator
$\mathcal{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B}}$
in the Banach space $\mathcal{B}\rightarrow\mathcal{B}$ and
$f(\omega)=\phi(\omega)$, then Eq. (24) can be written as

$\chi(\omega)=f(\omega)+\mathcal{B}\big{(}\chi(\omega)\big{)}.$ (27)

Let the solution be given by

$\chi(\omega)=\sum_{i=0}^{\infty}\chi_{i}(\omega),$ (28)

and decomposw the operator $\mathcal{B}\big{(}\chi(\omega)\big{)}$ according
to

$\displaystyle\mathcal{B}\Bigg{(}\sum_{i=0}^{\infty}\chi_{i}(\omega)\Bigg{)}$
$\displaystyle=$
$\displaystyle\mathcal{B}\big{(}\chi_{{}_{0}}(\omega)\big{)}+\sum_{i=1}^{\infty}\Bigg{[}\mathcal{B}\Bigg{(}\sum_{j=0}^{i}\chi_{j}(\omega)\Bigg{)}$
(29) $\displaystyle-$
$\displaystyle\mathcal{B}\Bigg{(}\sum_{j=0}^{i-1}\chi_{j}(\omega)\Bigg{)}\Bigg{]},$

which is equivalent to

$\displaystyle\sum_{i=0}^{\infty}\chi_{i}(\omega)$ $\displaystyle=$
$\displaystyle f(\omega)+\mathcal{B}\big{(}\chi_{{}_{0}}(\omega)\big{)}$ (30)
$\displaystyle+$
$\displaystyle\sum_{i=1}^{\infty}\\!\Bigg{[}\mathcal{B}\Bigg{(}\\!\sum_{j=0}^{i}\chi_{j}(\omega)\Bigg{)}\\!-\\!\mathcal{B}\Bigg{(}\\!\sum_{j=0}^{i-1}\chi_{j}(\omega)\Bigg{)}\\!\Bigg{]}.$
(31)

Defining the following recursion relation:

$\displaystyle\chi_{{}_{0}}(\omega)$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=}$
$\displaystyle f(\omega),$ (32) $\displaystyle\mathcal{\chi_{{}_{1}}}(\omega)$
$\displaystyle=$ $\displaystyle\mathcal{B}\bigg{(}\chi_{{}_{0}}(t)\bigg{)},$
(33) $\displaystyle\chi_{m+1}(\omega)$ $\displaystyle=$
$\displaystyle\mathcal{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B}}\bigg{(}\chi_{{}_{0}}(\omega)+\cdots+\chi_{m}(\omega)\bigg{)}$
(34) $\displaystyle-$
$\displaystyle\mathcal{\mathcal{B}}\bigg{(}\\!\\!\chi_{{}_{0}}(\omega)+\cdots+\chi_{m-1}(\omega)\\!\\!\bigg{)};\,\,m\geqslant
1,$

then

$\chi(t)=f(\omega)+\sum_{i=1}^{\infty}\chi_{i}(\omega).$ (35)

Thus the solution is recursively calculated through each of the
$\chi_{i}(\omega)$. The convergence of the summation to be one of the set of
possible solutions is guaranteed in view of the Banach fixed-point theorem
Daftardar-Gejji and Jafari (2006). Since $\chi(\omega)$ is complex, Eq. (35)
can be written as a system of equations for its real and imaginary parts. In
particular, the position correlation function is written as an integral of its
imaginary part with poles located in the lower half plane with values given by
the roots of the denominator of the imaginary part of
$\widetilde{\chi}(\omega)$ as in the HO Grabert  _et al._ (1984).

Apart from $\chi_{{}_{0}}(\omega)$, functions of superior order requires one
to know $\widehat{\sigma^{2}}_{{}_{\\!\mathrm{F}}}(\omega^{\prime})$ which
involves the FT of the multidimensional integral given by Eq. (19). Using the
method described in Ref. [2] allows one to write down the noise correlation in
terms of the numerical position correlation function. Once the standard
deviation and its FT are calculated then it can be fitted as a given function
and subsequently used in the integral depicted in Eq. (26). In general, the
different terms will have among others, large polynomials in their
constitutive equations because the zero order term will appear in each of the
$\chi_{i}(\omega)$. As a collateral result, once the noise correlation is
calculated then the diffusion constant of the system can also be determined as
is exemplified in Colmenares (2019) for the parabolic potential.

Having known $\chi(\omega)$ then $R(t)$ is given by its inverse, that is,

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}R(t)=\frac{1}{2\,\pi}\,\int_{-\infty}^{\infty}d\omega\,\chi(\omega)\,\mbox{e}^{-\imath\,\omega\,t}.}$
(36)

This is the usual method found in textbooks to get the response function. In
the next section the application of the Banach operator method to solve the
response function ODE will be shown.

##  III Response function ODE and Duffing equation

Equation (21) is a highly non linear damped non autonomus Duffing equation
(DE) for which the literature about solution methods and application for
specific systems is very rich. A good review of the application to control
systems is found in Chap. 5 of Kovacic and Brennan’s book Kovavic and Brennan
(2011) and approximated methods of solutions in Nayfeh’s Nayfeh (1973). The
methods are primarily based on perturbation techniques of some small parameter
Nayfeh (1973). Although they have been formally deducted, its algebra is
cumbersome and lengthy. Instead, the operator method discussed above will be
apply.

The conversion of the ODE into the functional form such as given by Eq. (27)
is simple. To do that, let $L_{tt}(\cdot)=\partial^{2}(\cdot)/\partial t^{2}$,
$L_{t}(\cdot)=\partial(\cdot)/\partial t$,
$L_{tt}^{-1}\,(\cdot)=\int_{0}^{t}dy\int_{0}^{y}ds(\cdot)$ and
$L_{t}^{-1}(\cdot)=\int_{0}^{t}dy\,(\cdot)$. Then, applying the operator
$L_{tt}^{-1}(\cdot)$ to the response function’s ODE, using the initial
conditions and that $L_{tt}^{-1}(\cdot)=\int_{0}^{t}dy\,(t-y)(\cdot)$, the
resulting functional equation is of the type of Eq. (27, with

$\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}4-\frac{\epsilon}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}F_{0}}}t\bigg{)}t,$
(37)
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{B}\big{(}R(t)\big{)}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\int_{0}^{t}dy\,\Bigg{\\{}\gamma\,R(y)+(t-y)\bigg{[}\eta\,R(y)}$
(38) $\displaystyle-$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha\bigg{(}F_{0}^{2}\,R^{3}(y)\\!+\\!3\,R(y)\,\sigma^{2}(y)\bigg{)}\bigg{]}\Bigg{\\}}},$

with Eq. (38) defining the operator $\mathcal{B}$ of the system. The
functional equation is a Volterra integral equation type such that the method
mentioned above can be used to find its solution. Each of the $R_{i}(t)$ are
calculated by the recursion formula shown above. The method of analysis
described in Ref. Kovavic and Brennan (2011) could be used to characterize the
dynamical aspects of the solution due to the nonlinearity.

Unfortunately, as in the solution of the susceptibility equation, we cannot
analytically go any further but numerically, because the Banach operators
requires the standard deviation, which in turn must also be determined
numerically.

##  IV Final Remarks

The calculation of the susceptibility associated to the QCLE with nonlinear
potentials using standard tools of stochastic equations and statistical
mechanics leads to an integral equation for this property. It has such a
complexity that numerical methods have to be used to solve it. However, the
proposed recursive method of Banach operators reduces in a large extension the
calculations bypassing any iterative technique in the solution. In any case,
the incorporation of nonlinear potentials and the tools already mentioned do
not allow getting a basic physical insight about the susceptibility just as
that found for the pure harmonic oscillator. Doing numerical calculations can
only elucidate it. This drawback opens the doors to find alternative routes to
tackle this problem in order to obtain at least simpler equations. In any
case, the equations are exact without any assumptions made in their
derivations.

The Banach operator method for functional equations is physically sound. It
shows how any function can be expanded up to a specified order so that
additional calculations whose contributions are negligible can be prevented.
The maximal error remainder test Al-Jawary and Abd-Al-Razaq (2016) can be used
to stop the numerical procedure. This is guaranteed by the Banach fixed-point
theorem.The derivation of the formulas was possible because the PDF or Green
function of the dynamics is completely known. In cases where the PDF is not
available, the theory of Diezemann Diezemann (2018) properly adapted to the
problem at hand, would be the route to find the different orders of the
response function.

The problem of classical Markovian systems is easier to handle because the
noise is delta correlated. In particular the standard deviation is analytical
Colmenares (2019).

A previous work Colmenares (2019) mentions the class of systems potentially
can be used with this formalism.

## Acknowledgments

The author is grateful to the referees for helpful suggestions.

## References

  * Grabert  _et al._ (1984) H. Grabert, U. Weiss, and P. Talkner, Z. Phys. B: Cond. Matt. 55, 87 (1984).
  * Colmenares (2019) P. J. Colmenares, Phys. Rev. E 100, 052117 (2019), The term inside the summation in Eq. (10d) must be $(\nu_{n}\exp[-\nu_{n}t])((\nu_{n}^{2}+\lambda_{1}^{2})(\nu_{n}^{2}+\lambda_{2}^{2})))^{-1}$.
  * Gardiner (1985) C. W. Gardiner,  _Handbook of Stochastic for Physics, Chemistry and the Natural Sciences_ , 2nd ed. (Springer-Verlag, Berlin, 1985).
  * Risken (1984) H. Risken,  _The Fokker-Planck Equation. Method of Solutions and Applications_ (Springer-Verlag, Berlin, 1984).
  * Landa and MacClintock (2000) P. S. Landa and P. V. E. MacClintock, Phys. Rep. 323, 1 (2000).
  * Diezemann (2018) G. Diezemann, in  _Nonlinear Dielectric Spectroscopy. Advances in Dielectrics_, edited by R. Richert (Springer, Berlin, 2018) p. 75.
  * Grond  _et al._ (2013) J. Grond, A. I. Streltsov, A. U. J. Lode, K. Sakmann, L. S. Cederbaum, and O. E. Alon, Phys. Rev. A 88, 023606 (2013).
  * Kryvohuz and Mukamel (2014) M. Kryvohuz and S. Mukamel, J. Chem. Phys. 140, 034111 (2014).
  * Corberi  _et al._ (2010) F. Corberi, E. Lippiello, A. Sarracino, and M. Zannetti, J. Stat. Mech. Theory Exp. , P04003 (2010).
  * Bouchaud and Biroli (2005) J-P. Bouchaud and G. Biroli, Phys. Rev. B 72, 064204 (2005).
  * Marconi  _et al._ (2008) U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Phys. Rep. 461, 111 (2008).
  * Colmenares (2018) P. J. Colmenares, Phys. Rev. E 97, 052126 (2018), There is a misprint in Eq. (12), $\xi(t)$ must be $\varphi_{\mathrm{v}}(t)$.
  * Caldeira and Leggett (1983) A. O. Caldeira and J. Leggett, Physica A 121, 587 (1983).
  * Grabert  _et al._ (1988) H. Grabert, P. Schramm, and G. L. Ingold, Phys. Rep. 168, 115 (1988).
  * Daldrop  _et al._ (2017) J. O. Daldrop, B. G. Kowalik, and R. R. Netz, Phys. Rev. X 7, 041065 (2017).
  * Lisý and Tóthová (2019) V. Lisý and J. Tóthová, Results Phys. 12, 1212 (2019).
  * Olivares and Colmenares (2016) W. Olivares and P. J. Colmenares, Physica A 458, 76 (2016).
  * Ingold (2002) G. L. Ingold, “Path integrals and their application to dissipative quantum systems,” in  _Coherent Evolution in Noisy Environments_ , Lecture Notes in Physics, Vol. 611, edited by A. Buchleiter and K. Hornberger (Springer, Berlin, 2002) Chap. 1, pp. 1–53.
  * Schramm  _et al._ (1985) P. Schramm, R. Jung, and H. Grabert, Phys. Lett. A 107, 385 (1985).
  * Reichl (1998) L. E. Reichl,  _A Modern Course in Statistical Physics_ , 2nd ed. (John Wiley & Sons, New York, 1998).
  * Daftardar-Gejji and Jafari (2006) V. Daftardar-Gejji and H. Jafari, J. Math. Anal. Appl. 316, 753 (2006).
  * Al-Jawary and Abd-Al-Razaq (2016) M. A. Al-Jawary and S. G. Abd-Al-Razaq, Int. J. Bas. Appl. Sci. 5, 115 (2016).
  * Kovavic and Brennan (2011) I. Kovavic and M. J. Brennan, “The Duffing Equation,” (John Wiley & Sons, U.S.A., 2011) 1st ed.
  * Nayfeh (1973) A. Nayfeh,  _Perturbation Methods_ , Pure & Applied Mathematics. A Wiley-Interscience Series of Texts, Monographs & Tracts (John Wiley & Sons, New York, 1973).