\section{Introduction}
 Let $H, A, B$ be  smooth classical observables on $\mathbb{R}^{2d}$ in the variables $X=(x,\xi)$.  The Poisson brackets
 is defined as $\{A,B\} =\partial_\xi A\cdot\partial_xB-\partial_xA\cdot\partial_\xi B$. So the classical  time evolution of $A$  determined by the Hamilton equation for $H$ is solution of the equation:
 \begin{eqnarray}}  \def\eea{\end{eqnarray}
 \frac {d}{dt}A(t) &=& \{A(t), H\} \\
 A(0) &=&A. \nonumber
 \eea
 The Weyl quantization $\widehat A$  of $A$  is defined as as the following operator:
  \begin{equation}}   \def\eeq{\end{equation}
  \widehat A f(x) :={\rm Op}^w_\hbar f(x)=(2\pi\hbar)^{-d} \int_{\mathbb{R}^{2d}} A\!\left(\frac{x+y}{2}, \xi\right)
{\rm e}^{i\xi\cdot(x-y)/\hbar} f(y) \,dy\, d\xi
\eeq
for any $f\in{\mathcal S}(\mathbb{R}^d)$.  \\
The quantum time evolution of the quantum observable $\widehat A$ must satisfy the Heisenberg equation
\begin{eqnarray}}  \def\eea{\end{eqnarray}
 \frac {d}{dt}\widehat A(t) &=& \frac{i}{\hbar}[ \widehat A(t),  \widehat H]\\
  \widehat A(0) &=&\widehat A
 \eea
 where $[\hat A, \hat B] = \hat A\hat B -\hat B\hat A$.\\
The Moyal bracket  of the observables $A, H$ is defined such that
\begin{equation}}   \def\eeq{\end{equation}
\frac{i}{\hbar}[ \widehat A,  \widehat H]={\rm Op}^w_\hbar(\{A,H\}_\circledast.
\eeq
Notice that it results from   the Weyl quantization  calculus with a small parameter $\hbar$  that we have
$$
\lim_{\hbar \searrow 0}\{A,H\}_\circledast = \{A,H\}.
$$

There is an equivalent definition  by introducing the Moyal product  $A\circledast B$ (see the next section)
such that
$$
({\rm Op}^w_\hbar A)({\rm Op}^w_\hbar B) = {\rm Op}^w_\hbar(A\circledast B).
$$
Then we have
$$
\{A,B\}_\circledast  = \frac{i}{\hbar}(A\circledast B-B\circledast A).
$$
These definitions make sense for $A, B\in {\mathcal S}(\mathbb{R}^{2d})$ and can be extended to suitable  classes of symbols with moderate growth. To be more explicite we introduce  the classes
${\mathbb{S}}_\delta^\mu$, for $\delta<1$, $\mu\in\mathbb{R}$. $A\in{\mathbb{S}}_\delta^\mu$ iff  $A\in C^\infty(\mathbb{R}^{2d})$ and   for any multiindex $\gamma\in\mathbb{N}^{2d}$
 we have:
$$
\vert\partial^\gamma_X A(X)\vert \leq C_\gamma\langle X\rangle^{\mu+\delta\vert\gamma\vert}
$$
Using Theorem A.1 in \cite{BR}, we can see that $A\circledast H$ is a smooth symbol if $H\in{\mathbb{S}}_\delta^\mu$ and $A\in{\mathbb{S}}^0_0$.
Our aim here is to prove the following result.
\begin{theorem}\label{MainThm} Assume that  $\hbar$  is fixed ($\hbar =1$).
Let be $H\in {\mathbb{S}}_\delta^\mu$ for some $\mu\in\mathbb{R}$ and $\delta<1$.  Assume that  for any $A\in{\mathbb{S}}^0_0$  we have
$ \{A, H\}_\circledast =  \{A, H\}$. Then $H(X)$ is a polynomial in $X=(x,\xi)$ of degree at most 2.
\end{theorem}
\begin{remark}
It is well known that if $H$ is a polynomial of degree at most 2 then $\{A, H\}_\circledast =  \{A, H\}$  for any $A\in{\mathbb{S}}^\nu_0$. I do not know any  reference  for a proof of the converse statement. The proof given here is a direct consequence of basic properties of the Weyl quantization.
\end{remark}
\begin{remark}
The usual proofs of the Groenewold-van Hove Theorem on the phase space $\mathbb{R}^{2d}$  concern  more general quantization procedures  but are restricted to polynomial symbols $A, H$.\\
       A quotation from \cite{vHove} p.66-67:\\
       "On \'etablit ensuite qu'une correspondance biunivoque entre grandeurs classiques et quantiques, ayant le  caract\`ere d'un isomorphisme entre alg\`ebres de Lie, existe entre les grandeurs repr\'esent\'ees par des polyn\^omes de degr\'e 0, 1, 2 en les variables $p_1,\cdots p_N,q_1,\cdots q_n$ mais ne peut \^etre \'etendue sans perdre ses propri\'et\'es essentielles \`a l'ensemble de toutes les grandeurs classiques"   \\
       The Theorem of Groenewold-van Hove is detailed  p.76 and the quadratic case p.87 of \cite{vHove}.\\
       Notice that the quadratic case is related with the metaplectic representation \cite{Gotay}.
\end{remark}

 \section{Weyl calculus}
 \subsection{Introduction to the Weyl quantization}
 In this section, we recall some  basic properties of the Weyl calculus (for more details see \cite{ho}).\\
 Weyl quantization start by quantization of exponent of linear forms $L_Y(X) = \sigma(Y,X)=\eta\cdot x-y\cdot\xi$
  with $X=(x,\xi)$, $Y=(y,\eta)$. Apart the usual properties asked for an admissible quantization, Weyl quantization is uniquely determined by imposing
  that the Weyl symbol of ${\rm e}^{i\widehat{L_Y}}$ is ${\rm e}^{iL_Y}$. Recall that $ \widehat T(Y):={\rm e}^{-\frac{i}{\hbar}\widehat{L_Y}}$   is the Weyl-Heisenberg translation operator by $Y$ in the phase space $\mathbb{R}^{2d}$.
  Then  for any observable $A$,  using a Fourier transform,  the  Weyl quantization $ A$ is defined  for any $\psi\in {\mathcal S}(\mathbb{R}^d)$, as
  $$
   \widehat A\psi = (2\pi)^{-d}\int_{\mathbb{R}^{2d}}{\tilde A}_\sigma(Y) T(Y)\psi dY
  $$
  where ${\tilde A}_\sigma(Y)=\int_{\mathbb{R}^{2d}} A(z){\rm e}^{-i\sigma(Y,z)}dz$ is the symplectic Fourier transform of $A$ (in the sense of distributions).
   So that the family $\{T(Y)\}_{Y\in\mathbb{R}^{2d}}$ is an over-complete   basis for operators   between the Schwartz spaces ${\mathcal S}(\mathbb{R}^d)$ and
   ${\mathcal S}'(\mathbb{R}^d)$

 \subsection{The Moyal Product}
We first recall the formal product rule for quantum observables
with Weyl
quantization. Let $A, B \in {\mathcal S}(\mathbb{R}^{2d})$. The Moyal product $C:=A\circledast B$ is the
observable
$C$ such that $\widehat{A}\cdot\widehat{B} = \widehat{C}$. Some computations with
the Fourier
transform give the following well known formula \cite{ho}
\begin{equation}}   \def\eeq{\end{equation}\label{moy1}
C(x,\xi) =
\exp\left(\frac{i\hbar}{2}\sigma(D_q,D_p;D_{q^\prime},D_{p^\prime})\right)A(q,p)
B(q^\prime,p^\prime)\vert_{(q,p)=(q^\prime,p^\prime)=(x,\xi)},
\eeq
where $\sigma$ is the symplectic bilinear form $\sigma((q,p), (q',p'))=p\cdot q'-p'\cdot q$ and
$D=i^{-1}\hbar\nabla$.
By expanding the exponential term in a formal power series in $\hbar$ we get
\begin{equation}}   \def\eeq{\end{equation}\label{prod2}
C(x,\xi) =
\sum_{j\geq
0}\frac{\hbar^j}{j!}\left(\frac{i}{2}\sigma(D_q,D_p;D_{q^\prime},D_{p^\prime}) \right)^j
A(q,p)B(q^\prime,p^\prime)\vert_{(q,p)=(q^\prime,p^\prime)=(x,\xi)}.
\eeq
So that $C(x,\xi)$ is a formal power series in $\hbar$ with coefficients
given by
\begin{equation}}   \def\eeq{\end{equation} \label{exprod}
C_j(A,B;x,\xi) = \frac{1}{2^j}\sum_{\vert\alpha+\beta\vert=j}
\frac{(-1)^{\vert \beta\vert}}{\alpha!\beta!}
(D^\beta_x\partial^\alpha_\xi A).( D^\alpha_x\partial^\beta_\xi B)(x,\xi).
\eeq
\medskip
Furthermore we need  a remainder estimates for the expansion of the Moyal product.\\
For every $N \geq 1$. we denote
\begin{equation}}   \def\eeq{\end{equation}
R_N(A,B;z):=
A\circledast B(z) - \sum_{0\leq j \leq N}\hbar^jC_j(z) .
\eeq
The following estimate is a particular case of Theorem A.1 in \cite{BR} see also Remark A.3.

\begin{lem}\label{lem:Moyalest}
Let be $A\in{\mathbb S}^{\mu_A}_\delta$ and  $B\in{\mathbb S}^{\mu_B}_\delta$  then for any $N\geq 1$, $\gamma\in\mathbb{N}^{2d}$, $M\geq M_0$  there exists $C_{N,\gamma,M}>0$(independent of $(A, B)$)
such that
\begin{eqnarray}}  \def\eea{\end{eqnarray}\label{remest1}
\vert\partial_z^\gamma R_N(A,B;z)\vert \leq C_{N,\gamma,M}\hbar^{N+1}\sum_{\substack{\vert\alpha+\beta\vert=N+1\\ \vert\mu+\nu\vert\leq M+\vert\gamma\vert}}\\\sup_{u, v\in\mathbb{R}^{2d}} (1+\vert u\vert^2+\vert v\vert^2)^{(M_0-M)/2}\vert\partial_u^{(\alpha,\beta)+\mu}A(z+u)\vert\vert\partial_v^{(\beta,\alpha)+\nu}B(z+v)\vert
\eea
In particular $ R_N(A,B,z)\in \mathbb S_\delta^{\mu_{AB}}$ for some $\mu_{AB}\geq \mu_A+\mu_B$.
\end{lem}
For proving  this Lemma  one assume first that $A, B\in{\mathcal{S}}(\mathbb{R}^{2d})$. For the general case we put $A_\varepsilon(X)={\rm e}^{-\epsilon\vert X\vert^2}A(X)$,
$B_\varepsilon(X)={\rm e}^{-\epsilon\vert X\vert^2}B(X)$ and pass to the limit for $\varepsilon\searrow 0$.
\section{Proof of Theorem(\ref{MainThm})}
 Here $\hbar=1$. It is enough to consider the test observables $ T_Y:={\rm e}^{-iL_Y}$ ($Y\in\mathbb{R}^{2d}$). \\
 We have
 $$
 \widehat T_Y\widehat H\widehat T_Y^* = [\widehat T_Y, \widehat H]\widehat T_Y^* + \widehat H
 $$
Using the assumption  of  Theorem(\ref{MainThm}) we get
\begin{equation}}   \def\eeq{\end{equation}\label{MPC1}
\frac{1}{i}(\{T_T^*, H\}\circledast T_Y)(X) = H(X+Y)-H(X), \forall X, Y\in\mathbb{R}^{2d}.
\eeq
Computing the Poisson bracket in \eqref{MPC1} gives
\begin{equation}}   \def\eeq{\end{equation}\label{MPC2}
(((y\cdot\partial_x H +\eta\cdot\partial_\xi H)T_Y^*)\circledast T_Y)(X) = H(X+Y) -H(X),\;\forall X, Y\in\mathbb{R}^{2d}.
\eeq
Our aim is to prove that \eqref{MPC2}  implies that $H(X)$ is a polynomial of degree at most 2. For that purpose we shall compute
the asymptotic expansion  as $Y\rightarrow 0$ of the left hand side  of \eqref{MPC2} and compare it with the Taylor expansion for  $H(X+Y)$ modulo
 $O(\vert Y\vert^4)$. From that we shall conclude that  all the third  order derivatives of $H$ vanish in $X$ hence the conclusion will follow.

 We  have
 $$
 \partial_x^\alpha\partial_\xi^\beta T_Y = i^{-\vert\alpha+\beta\vert}\eta^\alpha y^\beta T_Y
 $$
 Let us denote by $C(X,Y)$ the left hand side in \eqref{MPC2}.
 So using Lemma(\ref{lem:Moyalest})   we have
 $$
 C(X,Y) = \sum_{0\leq j\leq 2}(C_j(X,Y) + O(\vert Y\vert^4),
 $$
 where
 \begin{eqnarray}}  \def\eea{\end{eqnarray}
 C_0(X,Y) &= & Y\cdot\nabla_XH(X)\\
 C_1(X,Y) &=& \frac{1}{2}Y\cdot \nabla_X^2H(X)Y,
 \eea
 where $\nabla_X^2H(X)$ is the Hessian matrix of $H$.\\
 Let us compute now $C_2(X,Y)$, which is   an homogeneous polynomials of degree 3 in $Y$.\\
 For simplicity let us consider the 1-D case. The same computation can clearly  be done for $d>1$.\\
 Using \eqref{prod2} we get with $Y=(y,\eta)$,
 \begin{equation}}   \def\eeq{\end{equation}\label{C2}
 C_2(X,Y) = \frac{1}{8}\left(y^3\partial_x^3H +\eta^3\partial_\eta^3H -y^2\eta\partial_\xi\partial_x^2H - y\eta^2\partial_\xi^2\partial_xH\right).
 \eeq
 According \eqref{MPC2}, $C_2(X,Y)$ must  coincide with the term of order 3  in $Y$ of the Taylor expansion in $X$ for $H(X+Y)-H(X)$.
 But this is possible only if $\partial_x^3H=\partial_\eta^3H =\partial_\xi\partial_x^2H =\partial_\xi^2\partial_xH=0$   for any $(x,\xi)\in \mathbb{R}^2$.
 So $H$ must be a polynomial of degree $\leq 2$. $\square$

 \section{Extension  to polynomials of arbitrary degree }
 The asymptotic expansion in $\hbar$ in the Moyal product suggest to introduce the following semi-classical approximations of the Moyal bracket:
 $$
 \{ A,B\}_{\circledast,m} = \{A,B\} +\hbar^2\{A,B\}_3 +\cdots +\hbar^{2m}\{A,B\}_{2m+1},
 $$
 where $\{A,B\}_j = \frac{i}{\hbar}(C_j(A,B)-C_j(B,A))$  (notation of  \eqref{exprod}).
 Notice that $\{A,B\}_j =0$ for $j$ even.\\
 It is clear that if $H$ is a polynomial of degree at most $2m+2$ then we have $ \{ A,H\}_{\circledast,m} = \{ A,H\}_{\circledast}$
 for any $A$.  Conversely we have
\begin{theorem}\label{thm:orderm} Assume $\hbar =1$. If for any $A\in {\mathbb S}^0_0$ we have $ \{ A,H\}_{\circledast,m} = \{ A,H\}_{\circledast}$
 then $H$  must be  a polynomial of degree at most $2m+2$.
 \end{theorem}
 {\em Proof}. Here we give a proof different from the case $m=0$, without connection with the Taylor   formula,  for simpler computations.\\
 Using Lemma \ref{lem:Moyalest} we have
 \begin{equation}}   \def\eeq{\end{equation}\label{comp1}
 T_Y^*\left( \{T_Y,H\}_{\circledast} - \{T_Y,H\}_{\circledast,m} \right)= {\mathcal O}(\vert Y\vert^{2m+3}),\; Y\rightarrow 0.
 \eeq
 Moreover from \eqref{exprod} we get:
 \begin{equation}}   \def\eeq{\end{equation}\label{comp2}
  T_Y^*\{T_Y,H\}_{2j+1} =\frac{1}{2^{j+1}}\sum_{\vert\alpha+\beta=2j+1}\frac{y^\alpha\eta^\beta}{\alpha!\beta!}\partial_x^\alpha\partial_\xi^\beta H
   \eeq
Using the assumption of Theorem \ref{thm:orderm}. and \eqref{comp1}  we get that
$$
 T_Y^*\{T_Y,H\}_{2m+3} = {\mathcal O}(\vert Y\vert^{2m+5})
 $$
 But  $T_Y^*\{T_Y,H\}_{2m+3} $ is an homogeneous polynomial of degree $2m+3$ in $Y$ so we get that this polynomial is 0 and from
 \eqref{comp2} we get that $ \partial_x^\alpha\partial_\xi^\beta H =0 $   for $\vert\alpha+\beta\vert=2m+3$. Then we can conclude that $H(X)$ is a polynomial of degree at most $2m+2$ in $X\in\mathbb{R}^{2d}$. $\square$.


\section*{Introduction}

\subsection*{Motivation}
Optical pump-probe experiments provide a powerful avenue to explore collective dynamics in correlated quantum materials. Several experiments have demonstrated dramatic changes in the terahertz and mid-IR reflectivity following photo-excitation. In experiments on  $\rm{YBa_2Cu_3O_{6.5}}$ cuprate superconductors, parametric resonances caused by internal Josephson plasma oscillations can induce extra "edge" features in the optical reflectivity \cite{Kaiser14,vonHoegen19,Marios20}. On the other hand, the phonon-polariton system $\rm SiC$\cite{Limonov17} and the bulk superconductor $\rm K_3 C_{60}$\cite{budden20, dolgirev21} provide examples where oscillations inside these materials result in significant reflectivity enhancement, possibly exceeding unity. A unified interpretation of these seemingly disparate pump induced features can be formulated in terms of Floquet theory under the assumption that collective excitations create a system with time-periodic properties. Previous analysis has demonstrated the existence of four types of drive induced features in reflectivity, depending on the relative strength of parametric driving and dissipation\cite{FloquetFresnel21}. Edge-like features occur from interference between different Floquet components of the transmission channels when dissipation is strong compared to the oscillation amplitude. In the opposite regime, strong amplification of reflectivity occurs when the parametric drive is not compensated by dissipation and the material exhibits a lasing instability. In fact, such reflectivity features can serve as reporters of a lasing instability, indicating that the effective Floquet medium can be used as a gain medium in a laser.

\begin{figure*}[t!]
    \centering
    \includegraphics[width=1.\textwidth]{Figure1.png}
    \caption{Schematic mechanism that leads to amplification of THz optical reflectivity following high frequency pumping. a) An ultrafast laser pulse (0.5eV, 150 fs) photoexcites electrons through direct dipole transitions between the valence and conduction bands. The photo-excited electrons generate pairs of phonons through the electron-phonon interaction squeezing the phonon field (see equation~(\ref{eq:electronphononInt})). b) Reflectivity amplification of pumped $\rm Ta_2NiSe_5$. Coherent phonon field fluctuations oscillate at twice the phonon frequency, $\omega_{ph}$, creating an effective THz Floquet medium. Parametric driving from the phonon oscillations can create pairs of photons at the signal and idler frequency once stimulated by the probe pulse. This enhances the reflectivity and also scatters counter-propagating light oscillating at the idler frequency. c) Relative change in the reflectivity as a function of frequency subsequent to photoexcitation.  The experimental results are shown at two different temperatures together with the theoretical fit which considers parametric amplification by a 9.4 THz oscillating field. }
    \label{fig:process}
\end{figure*}

Here, we develop a theory to explain the experimentally observed pump-induced terahertz amplification of the reflectivity in layered quasi-1D semiconductor $\rm{Ta_2 Ni Se_5}$, a potential excitonic insulator material. An excitonic insulator phase can emerge when bound electron-hole pairs form a condensate in a semi-metallic phase opening up a gap and causing a phase transition to an insulating state\cite{Jerome1967}. In $\rm{Ta_2 Ni Se_5}$, alternating chains of Ta and Ni are aligned along the axis, forming a sheet in the $\textit{ac}$-plane, and multiple layers are stacked along the $\textit{b}$-axis (see Fig. 1A). An exchange of electron between adjacent Ta and Ni chains along the $\textit{a}$-axis creates an exciton across the chain (along the $\textit{c}$-axis) \cite{kaneko2013,Werdehausen2018,Sugimoto2018,Mazza2020,Lu2017,Larkin2017,Larkin2018,Mor2018,Seki2014}. Previous studies have reported a second order monoclinic-orthorhombic phase transition at $\rm{T_c}$ = 326 K. This leads to an insulating gap below the critical temperature which grows to be 160 meV at low temperatures. The pump-probe experiment demonstrates that above-gap frequency optical excitation induces reflectivity amplification which is maximal at $4.7$ THz, corresponding to an IR-active $\rm{B_{3u}}$ phonon mode confirmed by previous infrared measurements \cite{Larkin2018,Larkin2017}. This turns out to be a surprising result given that the high-frequency pump cannot directly excite the collective THz modes. We show phenomenologically and using ab-initio calculations that this dramatic down conversion is caused by a nontrivial interplay of electron-phonon interactions and phonon non-linearities.

The multi-step process described in this article, converting high frequency pumping to THz amplification of reflectivity, is schematically illustrated in figure~\ref{fig:process}(a) and is outlined as follows: a) The pump excites electrons from the valence bands to the conduction bands through direct dipole transitions. b) The photo-excited electrons exhibit an unusual quadratic coupling to the $4.7$ THz phonons (see Results):
\begin{equation}
    H_{el - ph} = \sum_k g_k n_{el,k} Q^2,
\label{eq:electronphononInt}
\end{equation}
where $Q$ is the IR-phonon displacement, $n_{el,k}$ the photo-excited electron occupation and $g_k$ the effective electron-phonon coupling.
The pump is non-resonant with the IR active terahertz phonons and thus does not directly initiate the enhanced reflectivity dynamics. Instead, IR-phonon pair generation occurs via a Raman process caused by pump induced changes in the electronic occupation. The result is that the expectation value of the phonon displacement is zero, $\expe{Q} = 0$, but the fluctuations are squeezed and coherently oscillate at twice the phonon frequency, $\expe{Q^2(t)} = \expe{Q^2}_0 + A \cos(2\omega_{\rm ph} t) $. c) The squeezed phonon oscillations and phonon nonlinearities create a Floquet material oscillating at $2 \omega_{\rm ph}  = 9.4 $ THz. d) Parametric resonances due to the oscillating field occur primarily around $\omega_{\rm ph}$, parametrically amplifying the reflectivity as depicted in Figure~\ref{fig:process} (b).

The theory developed in this article addresses two inter-connected questions: 1) Can a Floquet material oscillating at 9.4 THz provide parametric amplification similar to experimental observations? 2) What is the origin of the 9.4 THz oscillation?.

\subsection*{Summary of results}
We begin our analysis by showing that the parametrically amplified reflectivity at $4.7$ THz is consistent with the presence of a coherently oscillating Raman mode at twice the frequency, $9.4$ THz. The theoretical fit captures the experimental amplification profile, plotted in figure~\ref{fig:process}(c). A schematic Hamiltonian describing parametric amplification of the optical reflectivity in the presence of a coherently oscillating mode at $\omega_d = 9.4$ THz is given by:
\begin{equation}
    H_{ampl.} = X(t) a_{1,k}^\dag a_{2,-k}^\dag + h.c.,
\end{equation}
where $X(t)$ is a Raman mode and $ \{ a^\dag_{1,k} , a^\dag_{2,-k} \} $ are the creation operators of photons with opposite momentum which may be associated with different phonon-polariton bands $\{1, 2\}$ at frequencies parametrically resonant with the drive, $\omega_{1} (k) + \omega_2 ( - k)  = \omega_d$. As shown in figure~\ref{fig:process}(b), Floquet drive can generate pairs of photons at the signal frequency of the incoming probe, $\omega_s$, and a photon oscillating at the idler frequency, $\omega_{\rm id} = \omega_d - \omega_s$. Reflectivity amplification is then caused by stimulated emission of photon pairs by the driven material. To compute the reflection coefficient, we use degenerate Floquet perturbation theory to construct the eigenmodes and solve the Fresnel equations at the boundary.

We then proceed to investigate the origin of the $9.4$ THz oscillation. We show that electron-phonon coupling between electronic bands involved in the photo-excitation and the 4.7 THz IR-phonon naturally leads to phonon squeezing. Phonon squeezing oscillations have twice the frequency of the phonon and act as a Raman mode. We find that this mechanism is generic and should apply to all IR-phonons present in the system. However, for phonon squeezing oscillations to create an effective Floquet material that can significantly enhance the reflectivity, the phonon needs to be both strongly coupled to the photo-excited bands and also have a strong coupling to the electric field.

We hence perform ab-intio calculations to determine both the magnitude of the electron phonon coupling and the IR activity of the different phonon modes. We find that, indeed, the $4.7$ THz IR-phonon is special, exhibiting strong electron-phonon coupling and appreciable IR activity as compared to other IR-phonons in the same frequency range. Moreover, we confirm that phonon fluctuations are responsible for the THz amplification. The theory shows that the predicted reflectivity amplification is a reporter of a lasing instability, implying that pumped $\rm Ta_2 Ni Se_5$ could, in principle, serve as a Floquet gain medium for THz lasing.

As a final intriguing aspect, we find that the $4.7$ THz mode is strongly coupled to the elusive and hotly debated order parameter of $\rm Ta_2 Ni Se_{5}$\cite{Baldini20,Sugimoto2018,Mazza2020,Golez2020,Sun2020,Subedi2020,kaneko2013}, which is thought to have some excitonic insulating character. A number of state-of-the-art experiments have been performed to shed light on the origin of the phase transition at $326$ K as well as what role the excitonic and lattice instabilities play in this phase transition process. However, the results show contradictory results \cite{Werdehausen2018,Werdehausenagain2021,Lu2017,Larkin2017,Larkin2018,Mor2017,Mor2018,Seki2014,Kim2020,Kim2021,Bretscher2021,Bretscheragain2021,Baldini20,Volkov2021,Pal2020,Tang2020,Chen2020,Seo2018,Wakisaka2009,Nakano2018,Saha2021,Lee2019,Suzuki2021} and the question whether $\rm Ta_2 Ni Se_{5}$ hosts an excitonic insulator phase remains an open question. Using DFT frozen phonon calculations, we find that the electron-phonon coupling for the $4.7$ THz phonon effectively vanishes in the high temperature orthorhombic phase. This indicates that the THz parametric amplification of reflectivity is mediated by the low temperature phase and is expected to be sensitive to the phase transition. This interplay between the order parameter and phase transition has been observed in pump-probe \cite{Werdehausen2018,Werdehausenagain2021}, near infrared \cite{Bretscher2021,Bretscheragain2021}, and time-resolved ARPES experiments \cite{Baldini20,Tang2020,Suzuki2021,Denis2021} , and has unraveled key ingredients in the physics of $\rm{Ta_2NiSe_5}$ in addition to providing a clearer understanding of the nature of the order parameter.  The link between the 4.7 THz phonon and the order parameter opens new avenues for using parametric amplification as a novel method to track order parameter dynamics.

\section*{Results}

\subsection*{Parametric amplification of reflectivity in pumped $\rm{Ta_2 Ni Se_5}$}
\label{sec:paraAmpl}

\subsubsection{Equilibrium reflectivity}
The reflectivity of a material is captured by the frequency dependent refractive index appearing in the Maxwell equations:
\begin{equation}
    \left( \frac{n(\omega)^2\omega^2}{c^2} - k^2 \right) E = 0
\end{equation}
where $E$ is the electric field and all information about phonons and other IR active modes is encoded in $n(\omega)$. We assume that the probe corresponds to an electromagnetic wave reflected from the sample at normal incidence. The propagation direction is along the $b$-axis of the crystal, which we refer to as the $y$-direction, whereas the electric field points in the $a$-direction which we refer to as the $z$-axis. The refractive index can be directly extracted experimentally from the complex reflection coefficient at normal incidence in an equilibrium system\cite{WOOTEN72}:
\begin{equation}
    n(\omega) = \frac{1 - r (\omega)}{1 + r (\omega)}.
\label{eq:indexRefl}
\end{equation}

\subsubsection{Eigenstates in Floquet state}
Once we have obtained the refractive index from the equilibrium reflectivity, we model the Floquet material as experiencing a parametric drive oscillating at frequency, $\omega_d$, which mixes signal and idler frequencies:
\begin{equation}
    E(t) =  e^{ i k y} \left( E_s e^{ - i \omega_s t} + E_{id}e^{i \omega_{id} t} \right).
\label{eq:ansatzTNS}
\end{equation}
where $\omega_s$ is the frequency of the incoming probe and $\omega_{id} = \omega_d - \omega_s$. Such mixing corresponds to degenerate perturbation theory in Floquet systems and the idler component is the nearest Floquet band contribution to $\omega_s$ which is responsible for parametric instabilities \cite{Marios20, Sho19, Eckardt15}. The oscillating mode is included phenomenologically through a time-periodic contribution to the electric permittivity:
\begin{equation}
    \delta \epsilon(t) = 2 A_{\rm drive} \cos(\omega_d t).
\end{equation}
Using the ansatz in equation~(\ref{eq:ansatzTNS}), the equations of motion in the Floquet state for the different oscillating components of the electric field become:
\begin{subequations}
\begin{align}
\left( \frac{n(\omega_s)^2\omega_s^2}{c^2} - k^2 \right) E_s + A_{dr} E_{id} =& 0 \\
\left( \frac{n(\omega_{id})^2\omega_{id}^2}{c^2} - k^2 \right) E_{id} + A_{dr} E_s  =& 0.
\end{align}
\label{eq:FloquetEigen}
\end{subequations}
To compute the reflectivity at normal incidence, we first find the allowed $k$ values for a given $\omega_s$. Due to the coupling of signal and idler components, two such $k$ values exist, associated with two transmission channels, both of which oscillate at signal and idler frequencies. The transmission channels correspond to eigenvectors of the Floquet equations of motion inside the material,
\begin{equation}
    E_{i} = t_i E_0 e^{i k_i y} \left(e^{  - i \omega_s t}  + \alpha_i e^{ i \omega_{id} t} \right),
\end{equation}
where $\alpha_i$ is the relative amplitude of the signal and idler component in the eigenvector of wave-vector $k_i$, $t_i$ is the transmission coefficient of the $i$-th channel and $E_0$ the amplitude of the incoming field.

\subsubsection{Floquet-Fresnel equations}
In a reflection problem, the eigenvalue equation~(\ref{eq:FloquetEigen}) enables computation of the transmitted wavevectors as a function of a fixed frequency set by the incoming light. However, the answer is given in terms of $k_i^2$ rather than $k_i$ and the correct root is chosen such that the field vanishes at infinity, $\mbox{Im} \{k_i \} > 0  $. The reflectivity is computed by solving the Floquet-Fresnel equations, matching electric and magnetic fields parallel to the surface both at frequency $\omega_s$ and at frequency $\omega_{id}$. Inside the material we have the electric field:
\begin{equation}
    E_{mat} = E_0 \sum_i t_i e^{i k_i y} \left( e^{- i \omega_s t} + \alpha_i e^{ i \omega_{id} t} \right) ,
\end{equation}
and for vacuum we have:
\begin{equation}
\begin{split}
    E_{vac} = &E_0 \left(e^{ i \omega_{s}/c y - i \omega_s t} + r_s e^{- i \omega_s/c y - i \omega_s t} \right) +\\
    &E_0 r_{id} e^{  i \omega_{id}t. + i \omega_{id} y }.
\end{split}
\end{equation}
Using the homogeneous Maxwell equations, $ \nabla \times E = - \partial_t B$ to compute the magnetic field and matching boundary conditions at $y = 0$, we obtain the Fresnel equations for the driven system:
\begin{subequations}
\begin{align}
1 + r_s =& t_1 + t_2, \\
1 - r_s =& \frac{k_1}{\omega_s} t_1 + \frac{k_2}{\omega_s} t_2, \\
r_{id} =& t_1 \alpha_1 + t_2 \alpha_2, \\
r_{id} =& \frac{ k_1}{\omega_{id}} t_1 \alpha_1 + \frac{k_2}{ \omega_{id} } t_2 \alpha_2.
\end{align}
\label{eq:FloquetFresnel}
\end{subequations}
To fit the data, we choose a drive at $9.4$ THz and allow for small changes in the static properties of the system such as a photo-induced conductivity stemming from photo-excited charge carriers. The fitted parameters are given in the Materials and Methods section.

\subsection*{ Phonon squeezing initiated by photoexcited electrons}
\label{sec:eltoPhon}

\subsubsection{Electron phonon coupling}
In this section, we develop a microscopic theory of the coupling of electronic bands to IR active phonons. In the dipolar gauge\cite{Jiajun2020}, the effective dipole of the phonon is linearly coupled to the electronic dipole. For two electronic bands with an allowed dipole transition, this coupling takes the form:
\begin{equation}
    H_{\rm el- ph} = Q \sum_k \lambda_k \left( \hat{c}^\dag_{1,k} \hat{c}_{2,k} + \hat{c}^\dag_{2,k} \hat{c}_{1,k} \right),
\end{equation}
where $Q$ is the phonon coordinate, $\hat{c}_1$ and $\hat{c}_2$ the annihilation operators of the two electronic bands and $\lambda_k$ the electron-phonon coupling matrix element as a function of momentum $k$. Using a Schrieffer-Wolff transformation, we "integrate out" the linear electron-phonon coupling, which is non-resonant due to the different energy scales between the IR-phonon and the electronic transition, to reveal the resonant non-linear coupling between the electron occupation number and the phonon squeezing operator, $Q^2$. In the materials and methods material, we show that this procedure leads to the effective coupling:
\begin{equation}
    H_{\rm el-ph, eff.} = \sum_k \frac{\lambda_k^2}{ \Delta_k}  \left( n_{1,k} - n_{2,k} \right) Q^2,
\label{eq:Squeezingcoupling}
\end{equation}
where $\Delta_k = E_{1,k} - E_{2,k}$ is the energy difference between the electronic states and $\hat{n}_{i,k} = \hat{c}^\dag_{i,k} \hat{c}_{i,k } $ the number operator.  The above equation applies to phonons that are coupled to independent pairs of electronic bands. If three or more bands are simultaneously coupled to a specific phonon the effective electron-phonon Hamiltonian is more complicated but with similar qualitative features, such as the  coupling of electron density to the square of the phonon coordinate.

\subsubsection{Phonon squeezing}

Once electrons have been photo-excited, the finite number of optically excited electrons quenches the frequency of the phonon:
\begin{equation}
    H_{\rm el-ph, eff.} = \frac{M f(t) Q^2}{2},
\label{eq:ParametricPhonon}
\end{equation}
where $M$ is the mass of the phonon. In equation~(\ref{eq:ParametricPhonon}), the parametric driving $f(t)$ comes from the photo-excited distribution of electrons that is strongly coupled to the phonon:
\begin{equation}
    f(t) = \sum_k \frac{2 \lambda_k^2}{M \Delta_k} \left( \expe{n_{1,k}}(t) - \expe{n_{2,k}(t) }\right),
\end{equation}
The photo-excited electron dynamics are fast and can be approximated as a delta function
in time. More generally we can approximate the photo-excited distribution as having a characteristic life time, $f(t) = f_0 \theta(t)e^{-t/t_{decay}} $, but this does not change our conclusions.

The above electron induced drive describes a Raman process that does not excite the phonon directly (i.e. $\expe{Q} = 0$). However, the squeezing operator starts oscillating at frequency equal to twice the phonon frequency, as shown in the materials and methods:
\begin{equation}
   \left( \partial_t^2 - \left( 2 \omega_{\rm ph } \right)^2  \right) \expe{Q^2} = - 4 f(t) \expe{Q^2}_0,
\end{equation}
where $\expe{Q^2}_0$ is the equilibrium fluctuations and $\omega_{\rm ph}$ is phonon frequency at zero momentum which is renormalized by the coulomb force, $\omega_{\rm ph}^2 = \omega_{\rm ph,0}^2 + \frac{Z^2}{\epsilon \epsilon_0 M }$.

To show that this phenomenon is related to squeezing, we expand $Q^2$ in terms of creation and annihilation operators, using $Q = \frac{a + a^\dag}{\sqrt{2 M \omega_{ph}}}$:
\begin{equation}
    Q^2(t) = \frac{1}{2 M \omega_{ph}} \left( a^\dag(t) a^\dag(t) + a(t) a(t) + a^\dag(t) a(t) + a(t) a^\dag(t) \right).
\end{equation}
Expectation values of $a(t) a(t)^\dag$ do not oscillate rapidly while the anomalous pairs $a^\dag(t) a^\dag(t)$ and $a(t) a(t)$ oscillate at twice the phonon frequency. As a result a state with phonon fluctuations, $Q^2$, that oscillate at twice the phonon frequency implies the existence of a condensate of phonon pairs $\expe{a^\dag(t) a^\dag(t)} \neq 0$.

\subsubsection{Floquet matter}

Coherent oscillations of long lived modes, such as the phonon squeezing oscillations, turn matter into Floquet matter with time-periodic properties through non-linear interactions with the rest of the material. Following the same steps as in references \cite{Sho19,Michael21}, we show in the materials and methods how a lattice anharmonicity  proportional to $Q^4$ and squeezing dynamics lead to the signal and idler coupling in equation~(\ref{eq:FloquetEigen}). For this model, we find that the signal idler coupling is proportional to the IR activity of the mode, $Z$, the phonon anharmonicity, $u$, and the amplitude of the oscillations, $B$ : $A_{\rm drive} \propto Z^2 u B$.

 \subsection*{Ab-initio calculations : the 4.7 THz IR-Phonon}
 \label{sec:DFT}
 \begin{figure*}
    \centering
    \includegraphics[trim = 0pt 0pt 0pt 0pt,clip, scale = 0.4]{finalPlotv3.pdf}
    \caption{Summary of ab-initio calculations. a) Optical matrix elements along the \textit{a}-axis calculated for the frequency region between 0.4 to 0.9 eV. The size of the circles is proportional to the magnitude of the matrix elements and indicates which electronic valence and conduction bands are involved in photo-absorption. b) IR activity of the most IR active phonons; the phonon at 4.7 THz which dominates electron-phonon interactions in that region is highlighted with the mode character shown in the inset. A complete list of IR active phonons in this region for TNS is given in materials and methods. c) and d) show the effect of electron-phonon coupling on the band structure in the monoclinic and orthorhombic phase respectively. The electron-phonon coupling is captured by calculating the electronic band-structures using DFT in the presence of either negative or positive phonon displacement (frozen phonon approximation). The phonon 4.7 THz in the monoclinic phase is strongly coupled to the electronic bands 2 and 3 which are highlighted in the figure. In the high temperature orthorhombic phase the effective electron-phonon coupling is reduced by an order of magnitude. This indicates that the electron-phonon interaction is very sensitive to the order parameter.}
    \label{fig:DFTsummary}
\end{figure*}
The above discussion is generic and in principle applies to every IR-phonon inside a given material. In this section we use DFT calculations to determine which phonons make the dominant contribution to the parametric amplification of reflectivity in $\rm Ta_2 Ni Se_5$.
We start our discussion by identifying which valence and conduction bands are involved in photo-absorption by evaluating the optical dipole transitions matrix elements. In Fig.~\ref{fig:DFTsummary}(a), we plot the momentum resolved optical contribution, for field polarization along the \textit{a}-axis, of each band to dipole allowed transitions within the experimentally relevant energy window between 0.33 and 0.9 eV set by the pump parameters. For a given valence(conduction) band and k-point, the optical contribution is defined by summing the square of the dipole transition matrix elements associated with the transition. We find that it is the first three conduction bands that are predominantly excited by the pump.

Turning our attention to phonons, we use ab-initio calculations to identify all IR active phonons. In particular, we find a number of phonons in the low temperature monoclinic phase between 4 and 5 THz shown in the supplementary material. In Fig.~\ref{fig:DFTsummary} (b) we plot the IR activity of phonons as a function of frequency and identify the 4.7 THz phonon which, as shown below, turns out to be the most strongly electron coupled phonon.

We compute the electron-phonon interaction between IR active phonons and electrons. This is to identify the phonons that dominate the interaction with the photo-excited electronic bands. To accomplish this we use the method of frozen phonons. In this approach, the electronic bands are recalculated with the lattice shifted along a phonon eigendisplacement. To quantify the coupling strength of a specific phonon to a particular electronic band we integrate the energy changes of the band in the presence of the phonon over the Brillouin zone as outlined in the supplementary material. We find that in the vicinity of $4.5$ THz, relevant to the experimental observables, the $4.7$ THz mode has roughly an order of magnitude larger electron-phonon coupling strength compared to the nearby phonons. In particular, as shown in Fig.~\ref{fig:DFTsummary}(c), the 2nd and 3rd conduction bands (in this paper we number the conduction bands from lowest to highest in energy) are significantly renormalised in the presence of the 4.7 THz phonon. Since the phonon is mostly coupled to two electronic bands, the frozen phonon calculations are consistent with equation~(\ref{eq:Squeezingcoupling}) leading to the two bands shifting in energy by an equal and opposite amount given by:
\begin{equation}
\begin{split}
   \left.E_{k,3}\right|_{\expe{Q} \neq 0} -  \left.E_{k,3}\right|_{\expe{Q} = 0} =&  - \left(\left.E_{k,2}\right|_{\expe{Q} \neq 0} -  \left.E_{k,2}\right|_{\expe{Q} = 0}\right)\\
   =& \expe{Q}^2 \frac{\lambda_k^2}{\Delta_k \left( 1- \frac{\omega_{ph,0}^2}{\Delta^2_k} \right)}
\end{split}
\end{equation}
where $E_{k,2}$ and $E_{k,3}$ is energy of conduction band two and three at momentum $k$, and $\Delta_k$ is the energy difference of the two bands in equilibrium. Therefore, in this case, calculating the energy shifts as a function of momentum in the frozen phonon approximation allows for direct computation of the electron-phonon interaction.

To summarize, we use ab-initio calculations to identify the electron bands excited by
pumping and to establish the IR-active phonons in $\rm Ta_2 Ni Se_5$. Subsequently, frozen phonon calculations allowed for the identification of the phonon with the dominant electron-phonon coupling to the excited electronic bands. This leads to the identification of the $4.7$ THz mode as responsible for the parametric amplification observed in experiments through the mechanism outlined in previous subsections. We note that we do not exclude the possibility of subdominant amplification in reflectivity spectra arising from other phonons at different frequencies. In the supplementary material, we provide more details of our DFT analysis of IR-phonons in $\rm Ta_2 Ni Se_5$, the frozen phonon calculations and the calculation of optical matrix elements.

\subsubsection{Connection to the order parameter}

Finally, we discuss the effects of the phase transition on the phonon squeezing process. In figure~\ref{fig:DFTsummary}, we compute the electron-phonon interaction of the $4.7$ THz phonon in the low temperature monoclinic phase while for the high temperature orthorhombic phase we compute the electron-phonon interaction of the phonon adiabatically connected to the $4.7$ THz phonon eigenstate (for details on the identification of the adiabatically connected phonon see supplementary material). The electron-phonon interaction effectively disappears in the high temperature phase providing evidence that electron-phonon coupling is mediated by the order parameter. In particular, we suggest that close to the phase transition the electron-phonon coupling is proportional to the order parameter, $ \lambda_k = \Phi B_k$, where $\Phi$ is the order parameter and $B_k$ a constant. This, in turn, suggests that parametric amplification can be used as a nontrivial probe to investigate order parameter dynamics. The microscopic reason for the strong dependence of the electron-phonon coupling to the order parameter is a very interesting question that could reveal new insights about the nature of the order parameter in $\rm Ta_2 Ni Se_5$ and will be addressed in a subsequent publication.

\section*{Discussion}

We have investigated the microscopic mechanism of amplification of THz optical reflectivity in $\rm{Ta_2NiSe_5}$ arising from high frequency optical pumping. We showed that strong electron-phonon coupling opens new pathways towards realizing THz parametric amplification through high frequency pumping. Ab-initio calculations highlight the importance of the $4.7$ THz IR-active phonon which is strongly coupled to electrons allowing for the amplification to manifest in the THz reflection spectrum.

Our theory indicates that choosing which electronic band to photoexcite selects the IR-phonon that is most strongly coupled to that electronic band. As a result, we can use different pumping frequencies in the same material to tune the frequency of the THz parametric amplification through mode-selective phonon squeezing.

Finally, we showed that the electron-phonon coupling is strongly dependent on the order parameter and becomes suppressed in the high temperature orthorhombic phase. This suggests that THz amplification of reflectivity can be used as a new probe to order parameter dynamics.

\section*{Materials and Methods}

\subsection*{Theoretical fit of the parametric reflectivity amplification in pumped $\rm Ta_2 Ni Se_5$}
\label{app:TNSfit}

As mentioned in the main text in equation~(\ref{eq:indexRefl}), in principle the complex refractive index can be directly computed by the complex reflectivity amplitude, $r(\omega)$. However,  small phase errors upon experimental extraction of $r(\omega)$ could lead to unphysical behaviour of the reflectivity. To overcome this complication, we fit the data by assuming that for an insulator like $\rm Ta_2NiSe_5$ the refractive index is real. As a result, we can instead use the absolute value of $|r| = \sqrt{R}$ and express the refractive index as:
\begin{equation}
    n(\omega) = \frac{1 + \sqrt{R}}{1- \sqrt{R}}.
\end{equation}
Upon parametric resonance and solving equation~(\ref{eq:FloquetEigen}) together with the Floquet-Fresnel boundary conditions in equation~(\ref{eq:FloquetFresnel}), assuming a purely real spectrum with no dissipation leads to highly divergent behavior on parametric resonance. As a result, we put back dissipation by including an imaginary component to the refractive index. In the driven case we argue that this is physical, since $n_{\rm driven} (\omega)^2 = n^2(\omega) + i \sigma_{\rm driven}/\omega$ and dissipation can arise through through a transient contribution to the conductivity by electrons excited across the gap of the insulator. To fit the data we choose a parametric drive at $\omega_d = 9.4$ THz, drive amplitude $A_{\rm drive} = 7.5 \frac{THz^2}{c^2}$, overall constant renormalization of the refractive, $n^2_{\rm drive} = n^2(\omega) + \delta n^2 $, with $\delta n^2 = 0.1 + 0.01 i $ and an overall Gaussian broadening function with standard deviation of $0.1$ THz.

\subsection*{Electron - phonon interaction }
\label{app:SchriefferWolff}

The Hamiltonian of two electronic bands coupled by an allowed direct dipole transition that is, in turn, coupled to a phonon is given by:
\begin{equation}
    H_{\rm el}  = \frac{\Delta}{2} \left(\hat{c}_1^\dag \hat{c}_1 - \hat{c}_2^\dag \hat{c}_2 \right) + \lambda \left( \hat{c}_1^\dag \hat{c}_2+ \hat{c}_2^\dag \hat{c}_1 \right) Q,
\label{eq:elecHamil}
\end{equation}
where $\Delta = E_1 - E_2$ is the difference in energy between the two bands, $\lambda$ is the coupling constant coming from the dipole-dipole interaction in the dipole gauge, $Q$ is the phonon coordinate and $\{\hat{c}_1, \hat{c}_2\} $ are the annihilation operators of the two electron bands. The IR-phonon quadratic hamiltonian is given by:
\begin{equation*}
    H_{\rm ph} = Z E Q + M \omega_{\rm ph,0}^2 \frac{Q^2}{2} + \frac{\Pi^2}{2 M },
\label{eq:phononHamil}
\end{equation*}
where $\Pi$ is the conjugate momentum of the phonon coordinate $Q$, $Z$ the effective coupling to the electromagnetic field, E, and $\omega_{\rm ph,0} $ the phonon frequency. Since the transition itself is not resonant with the phonon mode, we decouple the linear electron phonon interaction perturbatively using a Schrieffer-Wolff transformation. This generates an interaction between the electron-hole pair density and the fluctuations of the phonon field which can be a resonant process. To perform this transformation, it is convenient to note that bilinear combinations of $\{ \hat{c}_1^\dag, \hat{c}_1, \hat{c}_2^\dag, \hat{c}_2 \}$ appearing in the Hamiltonian, obey SU(2) commutation relations by making the following identification:
\begin{subequations}
\begin{align}
    S^x = \frac{\hat{c}_1^\dag \hat{c}_2 + \hat{c}_2^\dag \hat{c}_1}{2}, \\
    S^y = -i \frac{\hat{c}_1^\dag \hat{c}_2 - \hat{c}_2^\dag \hat{c}_1}{2}, \\
    S^z = \frac{\hat{c}_1^\dag \hat{c}_1 - \hat{c}_2^\dag \hat{c}_2}{2},
\end{align}
\label{eq:spin}
\end{subequations}
where their commutators are $\left[S^x, S^y\right]  = i S^z$ and its cyclic permutations. In terms of the spin operators, we separate the Hamiltonian into a non-interacting and an  interacting part:
\begin{align}
    H_0 =& \Delta S^z + M \omega_{ph,0}^2 \frac{ Q^2}{2} + \frac{\Pi^2}{2 M }, \\
    V =& 2 \lambda S^x Q
\end{align}
To remove the interacting part $V$ to linear order in $\lambda$, we consider a unitary transformation of the type:
\begin{align}
    U =& {\rm Exp} \{ i A\}, \\
    A =& \alpha Q S^y + \beta \Pi S^x
\end{align}
The Schrieffer-Wolff expansion is given by:
\begin{equation}
    U H_{total}U^\dag = H_0 + V - i \left[H_0, A \right] - i \left[V , A \right] - \frac{1}{2} \left[\left[H_0, A\right],A \right].
\end{equation}
The parameters $\alpha$ and $\beta$ are found such that:
\begin{equation}
    V = i \left[ H_0 , A\right]
\end{equation}
Matching linear terms in $\Pi$ and $Q$, leads to the parameters:
\begin{subequations}
\begin{align}
    &\beta =  \frac{\alpha}{M \Delta }, \\
    -&\Delta \alpha + 2 \lambda + M \omega^2_{\rm ph,0} \beta = 0, \\
    \Rightarrow &\alpha = \frac{2 \lambda}{\Delta \left( 1 - \frac{\omega_{\rm ph,0}^2}{\Delta^2} \right)},
\end{align}
\end{subequations}
which confirms that this perturbation theory can be carried out as long as $\omega_{\rm ph, 0}$ is off-resonant with the transition energy $\Delta$.  The effective electron-phonon interaction after the Schrieffer-Wolff transformation is given by:
\begin{equation}
\begin{split}
    H_{\rm eff} =& - i\left[V, A \right] - \frac{1}{2} \left[\left[H_0, A\right] ,A\right] = - \frac{i}{2} \left[V, A\right],\\
    =& \frac{2 \lambda^2}{\Delta \left(1 - \frac{\omega^2_{\rm ph,0}}{\Delta^2} \right)}Q^2 S^z + \frac{2 \lambda^2}{\Delta \left( 1 - \frac{\omega^2_{\rm ph,0} }{\Delta^2} \right)} \left(S^x\right)^2,
\end{split}
\end{equation}
where the second term does not depend on the phonons. The residual electron phonon interaction is thus give by:
\begin{equation}
    H_{\rm el-ph, eff} = \frac{\lambda^2}{\Delta \left( 1 - \frac{\omega_{\rm ph, 0}^2}{\Delta^2} \right)} Q^2 \left( n_1 - n_2 \right)
\end{equation}
where $n_1 = \hat{c}_1^\dag \hat{c}_1$ is the occupation number of mode 1. This implies that if either mode 1 or mode 2 are photo-excited the phonon can be squeezed.

The linear term in the electric field, $Z E Q$, appearing in equation~\ref{eq:phononHamil}, is also transformed by the Schrieffer-Wolff transformation. It gives rise to a term $\beta Z E S^x$ which provides a small renormalization of the dipole transition amplitude between the bands and does not affect our discussion.

The above result can be generalized to an arbitrary number of independent pairs of electronic states. For example, for two electronic states the above result in the limit of $\frac{\omega_{ph,0}^2}{\Delta_k^2} \rightarrow 0 $ is generalized to:
\begin{equation}
     H_{\rm el - ph, eff} = \sum_k \frac{\lambda_k^2}{\Delta_k } Q^2 \left( n_{1,k} - n_{2,k}\right),
\label{eq:elphonEff}
\end{equation}
This expression is the one used in the main text.

\subsection*{Phonon squeezing}
\label{app:phononSqueeze}
Using a Hartree-Fock type approximation on the effective electron-phonon Hamiltonian in equation~(\ref{eq:elphonEff}), we derive an effective Hamiltonian for the phonon system only:
\begin{equation}
    H_{\rm ph} = Z E Q + M \left( \omega_{\rm ph, 0 }^2 + f(t) \right) \frac{Q^2}{2} + \frac{\Pi^2}{2 M }
\label{eq:phononEff}
\end{equation}
where $Z$ is the IR activity of the phonon mode, $\omega_{\rm ph,0}$ the bare phonon frequency and the effective parametric drive, $f(t)$, is given by photo-excited electron density coupled to the phonon mode:
\begin{equation}
    f(t) = \sum_k \frac{2 \lambda_k^2}{M \Delta_k } \left( \expe{n_{1,k}} - \expe{n_{2,k} }\right),
\end{equation}
Due to the fast dynamics of electrons, $f(t)$ acts as an impulsive delta function like parametric drive. Such a drive is not periodic but it can linearly excite phonon fluctuations $\expe{Q^2}$ which will oscillate in time. To show this, we compute the equations of motion for fluctuations of the phonon field:
\begin{subequations}
\begin{align}
\partial_t \expe{Q^2} =& \frac{\expe{ \Pi Q+ Q \Pi}}{M} , \\
\begin{split}
    \partial_t \expe{ \Pi Q + Q \Pi}=& -2 M \left( \omega^2_{\rm ph,0} + f(t) \right) \expe{ Q^2} + \\& 2\frac{\expe{ \Pi^2}}{M} - 2 Z \expe{E Q} ,
\end{split}\\
\begin{split}
\partial_t \expe{\Pi ^2} =& - M \left( \omega^2_{\rm ph, 0} + f(t) \right)\expe{ \left( \Pi Q+ Q \Pi \right) } - \\&Z \expe{E \Pi + \Pi E} .
\end{split}
\end{align}
\end{subequations}
At $k = 0$, we can use Maxwell's equations to remove the electric field dependence, $E = \frac{Z}{\epsilon \epsilon_0 } Q $. Performing this substitution simplifies the equations of motion,
\begin{subequations}
\begin{align}
\partial_t \expe{Q^2} =& \frac{\expe{ \Pi Q + Q \Pi}}{M} , \\
    \partial_t \expe{ \Pi Q + Q \Pi}=& -2 M \left( \omega^2_{\rm ph} + f(t) \right) \expe{ Q^2} +  2\frac{\expe{ \Pi^2}}{M} ,\\
\partial_t \expe{\Pi ^2} =& - M \left( \omega^2_{\rm ph} + f(t) \right) \expe{  \left( \Pi Q+ Q \Pi \right) }.
\end{align}
\label{eq:fluceom}
\end{subequations}
where $\omega_{\rm ph}^2 = \omega^2_{\rm ph,0} + \omega_{\rm pl, phonon}^2$ is the frequency of the phonon at $k = 0$ which differs from the bare frequency by the phonon plasma frequency given by $\omega^2_{\rm pl, phonon} = \frac{Z^2}{\epsilon \epsilon_0 M }$. Being perturbative in the drive $f(t)$, we expand the phonon fluctuations,
\begin{equation}
    \expe{Q^2} = \expe{Q^2}_0 + \expe{Q^2}_1,
\end{equation}
where $\expe{Q^2}_0$ is the thermal expectation value and $\expe{Q^2}_1 \propto f(t)$. To linear order in $f(t)$, the equations of motion imply, $\expe{\Pi^2}_1/ M = M \omega^2_{\rm ph} \expe{Q^2}_1 + M f(t) \expe{Q^2}_0 + \mathcal{O}(f^2) $. Finally, combining equations in (\ref{eq:fluceom}) we find that:
\begin{equation}
   \left( \partial_t^2 + 4 \omega_{\rm ph}^2 \right) \expe{Q^2}_1 = - 4 f(t) \expe{Q^2}_0.
\end{equation}
This result shows that phonon fluctuations are linearly driven by photo-excitation and oscillate at twice the phonon frequency, $2 \omega_{\rm ph}$. These coherent oscillations of phonon fluctuations behave as a Raman mode.

\subsection*{Floquet matter from squeezing dynamics of phonons}
\label{app:IRphonon}

Material properties such as the electric permittivity, become time-periodic in the presence of oscillating fields through interactions. Here, we demonstrate how lattice potential anharmonicities lead to a time-periodic index of refraction:

We consider a phonon system with a $Q^4$ anharmonicity for the IR-phonon with a Hamiltonian:
\begin{equation}
     H_{\rm ph} = Z E Q + M \left( \omega_{\rm ph, 0 }^2 + f(t) \right) \frac{Q^2}{2} + \frac{\Pi^2}{2 M } + u Q^4
\label{eq:NonlinHamil}
\end{equation}
The equations of motion for the phonon given by the Hamiltonian in equation~(\ref{eq:NonlinHamil}) is
\begin{equation}
    \left( \partial_t^2 + \gamma \partial_t + \omega_{\rm ph,0}^2 + 4 uQ^2 \right) Q = \frac{Z}{M} E.
\end{equation}
Using a gaussian ansatz for the phonons, we can linearize the above equation as:
\begin{equation}
     \left( \partial_t^2 + \gamma \partial_t + \omega_{\rm ph,0}^2 + 12 u \expe{Q^2}(t) \right) Q = \frac{Z}{M} E,
\label{eq:phononeom}
\end{equation}
where the fluctuations $\expe{Q^2} = \expe{Q^2}_0 + A \cos( 2 \omega_{ph} t )$.
The phonon mode appears in the Maxwells equations as:
\begin{equation}
    \left( \frac{1}{c^2}\partial_t^2 - k^2 \right) E = - Z \partial^2_t Q.
\end{equation}
To find the effective signal idler mixing presented in equation~(\ref{eq:FloquetEigen}), we expand the equations of motion in signal and idler contributions:
\begin{equation}
    Q = Q_{s} e^{-i \omega_s t} + Q_{\rm id} e^{i \omega_{\rm id} t}
\end{equation}
Equation~(\ref{eq:phononeom}) becomes:
\begin{equation}
\begin{split}
&\begin{pmatrix} Q_s \\ Q_{id} \end{pmatrix} = \begin{pmatrix} \frac{Z}{ \omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 } && 0 \\ 0 &&\frac{Z}{ \omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2}  \end{pmatrix} \cdot \begin{pmatrix} E_{s} \\ E_{\rm id}  \end{pmatrix} + \\
&\frac{Z u A}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2 \right) } \begin{pmatrix} E_{\rm id} \\ E_s \end{pmatrix} .
\end{split}
\label{eq:Floqphonon}
\end{equation}
Substituting equation~\ref{eq:Floqphonon} in Maxwells equation we find the equations of motion for the signal and idler component of the electric field to be:
\newline
\begin{subequations}
\begin{align}
\left( \frac{n^2_{eq.} (\omega_s) }{c^2} \omega_s^2 - k^2 \right) E_s+  A_{\rm drive, s} (\omega_s, \omega_{\rm id} ) E_{\rm id}=&0 ,\\
\left( \frac{n^2_{eq.} (\omega_{\rm id}) }{c^2} \omega_s^2 - k^2 \right) E_s+  A_{\rm drive, id} (\omega_s, \omega_{\rm id} ) E_{\rm id}=&0
\end{align}
\end{subequations}
where the signal and idler driving amplitudes, $A_{\rm drive, s}$ and $A_{\rm drive, id}$ are given by:
\begin{subequations}
\begin{align}
\begin{split}
&A_{\rm drive,s} =\\
&\frac{Z^2 u  A \omega_s^2}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2 \right) } ,
\end{split} \\
\begin{split}
     &A_{\rm drive,id} = \\
     &\frac{Z^2 u  A \omega_{\rm id}^2}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id}- \omega_{\rm ph}^2 \right) } .
\end{split}
\end{align}
\end{subequations}

\subsection*{DFT calculations of IR-phonons in $\rm Ta_2 Ni Se_5$}
The phonon spectrum at $\Gamma$ has been computed using the Density Functional Perturbation Theory routines of the VASP code \cite{Vasp1,Vasp2,Vasp3,Vasp4} using the vdW-opt88 functional \cite{Klime2011,Klimes2010} on a 48x4x6 mesh. The Born effective charges have been computed on a similar k-mesh, the IR-activity has been computed using phonopy \cite{phonopy} and the phonon spectroscopy code by Skelton et. al \cite{IR_activity}. The electronic bandstructure plots have been obtained using a 24x4x6 k-mesh and a 320 eV cutoff using the standard PBE functional.

\subsubsection{Phonons}
For the monoclinic geometry there are 21 symmetry allowed IR-active modes, 11 B$_\text{u}$ and 10 A$_\text{u}$ modes.Of these modes, only four are in a frequency range near 4.5THz, listed in table~\ref{tab: IR_phonons}.

\begin{table}[]
    \centering
    \begin{tabular}{c|c|c}
    Phonon & Symmetry & Frequency (THz) \\ \hline
    21  &  B$_\text{u}$  &  4.238  \\
    22  &  A$_\text{u}$  &  4.314  \\
    25  &  A$_\text{u}$  &  4.634  \\
    26  &  B$_\text{u}$  &  4.699  \\
    \end{tabular}
    \caption{IR-active phonon modes with a frequency close to 4.5 THz. }
    \label{tab: IR_phonons}
\end{table}

In the main text we identify phonon 26 which has a frequency of 4.699 THz as the relevant for phonon for parametric amplification. Therefore, we have investigated how this phonon evolves across the structural phase transition. To do this we have linearly interpolated between the monoclinic and high temperature orthorhombic atomic and lattice geometry using an interpolation parameter d.

\begin{equation}
\begin{split}
    \vec{v}_\text{i}(d) &= \vec{v}_\text{o} + d \cdot \vec{v}_\text{t} \\
    \vec{v}_\text{lat,i}(d) &= \vec{v}_\text{lat,o} + d \cdot \vec{v}_\text{lat,t} \quad , \\
\end{split}
\end{equation}
with the transition vector for the atomic configuration and the lattice defined as
\begin{equation}
\begin{split}
\vec{v}_\text{t} &= \vec{v}_{\text{m}} - \vec{v}_{\text{o}} \\
\vec{v}_\text{lat,t} &= \vec{v}_\text{lat,m} - \vec{v}_\text{lat,o} \quad .\\
\end{split}
\end{equation}
The subscripts m denotes the monoclinic geometry and o the orthorhombic geometry.
Because the phonon eigenmodes form an orthonormal basis for the movement of all atoms we can trace how the monoclinic 4.7 THz phonon hybridizes along the phase transition into the orthorhombic phonons, by calculating the overlap of the monoclinic phonon eigenvector with the phonon eigenvectors along the transition. The result is shown in figure \ref{fig: IR_phonon_overlap}. It shows that the 4.7 THz phonon predominantly disperses into two orthorhombic phonons at 4.38 THz and 4.56 THz.
\begin{figure}
    \centering
    \includegraphics[scale=0.42]{fatbands_IRmode26.pdf}
    \caption{Overlap of the 4.7 THz monoclinic phonon with the phonons calculated along the orthorhombic to monoclinic phase transition. The dot size is proportional to the overlap. The figure shows that the 4.7 THz phonon predominantly hybridizes into two orthorhombic phonons at frequencies of 4.38 THz and 4.56 THz. Note that the orthorhombic structure exhibits two phonon instabilities with complex frequency which are displayed as negative energies here.}
    \label{fig: IR_phonon_overlap}
\end{figure}

\subsubsection{Frozen phonon calculations}
To estimate the electron phonon coupling of the four IR-active phonon modes displayed in table \ref{tab: IR_phonons}, we show frozen phonon band structure calculations for these modes. As discussed in the main text, only the strongly amplified IR mode should display a considerable coupling to the electronic bands.
Upon displacing the atomic positions along the eigendisplacement of the four phonon modes and recalculating the electronic bandstructure, only mode 26 shows a strong modulation of the conduction bands, indicating a strong coupling of this mode to the electronic degrees of freedom upon excitation (see Fig. \ref{fig: IR_phonons}).

\begin{figure}
    \centering
    \includegraphics[scale=0.42]{PBE_bands.pdf}
    \caption{Frozen phonon bandstructure of all possible IR active active phonons in the monoclinic phase between 4.2 THz and 4.7 THz using the PBE functional. Phonon 26 shows she strongest coupling to the bandstructure upon displacement along its eigenmode. Most importantly, the second and third conduction band are modulated by this mode and are displaced in opposite directions.}
    \label{fig: IR_phonons}
\end{figure}

To quantify the strength of the electron-phonon coupling $\lambda_i$,  we define

\begin{equation}\label{EQ: e-ph coupling}
    \lambda_i = \sum_k  w_k | \,\, (E_{i,k})\bigr\rvert_{\langle Q \rangle = 0} - (E_{i,k})\bigr\rvert_{\langle Q \rangle \neq 0} \,\, | ,
\end{equation}
as the difference of the Kohn-Sham energies at the atomic equilibrium position $\langle Q \rangle = 0$ and the frozen phonon displaced position $\langle Q \rangle \neq 0$. k runs over all k-Points along the chosen Z-G path near the bandedge and $w_k$ is the symmetry weight of the k-Point. Computing this coupling for the band structures shown in figure \ref{fig: IR_phonons} one obtains for the first four conduction bands the values presented in table \ref{tab: IR_couplings}. It shows that the strongest coupling of the conduction bands is to phonon mode 26. For this mode the second and third conduction band are shifted into opposite directions along the Z-$\Gamma$ near the band edge which is consistent with the above proposed two band model. Therefore, we identify this mode as the relevant mode for the parametric amplification process.

\begin{table*}[]
    \centering
    \begin{tabular}{c|c|c|c|c}
         & Phonon21 (meV) & Phonon22 (meV) & Phonon25 (meV) & Phonon26 (meV) \\ \hline
           &         &         &        &         \\
cond. Band 1 &  1.197 &       2.762 &       3.321 &       8.171 \\
cond. Band 2 &  1.568 &       2.378 &       1.324 &      24.953 \\
cond. Band 3 &  5.363 &       3.391 &       6.473 &      34.436 \\
cond. Band 4 &  1.329 &       5.657 &      17.629 &      11.814 \\
    \end{tabular}
    \caption{Coupling as defined by equation \eqref{EQ: e-ph coupling} for the bandstructure presented in figure \ref{fig: IR_phonons} . Only states near the bandedge along Z-$\Gamma$ have been considered for the summation in equation \ref{EQ: e-ph coupling}. One can see that conduction bands 2 and 3 shift in opposite directions upon displacing along phonon mode 26 (see table \ref{tab: IR_phonons}). This is consistent with the parametric driving mechanism presented in the main text.}
    \label{tab: IR_couplings}
\end{table*}

\subsubsection{Optical transitions}
\begin{figure}
    \centering
    \includegraphics[scale=0.35]{PBE_fatbands_equ_05ev_allDir.pdf}
    \caption{Calculation of the dipole matrix elements for all three cartesian directions. We filtered only transitions in the experimentally relevant energy range between 0.33 eV and 0.9 eV. The size of the dots shows the size of the dipole matrix element at a given k-Point and band.}
    \label{fig: opticalTransitions}
\end{figure}
In the experiment, parametric amplification has been pumped with a 0.5 eV laser. In the above section we have identified a phonon mode which couples strongly to the second and third conduction band and might be responsible for the parametric amplification process. To fully justify this claim we need to show that using the 0.5 eV pump electrons are indeed excited into these strongly coupled conduction band states. Here we have computed the optical transition matrix elements and plotted them as dots on the bandstructure, where the dot size for the valence and conduction bands corresponds to
\begin{equation}
\begin{split}
    M_{\lambda, v,k} &= \sum_c M_{\lambda,v,c,k} \\
    M_{\lambda, c,k} &= \sum_v M_{\lambda,v,c,k} \\
\end{split}
\end{equation}
with $M_{\lambda,v,c,k} = |\langle \phi_{v,k} | (-i) \frac{d}{d k_\lambda} | \phi_{c,k} \rangle|^2 $.
$k$ labels the k-point, $v$ the valence, $c$ the conduction band involved in the optical transition and $\lambda$ the Cartesian directions. \newline
The result for all three Cartesian directions is shown in figure \ref{fig: opticalTransitions}. We have have filtered only the experimentally relevant transitions between 0.33 eV and 0.9 eV transition energy. Such a large energy window should compensate for any underestimation of the electronic gap given by the PBE functional. One sees that optical transitions in x-directions are strongest for the first 3 conduction band states. Thus, we can conclude that upon pumping with 0.5 eV there are many electrons excited into the electronic states that couple strongly to phonon 26 and can lead to the parametric amplification process.

\label{app:DFT}


\section{introduction}
Black hole is one of the most prominent predictions of general relativity. With the development of science and technology, the gravitational wave generated by the collision of two black holes was detected~\cite{Abbott2016,Cervantes2016,Krolak2021} and the shadow of black holes was revealed by Event Horizon Telescope collaboration~\cite{Akiyama2019}, which directly confirmed the existence of black holes, greatly stimulated people's interest in the physics of black holes. However, spacetime singularity is one of the most intriguing problems. According to Penrose and Hawking's singularity theorems, every black hole inevitably contains a singularity within the framework of Einstein's general relativity as long as some certain conditions are satisfied~\cite{Hawking1970}. At the singularity of spacetime, physical quantities and geometry diverge, and the laws of physics lose their power completely.
To avoid the singularity problem, Bardeen~\cite{Bardeen1968} first proposed a regular black hole model named the Badeen black hole, which was later proved to be the exact solution of the Einstein's field equations coupling with nonlinear electrodynamics~\cite{AyonBeato2000}. Subsequently, a number of other solutions to regular black holes such as Hayward black hole~\cite{Hayward2006}, Ay\'{o}n-Beato and Garc\'{\i}a black hole~\cite{AyonBeato1998} and  Simpson-Visser black holes~\cite{Simpson2019} have been proposed, and the rotating solutions of these black holes has been obtained in~\cite{Toshmatov2014,Abdujabbarov2016,Mazza2021}.

Up to now, Bardeen black hole has attracted extensive attention. Bret\'{o}n et al studied the stability of black holes with nonlinear electromagnetic fields~\cite{Breton2015} and Quasinormal modes was investigated by Fernando and Correa~\cite{Fernando2012}. Gravitational lensing by black holes with nonlinear electromagnetic fields was investigated in~\cite{Eiroa2011,Ghaffarnejad2016}. By analyzing the behavior of the effective potential, Zhou et al investigated the time-like and null geodesic structures~\cite{Zhou2012}. Abdujabbarov et al studied the Shadow of rotating regular black holes, and found that magnetic charge decreases the size of the shadow~\cite{Abdujabbarov2016}. Gao et al studied the bound orbits around the Badeen black hole~\cite{Gao2020} and greybody factor are studied in~\cite{Sharif2019,AmaTulMughani2021,Sharif2022}. The thermodynamic properties of the Bardeen black holes are discussed in~\cite{Man2014}. Some properties of Bardeen black holes in other modified theories have also been investigated~\cite{Mehdipour2016,Ghaderi2018,Saleh2018,Saleh20182,Pourhassan2019,SinghDV2020,SinghBK2020,Kumar2022}.

The motion of bodies in the vicinity of black holes is an important problem of relativistic astrophysics.
To solve the motion of the bodies in the strong gravitational fields, the post-Newtonian (PN) approximations have been extensively applied, and various analytical solutions have been obtained, including the first and higher PN effects of the mass~\cite{Brumberg1972,Soffel1987,Soffel1989,KlionerKopeikin1994,KopeikinEfroimskyKaplan2012,DamourSchafer1988,SchaferWex1993,MGS2004,Boetzel2017,Cho2018,SoffelHan2019,YangLin2020a,YangLin2020b}, the leading order~(1.5PN) spin effects~\cite{Wex1995,GPV1998a,GPV1998b,GPV1998c,KonigsdorfferGopakumar2005,KMG2005,KonigsdorfferGopakumar2006b,GopakumarSchafer2011,BMFB2013,GergelyKeresztes2015,Mikoczi2017} on the general motion, as well as the leading order~(2PN) mass quadrupole effects on the circular motion~\cite{Poisson1998}.

In our previous work~\cite{YangLin2019,YangJiangLin2022}, we have derived the quasi-Keplerian motion for the neutral test particle in the Reissner-Nordstr\"{o}m spacetime, and that for the charged test particle. We obtained the effects of the spin-induced quadrupole on the equatorial motion of a test particle in the Kerr spacetime~\cite{YangLin2020c} and the next-to-leading order spin-orbit effects on the test particle's generally inclined motion~\cite{YangLin2021}. In this work, we will derive the quasi-Keplerian motion of a test particle in the gravitational field of regular Bardeen black hole under the standard spherical coordinates. The magnetic charge effects of the Bardeen black hole on the quasi-Keplerian motion including the orbital perihelion precession are displayed explicitly.

The rest of this paper is organized as follows. Section \ref{sec:2nd} introduces the standard metric of Bardeen black hole in the 2PN approximation, the corresponding Lagrangian, as well as the orbital energy and angular momentum. In Section \ref{sec:3rd} we present the detail derivation on the 2PN solution for the quasi-Keplerian motion. Summary is given in section \ref{sec:4th}.

\section{The 2PN Lagrangian, energy, angular momentum in Bardeen spacetime}\label{sec:2nd}

Assume the regular Bardeen black hole is located at the coordinate origin. In the standard spherical coordinates, the exact metric of Bardeen black hole can be
expressed as
\begin{eqnarray}
&& ds^{2}=-f(r)dt^{2}+f(r)^{-1}dr^2+r^2d\theta^2+r^2\sin \theta^2 d\phi^2,  \label{HarmonicKerr}
\end{eqnarray}
where the lapse function $f(r)$ is given by
\begin{eqnarray}
&& f(r)=1-\frac{2mr^2}{(r^2+g^2)^{3/2}}~,
\end{eqnarray}
here $m$ denotes the standard gravitational mass parameter, while $g$ corresponds to the magnetic charge parameter measured in units of $m$. When $g=0$, this metric behaves as a Schwarzschild metric. When $g^2\leq\frac{16}{27}m^2$, a no-singularity black hole is described. In particular, when $g^2=\frac{16}{27}m^2$, single horizon solution is allowed, i.e., regular extremal black hole as shown in Fig~\ref{Fig1}. The gravitational constant and the speed of light in vacuum are set as 1.

\begin{figure}[h!]
\centering
\begin{subfigure}[t]{0.575\textwidth}
\centering
\includegraphics[width=9cm,height = 5.9cm]{blackhole.eps}
\end{subfigure}
\caption{ Horizons of the Bardeen black hole with different values of $g$.}~\label{Fig1}
\end{figure}

According to the theory of the PN approximation, this metric can be expanded into the power of $\frac{m}{r}$ with $r\!=\!|\bm{x}|$. Keeping all terms to the 2PN order for the dynamics of the massive particle, we have
\begin{eqnarray}
&& g_{00}= -1 +\frac{2m}{r}-\frac{3 g^2m}{r^3}~,\label{eq:metric-1nd}\\
&& g_{0i}= 0~,\label{eq:metric-2nd}\\
&& g_{ij}=\delta_{ij}+\Big(\frac{2m}{r}+\frac{4m^2}{r^2}\Big)\frac{x^ix^j}{r^2} ,\label{eq:metric-3nd}
\end{eqnarray}
where $i$ and $j$ running from 1 to 3, $\epsilon_{ijk}$ is Levi-Civita symbol~\cite{Weinberg1972}.

With the 2PN metric Eqs.\,(\ref{eq:metric-1nd})-(\ref{eq:metric-3nd}), the Lagrangian of the test particle can be written as
\begin{eqnarray}
&& {\rm L} = \frac{1}{2}\bm{v}^{2}+\frac{m}{r}+\frac{1}{8}\bm{v}^{4}
+\frac{1}{2}\frac{m}{r}\bm{v}^{2}+\frac{1}{2}\frac{m^2}{r^2}+\frac{m(\bm{v}\!\cdot\!\bm{x})^2}{r^3}+\frac{1}{2}\frac{m^3}{r^3}\Big(1-\frac{3g^2}{m^2}\Big)\nn\\
&&~~~~~~+\frac{1}{16}\bm{v}^{6}+\frac{3}{4}\frac{m^2}{r^2}\bm{v}^{2}+\frac{3}{8}\frac{m}{r}\bm{v}^{4}
+\frac{3m^2(\bm{v}\!\cdot\!\bm{x})^2}{r^4}+\frac{m(\bm{v}\!\cdot\!\bm{x})^2\bm{v}^{2}}{2r^3}~,\label{eq:Lagrangian}
\end{eqnarray}
where $\bm{v}$ denotes the the test particle's velocity.
According to this Lagrangian, the energy $\mathcal{E}$ and the angular momentum $\mathcal{J}$ of the equatorial motion can be obtained as
\begin{eqnarray}
&& \mathcal {E}
=\frac{1}{2}\bm{v}^{2}-\frac{m}{r}+\frac{3}{8}\bm{v}^{4}+\frac{1}{2}\frac{m}{r}\bm{v}^{2}-\frac{1}{2}\frac{m^2}{r^2}+\frac{m(\bm{v}\!\cdot\!\bm{x})^2}{r^3}
-\frac{1}{2}\frac{m^3}{r^3}\Big(1-\frac{3g^2}{m^2}\Big)\!\nn\\
&& ~~~~~~+\frac{5}{16}\bm{v}^{6}+
\frac{3}{4}\frac{m^2}{r^2}\bm{v}^{2}+\frac{9}{8}\frac{m}{r}\bm{v}^{4}
+\frac{3m^2(\bm{v}\!\cdot\!\bm{x})^2}{r^4}+\frac{3}{2}\frac{m(\bm{v}\!\cdot\!\bm{x})^2\bm{v}^{2}}{r^3}~,\label{eq:total-energy}\\
&& \mathcal{J}
= |\bm{x}\!\times\!\bm{v}|\Big[1+\frac{1}{2}\bm{v}^{2}
+\frac{m}{r}+\frac{3}{8}\bm{v}^{4}+\frac{3}{2}\frac{m^2}{r^2}+\frac{3}{2}\frac{m}{r}\bm{v}^{2}+\frac{m(\bm{v}\!\cdot\!\bm{x})^2}{r^3}\Big]~.
\label{eq:total-angular momentum}
\end{eqnarray}

\section{The quasi-Keplerian motion in the 2PN approximation}\label{sec:3rd}

We deduce the analytical solution for the quasi-Keplerian motion according to the same procedure given by Brumberg~\cite{Brumberg1972}.

The trajectory of the test particle in the equatorial plane of the Bardeen black hole can be expressed as
\begin{equation}
\bm{x} = r(\cos\phi\,\bm{e}_{x}+\sin\phi\,\bm{e}_{y})~,\label{eq:x1}
\end{equation}
where $\phi$ is the azimuthal angle. $\bm{e}_{x}$ and $\bm{e}_{y}$ are the unit vectors of the $x$-axis and $y$-axis.

We can re-write the orbital energy and angular momentum in Eqs.~(\ref{eq:total-energy})-(\ref{eq:total-angular momentum}) as
\begin{eqnarray}
&& \hskip -0cm \mathcal {E}=
\frac{1}{2}(\dot{r}^2\!+\!r^2\dot{\phi}^2)\!-\!\frac{m}{r}\!+\!\frac{3}{8}(\dot{r}^2\!+\!r^2\dot{\phi}^2)^2\!+\!\frac{1}{2}\frac{m}{r}(\dot{r}^2\!+\!r^2\dot{\phi}^2)
\!-\!\frac{m^2}{2r^2}\!+\!\frac{m}{r}\dot{r}^2\!-\!\frac{m^3}{2r^3}\Big(1\!-\!\frac{3g^2}{m^2}\Big)\nn\\
&& \hskip 1cm +\frac{5}{16}(\dot{r}^2\!+\!r^2\dot{\phi}^2)^3\!+\!\frac{3}{4}\frac{m^2}{r^2}(\dot{r}^2\!+\!r^2\dot{\phi}^2)\!
+\!\frac{9}{8}\frac{m}{r}(\dot{r}^2\!+\!r^2\dot{\phi}^2)^2
\!+\!\frac{3m^2}{r^2}\dot{r}^2\!+\!\frac{3}{2}\frac{m}{r}\dot{r}^2(\dot{r}^2\!+\!r^2\dot{\phi}^2) ,\label{eq:Brumberg-total-energy}\\
&& \hskip -0cm\mathcal{J}=r^2\dot{\phi}
\Big[1\!+\!\frac{1}{2}(\dot{r}^2\!+\!r^2\dot{\phi}^2)\!+\!\frac{m}{r}\!+\!\frac{3}{8}(\dot{r}^2\!+\!r^2\dot{\phi}^2)^2\!+\!
\frac{3}{2}\frac{m}{r}(\dot{r}^2\!+\!r^2\dot{\phi}^2)\!+\!\frac{3}{2}\frac{m^2}{r^2}\!+\!\frac{m}{r}\dot{r}^2\Big]~,\label{eq:Brumberg-total-angular momentum}
\end{eqnarray}
where the dot means the derivative with respect to the time. From these two expressions we can obtain
\begin{equation}
r^4\dot{\phi}^2=\mathcal{J}^2\Big(1\!-\!2\,\mathcal {E}\!-\!\frac{4m}{r}\!+\!3\,\mathcal {E}^2\!+\!\frac{4m^2}{r^2}\!+\!8\,\mathcal {E}\frac{m}{r}\Big)~,\label{eq:r4dottheta}
\end{equation}
and
\begin{equation}
 \dot{r}^2=A+\frac{B}{r}+\frac{C}{r^2}+\frac{D}{r^3}+\frac{E}{r^4}~,\label{eq:rdot20}
\end{equation}
with
\begin{eqnarray}
&& A=2\mathcal {E}\Big(1\!-\!\frac{3}{2}\,\mathcal {E}\!+\!2\,\mathcal {E}^2\Big)~,\\
&& B=2m\Big(1\!-\!6\,\mathcal {E}\!+\!9\,\mathcal {E}^2\Big)~,\\
&& C=-\mathcal{J}^2\Big(1\!-\!2\,\mathcal{E}\!+\!\frac{8m^2}{\mathcal{J}^2}\!-\!24\frac{m^2\mathcal {E}}{\mathcal{J}^2}\!+\!3\,\mathcal {E}^2\Big)~,\\
&& D=6m \mathcal{J}^2\Big[1\!-\!2\,\mathcal{E}\!+\!\frac{4}{3}\frac{m^2}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\\
&& E=-12m^2 \mathcal{J}^2~.
\end{eqnarray}
Throughout this paper, the signs ``$\pm$" denote the anti-clockwise motion and the clockwise motion, respectively.

utilizing the relation
\begin{equation}
 \dot{r}^2=\Big[\frac{d(1/r)}{d\phi}\Big]^2(r^4\dot{\phi}^2)~,\label{eq:rdot2-relation}
\end{equation}
and plugging Eqs.~(\ref{eq:r4dottheta})-(\ref{eq:rdot20})  into (\ref{eq:rdot2-relation}), we can
express the radial equation in the form
\begin{eqnarray}
&& \Big[\frac{d(1/r)}{d\phi}\Big]^2=A'+\frac{B'}{r}+\frac{C'}{r^2}+\frac{D'}{r^3}+\frac{E'}{r^4}~,\label{eq:the-radial-equation1}
\end{eqnarray}
with
\begin{eqnarray}
&& A'=\frac{2\mathcal {E}}{\mathcal{J}^2}\Big(1\!+\!\frac{1}{2}\,\mathcal {E}\Big)~,\\
&& B'=\frac{2m}{\mathcal{J}^2}~,\\
&& C'=-1~,\\
&& D'=2m\Big(1\!-\!\frac{3g^2}{2\mathcal{J}^2}\Big)~,\\
&& E'=0~.
\end{eqnarray}
Since the right hand side of Eq.~(\ref{eq:the-radial-equation1}) is a fourth-order polynomial in $r^{-1}$, we can further re-write it as
\begin{eqnarray}
&& \Big[\frac{d(1/r)}{d\phi}\Big]^2=\Big[\frac{1}{r}\!-\!\frac{1}{a_{r}
(1+e_{r})}\Big]\Big[\frac{1}{a_{r}(1-e_{r})}\!-\!\frac{1}{r}\Big]\Big(C_{1}+\frac{C_{2}}{r}+\frac{C_{3}}{r^2}\Big)~.\label{eq:the-radial-equation2}
\end{eqnarray}
Comparing the coefficients between Eq.~(\ref{eq:the-radial-equation1}) and Eq.~(\ref{eq:the-radial-equation2}), we have
\begin{eqnarray}
&& a_{r}=\frac{m}{-2\mathcal{E}}\Big[1\!+\!\frac{3}{2}\mathcal{E}\!+\!\frac{1}{4}\mathcal{E}^2\!+\!8\frac{m^2\mathcal{E}}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\label{eq:aPN_EJ}\\
&& e_r^2=1\!+\!\frac{2\mathcal{E}\mathcal{J}^2}{m^2}\!-\!\mathcal{E}\Big(8\!+\!7\frac{\mathcal{E}\mathcal{J}^2}{m^2}\Big)
\!+\!\,\mathcal{E}^2\Big[16\frac{\mathcal{E}\mathcal{J}^2}{m^2}\!-\!20\Big(1\!-\!\frac{3g^2}{5m^2}\Big)\!-\!32\frac{m^2}{\mathcal{E}\mathcal{J}^2}
\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\label{eq:ePN_EJ}\\
&& C_{1}=1\!-\!\frac{4m^2}{\mathcal{J}^2}\!-\!8\frac{m^2\mathcal{E}}{\mathcal{J}^2}\!-\!\frac{16m^4}{\mathcal{J}^4}\Big(1\!-\!\frac{3g^2}{8m^2}\Big)~,\label{eq:C1}\\
&& C_{2}=-2m\Big(1\!-\!\frac{3g^2}{2\mathcal{J}^2}\Big)~,\label{eq:C2}\\
&& C_{3}=0~.\label{eq:C3}
\end{eqnarray}
It can be seen from Eq.~(\ref{eq:the-radial-equation2}) that $r_{\pm}=a_{r}(1\pm e_{r})$ represent the maximal and minimal values for $r$. Hence, $a_{r}$ and $e_{r}$ can be regarded as the semi-major axis and the eccentricity of the quasi-Keplerian orbit.

The solution of Eq.~(\ref{eq:the-radial-equation2}) can be written as:
\begin{eqnarray}
&& r=\frac{a_{r}(1-e_{r}^2)}{1+e_{r}\cos{f}}~,\label{eq:rBrumberg}
\end{eqnarray}
with $f$ being the true anomaly for the quasi-Keplerian orbit and satisfying
\begin{eqnarray}
&& \Big(\frac{df}{d\phi}\Big)^2=C_{1}+\frac{C_{2}}{r}+\frac{C_{3}}{r^2}~.\label{eq:f_phi2}
\end{eqnarray}

Substituting Eqs.~(\ref{eq:C1})-(\ref{eq:rBrumberg}) into Eq.~(\ref{eq:f_phi2}), we have
\begin{eqnarray}
&& \frac{df}{d\phi}=F\Big\{1\!-\! \frac{m^2}{\mathcal{J}^2}\Big[1\!+\!2\mathcal{E}\!+\!10\frac{m^2}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{20m^2}\Big)\Big]e_{r} \cos f
\!-\!\frac{m^4}{4\mathcal{J}^4}e_{r}^2 \cos 2f\Big\}~,\label{eq:dfdphi}
\end{eqnarray}
with
\begin{eqnarray}
F=\Big[1\!-\!\frac{3m^2}{\mathcal{J}^2}\!-\!\frac{13}{2}\frac{m^2\mathcal{E}}{\mathcal{J}^2}\!-\!\frac{67}{4}\frac{m^4}{\mathcal{J}^4}\Big(1\!-\!\frac{18g^2}{67m^2}\Big)
\Big]~.\label{eq:F}
\end{eqnarray}

Making integration of Eq.~(\ref{eq:dfdphi}), we can obtain
\begin{equation}
\phi \Big(\frac{2\pi}{\Phi}\Big)=f+\frac{m^2}{\mathcal{J}^2}\Big[1\!+\!2\mathcal{E}\!+\!10\frac{m^2}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{20m^2}\Big)\Big]e_{r} \sin f\!+\!\frac{3}{8}\frac{m^4}{\mathcal{J}^4}e_r^2\sin 2f~,\label{eq:phi_f}
\end{equation}
with
\begin{equation}
 \Phi=2\pi\Big[1\!+\!3\frac{m^2}{\mathcal{J}^2}\!+\!\frac{15}{2}\frac{m^2\mathcal{E}}{\mathcal{J}^2}\!+\!\frac{105}{4}\frac{m^4}{\mathcal{J}^4}
\Big(1\!-\!\frac{6g^2}{35m^2}\Big)\Big]~.\label{eq:Phi}
\end{equation}

Finally, we derive the time dependence of the quasi-Keplerian motion.
Combining Eqs.~(\ref{eq:r4dottheta}) and (\ref{eq:dfdphi})-(\ref{eq:F}), we have
\begin{eqnarray}
&& \hskip -0cm r^2\dot{f}=\mathcal{J}\Big\{1\!-\!\mathcal {E}\!-\!2\frac{m}{r}\!-\!3\frac{m^2}{\mathcal{J}^2}
\!+\!\mathcal {E}^2\!+\!2\mathcal {E}\frac{m}{r}\!-\!\frac{7}{2}\frac{m^2\mathcal{E}}{\mathcal{J}^2}\!
+\!6\frac{m}{r}\frac{m^2}{\mathcal{J}^2}\!-\!\frac{m^4}{\mathcal{J}^4}
\Big(\frac{67}{4}\!-\!\frac{9g^2}{2m^2}\Big)\nn\\
&& \hskip 1.5cm -\frac{m^2}{\mathcal{J}^2}\Big[1\!+\!\mathcal {E}\!-\!2\frac{m}{r}\!+\!\frac{m^2}{\mathcal{J}^2}\Big(7\!-\!\frac{3g^2}{2m^2}\Big)\Big]e_{r} \cos f
\!-\!\frac{1}{4}\frac{M^4}{\mathcal{J}^4} e_{r}^2\cos 2f\Big\}.\label{eq:dfdt}
\end{eqnarray}

Introducing the post-Newtonian eccentric anomaly $u$ by the relations
\begin{equation}
\hskip -0cm \sin f\!=\!\frac{(1\!-\!e_{r}^2)^{\frac{1}{2}}\sin u}{1\!-\!e_{r}\cos u};\,\,\cos f\!=\!\frac{\cos u\!-\!e_{r}}{1\!-\!e_{r}\cos u};\,\,
 f\!=\!2 \!\arctan\!\Big(\sqrt{\frac{1\!+\!e_{r}}{1\!-\!e_{r}}}\tan \frac{u}{2}\Big),
\label{eq:u-relation}
\end{equation}
we have
\begin{eqnarray}
\frac{df}{dt}=\frac{(1-e_{r}^2)^{1/2}}{1-e_{r}\cos u}\frac{du}{dt}~,\label{eq:dfdu}
\end{eqnarray}
and we can formulate the orbit given in Eq.~(\ref{eq:rBrumberg}) in terms of $u$ as
\begin{equation}
 r=a_{r}(1-e_{r} \cos u)~.\label{eq:r}
\end{equation}

Integrating Eq.~(\ref{eq:dfdt}) and making use of Eqs.~(\ref{eq:u-relation})-(\ref{eq:r}), we can achieve the final piece of the 2PN closed-form solution for the equatorial motion in Bardeen spacetime.
\begin{equation}
 t\Big(\frac{2\pi}{{\rm T}_u}\Big)  = u-e_t \sin{u}+ \frac{30m \,\mathcal{E}^2 }{\sqrt{-2\,\mathcal{E}\mathcal{J}^2}} (f-u)~,\label{eq:nt}
\end{equation}
with ${\rm T}_{\!u}$ being the period for the eccentric anomaly $u$ of the quasi-Keplerian motion
\begin{equation}
 {\rm T}_u=\frac{2\pi m}{(-2\mathcal{E})^{\frac{3}{2}}}\Big[1\!-\!\frac{15}{4}\mathcal{E}\!-\!\frac{105}{32}\mathcal{E}^2\!+\!\frac{30m\mathcal{E}^2}{\sqrt{-2\mathcal{E}\mathcal{J}^2}}\Big]~,
\end{equation}
and $e_t$ being the time eccentricity
\begin{equation}
 e_t=e_r\Big[1\!+\!6\,\mathcal{E}\!+\!27\,\mathcal{E}^2\!-\!\frac{30m\mathcal{E}^2}{\sqrt{-2\mathcal{E}\mathcal{J}^2}}\!+\!8\frac{m^2\mathcal{E}}{\mathcal{J}^2}
\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~.
\end{equation}

In the literature, one usually replace the true anomaly $f$ in the formula of quasi-Keplerian equation with another true anomaly $\upsilon$, which requires that $\sin \upsilon$ contribution vanishes at each PN order in $\phi (\frac{2\pi}{\Phi})$~\cite{MGS2004,KonigsdorfferGopakumar2005,THS2010}.

Following the same method given in Ref~\cite{THS2010}, we set
\begin{equation}
 \upsilon=2 \arctan \Big(\sqrt{\frac{1+e_{\phi}}{1-e_{\phi}}}\tan \frac{u}{2}\Big)~,\label{eq:upsilon_u}
\end{equation}
with
\begin{equation}
 e_\phi=e_r(1+\epsilon\, c_1+\epsilon^2\, c_2)~,
\end{equation}
differing from the radial eccentricity $e_r$ by some 1PN and 2PN level corrections $c_1$ and $c_2$. Here $\epsilon$ only denotes the PN order and does not have any value. Eliminating $u$ in Eq.~(\ref{eq:u-relation}) with the help of Eq.~(\ref{eq:upsilon_u}), we have~\cite{THS2010}
\begin{equation}
\hskip -0.05cm f\!=\!\upsilon+\epsilon\, c_1\frac{e_r}{e_r^2\!-\!1}\sin\upsilon
+\epsilon^2\,\Big[\Big(c_2\!-\!c_1^2\frac{e_r^2}{e_r^2\!-\!1}\Big)\frac{e_r}{e_r^2
\!-\!1}\sin \upsilon+\frac{c_1^2}{4}\frac{e_r^2}{(e_r^2\!-\!1)^2}\sin 2\upsilon\Big].
\end{equation}
Substituting this result into Eq.~(\ref{eq:phi_f}) and requiring the $\sin \upsilon$ term to vanish in $\phi (\frac{2\pi}{\Phi})$, we can obtain
\begin{eqnarray}
&& c_1=-2\mathcal{E}~,\\
&& c_2=-\!18\mathcal{E}\frac{m^2}{\mathcal{J}^2}\Big(1\!-\!\frac{g^2}{6m^2}\Big)~,
\end{eqnarray}
which lead to
\begin{equation}
 e_{\phi} = e_r \Big[1\!-\!2\mathcal{E}\!-\!18\frac{m^2\mathcal{E}}{\mathcal{J}^2}\Big(1\!-\!\frac{g^2}{6m^2}\Big)\Big]~,
\end{equation}
\begin{equation}
 \phi \Big(\frac{2\pi}{\Phi}\Big)= \upsilon \!+\!\frac{1}{8}\frac{m^4}{\mathcal{J}^4}e_r^2\sin 2\upsilon~.
\end{equation}

With the true anomaly $\upsilon$, we can re-express the time dependance of the quasi-Keplerian motion Eq.~(\ref{eq:nt}) in the form of
\begin{equation}
 t\Big(\frac{2\pi}{{\rm T}_u}\Big)  = u-e_t \sin{u}+\frac{30m \,\mathcal{E}^2 }{\sqrt{-2\,\mathcal{E}\mathcal{J}^2}} (\upsilon-u)~.
\end{equation}

Notice that all the formulas in this work are valid up to the 2PN accuracy.

\section{Summary}\label{sec:4th}
We start with the 2PN metric of the Bardeen black hole in the standard coordinates, calculate the corresponding orbital energy and angular momentum of the test particle in the equatorial plane, and then through an iterative method and function fitting method, derive the 2PN solution for the quasi-Keplerian motion in the Bardeen spacetime. We obtain two slightly different but equivalent formulations in the 2PN approximation.

The results are summarized as follows.

The first formulation can be expressed as
\begin{eqnarray}
&& \bm{x} = r(\cos\phi\,\bm{e}_{x}+\sin\phi\,\bm{e}_{y})~,\nn\\
&& r=a_{r}(1-e_{r} \cos u)~,\nn\\
&& \phi \Big(\frac{2\pi}{\Phi}\Big)=f+ N_0\sin f+N_1\sin 2f~,\nn\\
&& f=2 \arctan \Big(\sqrt{\frac{1+e_{r}}{1-e_{r}}}\tan \frac{u}{2}\Big)~,\nn\\
&&  t\Big(\frac{2\pi}{{\rm T}_u}\Big)  = u-e_t \sin{u}+ N_2 (f-u)~,\nn
\end{eqnarray}
and the second formulation can be expressed as
\begin{eqnarray}
&& \bm{x} = r(\cos\phi\,\bm{e}_{x}+\sin\phi\,\bm{e}_{y})~,\nn\\
&& r=a_{r}(1-e_{r} \cos u)~,\nn\\
&& \phi \Big(\frac{2\pi}{\Phi}\Big)= \upsilon + N_3 \sin 2\upsilon~,\nn\\
&& \upsilon=2 \arctan \Big(\sqrt{\frac{1+e_{\phi}}{1-e_{\phi}}}\tan \frac{u}{2}\Big)~,\nn\\
&&  t\Big(\frac{2\pi}{{\rm T}_u}\Big)  = u-e_t \sin{u}+ N_2 (\upsilon-u)~,\nn
\end{eqnarray}
where
\begin{eqnarray}
&& a_{r}=\frac{m}{-2\mathcal{E}}\Big[1\!+\!\frac{3}{2}\mathcal{E}\!+\!\frac{1}{4}\mathcal{E}^2\!+\!8\frac{m^2\mathcal{E}}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\nn\\
&& e_r^2=1\!+\!\frac{2\mathcal{E}\mathcal{J}^2}{m^2}\!-\!\mathcal{E}\Big(8\!+\!7\frac{\mathcal{E}\mathcal{J}^2}{m^2}\Big)
\!+\!\,\mathcal{E}^2\Big[16\frac{\mathcal{E}\mathcal{J}^2}{m^2}\!-\!20\Big(1\!-\!\frac{3g^2}{5m^2}\Big)\!-\!32\frac{m^2}{\mathcal{E}\mathcal{J}^2}
\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\nn\\
&& e_t=e_r\Big[1\!+\!6\,\mathcal{E}\!+\!27\,\mathcal{E}^2\!-\!\frac{30m\mathcal{E}^2}{\sqrt{-2\mathcal{E}\mathcal{J}^2}}\!+\!8\frac{m^2\mathcal{E}}{\mathcal{J}^2}
\Big(1\!-\!\frac{3g^2}{8m^2}\Big)\Big]~,\nn\\
&& e_{\phi} = e_r\Big[1\!-\!2\mathcal{E}\!-\!18\frac{m^2\mathcal{E}}{\mathcal{J}^2}\Big(1\!-\!\frac{g^2}{6m^2}\Big)\Big]~,\nn\\
&& \Phi=2\pi\Big[1\!+\!3\frac{m^2}{\mathcal{J}^2}\!+\!\frac{15}{2}\frac{m^2\mathcal{E}}{\mathcal{J}^2}\!+\!\frac{105}{4}\frac{m^4}{\mathcal{J}^4}
\Big(1\!-\!\frac{6g^2}{35m^2}\Big)\Big]~,\nn\\
&&  N_0= \frac{m^2}{\mathcal{J}^2}\Big[1\!+\!2\mathcal{E}\!+\!10\frac{m^2}{\mathcal{J}^2}\Big(1\!-\!\frac{3g^2}{20m^2}\Big)\Big]
\Big(1\!+\!\frac{2\mathcal{E}\mathcal{J}^2}{m^2}\Big)^{\!\frac{1}{2}}~,\nn\\
&&  N_1= \frac{3}{8}\frac{m^4}{\mathcal{J}^4}\Big(1+\frac{2\mathcal{E}\mathcal{J}^2}{m^2}\Big)~,\nn\\
&&  N_2 =\frac{30m \,\mathcal{E}^2 }{\sqrt{-2\,\mathcal{E}\mathcal{J}^2}}~,\nn\\
&&  N_3= \frac{1}{8}\frac{m^4}{\mathcal{J}^4}\Big(1+\frac{2\mathcal{E}\mathcal{J}^2}{m^2}\Big)~,\nn\\
&& {\rm T}_u=\frac{2\pi m}{(-2\mathcal{E})^{\frac{3}{2}}}\Big(1\!-\!\frac{15}{4}\mathcal{E}\!-\!\frac{105}{32}\mathcal{E}^2
+\frac{30m\mathcal{E}^2}{\sqrt{-2\mathcal{E}\mathcal{J}^2}}\Big)~.\nn
\end{eqnarray}
In the formulations, $a_{r}$, $e_{r}$ and $u$ denote the semi-major axis, the eccentricity, the eccentric anomaly of the quasi-Keplerian motion in the PN approximation, respectively. $f$ and $\upsilon$ are two slightly different definitions of the true anomaly. ${\rm T}_u$ denotes the orbital period.
Notice that the effects of the magnetic charge on the test particle's motion are characterized by the terms containing $g^2$ in the above formulas.

The test particle's orbital perihelion precession can be calculated by \begin{eqnarray}
\Delta \phi \!\equiv\! \Phi-2\pi=6\pi\frac{m^2}{\mathcal{J}^2}\!+\!15\frac{\pi m^2\mathcal{E}}{\mathcal{J}^2}\!+\!\frac{105}{2}\frac{\pi m^4}{\mathcal{J}^4}
\Big(1\!-\!\frac{6g^2}{35m^2}\Big)~.\label{deltaphi}
\end{eqnarray}
When the magnetic charge of the Bardeen black hole vanishes ($g=0$), Eq.~(\ref{deltaphi}) can be simplified to the 2PN Schwarzschild periastron precession~\cite{DamourSchafer1988,MGS2004}.

Finally, it can be observed that the test particle's orbital period is independent on the the magnetic charge up to the 2PN order.

\section*{ACKNOWLEDGEMENT}
This work was supported in part by the National Natural Science Foundation of China (Grant Nos. 11973025 and 11947404).