# Squeezed Coherent States of Motion for Ions Confined in Quadrupole and
Octupole Ion Traps

Bogdan M. Mihalcea  Natl. Inst. for Laser, Plasma and Radiation Physics
(INFLPR),
Atomiştilor Str. Nr. 409, 077125 Măgurele, Romania<EMAIL_ADDRESS>
###### Abstract

Quasiclassical dynamics of trapped ions is characterized by applying the time
dependent variational principle (TDVP) on coherent state orbits, in case of
quadrupole and octupole combined (Paul and Penning) and radiofrequency (RF)
traps. A dequantization algorithm is proposed, by which the classical Hamilton
(energy) function associated to the system results as the expectation value of
the quantum Hamiltonian on squeezed coherent states. We develop such method
and particularize the quantum Hamiltonian for a combined and for a RF trap,
with axial symmetry and a RF anharmonic electric potential. We also build the
classical Hamiltonian functions for the particular traps we considered, and
find the classical equations of motion.

###### keywords:

Trapped ion, squeezed coherent states, anharmonic electric potential, Time-
Dependent Variational Principle , octupole combined ion trap.

###### PACS:

02.20.-a, 03.65.-w, 37.10.Ty

††journal: Journal of LaTeX Templates

##  1 Introduction

Recent advances in quantum optics [1] enable trapping of single particles or
atoms [2, 3], while progress in quantum engineering techniques allows
preparing these particles in well-defined quantum states [4, 5, 6]. Quantum
engineering using ion traps offers the possibility to prepare stable quantum
states by precise control of the interactions between a quantum system
(trapped ions) and the environment [7, 8, 9], while investigation of
nonclassical states of spin systems coupled to a harmonic oscillator offers
the possibility to investigate fundamental quantum phenomena, such as the
mechanisms responsible for decoherence and the quantum-classical transition
[10, 11].

The quantum time-dependent harmonic oscillator has been intensively used to
describe the dynamics of many physical systems [12, 13, 14, 15]. Quantum
dynamics of harmonic oscillators is obtained by the use of the so-called
Peremolov’s generalized coherent states [16, 17] of the Lie algebra associated
to the SU(1,1) group [18]. Quantum dynamics in a 3D RF ion trap characterized
by a uniform magnetic field and a time-dependent quadrupole electric
potential, can be reduced to the solution of the time-dependent quantum
oscillator equations, obtained by separating the axial and radial motion from
the Schrödinger equation [19]. In addition, the coherent state approach leads
to quantum solutions that are explicitly constructed as functions of the
classical trajectories on the phase space [17]. The properties of a RF
(Paul)-trap with a super-imposed magnetic field (combined trap) are presented
in [20], and it was shown that the regions of stability are significantly
larger than those for a Paul trap for both positive and negative charged ions.
Ion dynamics in a radiofrequency (RF) octupole trap was described in [21, 22],
demonstrating confinement of ions in a RF anharmonic electric potential and
characterizing the stability of this nonlinear parametric oscillator.
Collective dynamics for systems of ions confined in quadrupole 3D traps with
cylindrical symmetry is characterized in [23]. Ion dynamics in a linear
combined trap has been investigated both theoretically and experimentally in
[24], and it was shown that the presence of a homogeneous magnetic field
superimposed over the applied DC and RF electric fields, leads to a set of
coupled Mathieu equations. Coherent states for a set of quadratic Hamiltonians
in the trap regime are constructed in Ref. [25], and then particularized to
the asymmetric Penning trap. A method of finding a set of generators that form
a closed Lie algebra, which creates a framework to characterize a general
quantum Hamiltonian is presented in [26]. Thus, the Lie algebra can be
extended to study the Hamiltonian of a bi-dimensional charged particle in
time-dependent electromagnetic fields, by exploiting the similarities between
the terms of these two Hamiltonians.

This paper characterizes the evolution of squeezed coherent states of motion
for ions confined in quadrupole and octupole combined and RF traps, using the
coherent state formalism [27, 16, 19] developed in [28, 18], and the time
dependent variational principle (TDVP) [29]. We study the bosonic realization
of the Lie algebra for the SU(1,1) group, and (generalized) coherent states in
the Fock space for a trapped ion [30]. The paper is organized as follows: In
Section 2 we propose a dequantization algorithm [31] that enables explicit
calculus of the quantum Hamilton function associated to an anharmonic
oscillator (ion) confined in a combined or RF trap, which describes an
algebraic model when the anharmonic part is a polynomial. Such model is linear
for 3D quadrupole ion traps that exhibit axial symmetry. Section 2 presents
the solutions of the Schrödinger equation and the quasienergy spectrum for the
model we suggest. The method suggested in [19, 32] is developed and
particularized in Section 3 for both combined (Paul and Penning) and RF traps,
considering a RF anharmonic electric potential. We build the classical
Hamilton functions for such particular traps and find the classical Hamilton
equations of motion for the anharmonic combined trap, which represents another
original result. We suggest that in the pseudopotential approximation case
(ideal RF trap), the points of minimum of the classical Hamiltonian describe
equilibrium configurations for trapped ions, of interest for implementing
quantum logic. The results are discussed in Section 4, as the model is
straightforward to extend for 2D ion traps.

##  2 Quantum dynamics of ions in axially symmetric quadrupole and octupole
ion traps

### Algebraic models for axially symmetric nonlinear quadrupole and octupole
traps

The dynamical group $Sp\left(2,{\mathbb{R}}\right)_{a}\otimes
Sp\left(2,{\mathbb{R}}\right)_{r}\otimes SO\left(2\right)$ associated to the
Hamilton function which describes the dynamics of an ion of mass $M$ and
electric charge $Q$, confined within a 3D RF ion trap that exhibits both
cylindrical and reflection symmetry, is the direct product between the axial
and radial symplectic groups, and the rotations group $SO\left(2\right)$
generated by the axial angular moment operator $L_{z}$ [28]. The Lie algebra
basis of $Sp\left(2,{\mathbb{R}}\right)_{j},\;j=a,r$, is spanned by the
generators $K_{0j},K_{1j}$ and $K_{2j}$. We introduce the infinitesimal
generators of the axial symplectic group $Sp\left(2,{\mathbb{R}}\right)_{a}$
defined as [28, 19]

$K_{0,1a}=\frac{M\omega_{a}}{4\hbar}\left[z^{2}\pm\frac{p_{z}^{2}}{M^{2}\omega_{a}^{2}}\right]\
,\;K_{2a}=\frac{i}{4\hbar}\left[2z\frac{\partial}{\partial z}+1\right]\ ,$ (1)

where $\omega_{a}/2\pi$ is the trap axial frequency. The commuting relations
are

$\left[K_{0a},K_{1a}\right]=iK_{2a}\ ,\;\left[K_{2a},K_{0a}\right]=iK_{1a}\
,\;\left[K_{1a},K_{2a}\right]=iK_{0a}\ .$ (2)

Using eqs. (1) we infer

$\left(K_{0a}+K_{1a}\right)\left(K_{0a}-K_{1a}\right)=iK_{2a}+\frac{1}{4\hbar^{2}}z^{2}p_{z}^{2}\
.$ (3)

The Casimir operator that determines the Bargmann indices $k$ is [28, 16, 19]

$C_{2a}=K_{0a}^{2}-K_{1a}^{2}-K_{2a}^{2}=-\frac{3}{16}I=k\left(k-1\right)I\ .$
(4)

where $I$ is the unit operator. The Bargmann indices that characterize axial
motion are $k_{a+}=\frac{1}{4}\ ,\;k_{a-}=\frac{3}{4}$. We use cylindrical
coordinates $\rho$ and $\theta$, with radial coordinates $x=\rho\cos\theta$
and $y=\rho\sin\theta$. The radial symplectic group generators for fixed
angular moment $l$ are [28, 19, 23]

$K_{0,1r}=\frac{1}{2\hbar\omega_{r}}\left[\frac{1}{2M}p_{\rho}^{2}\pm\frac{M\omega_{r}^{2}}{2}\rho^{2}+\frac{\hbar^{2}}{2M}\left(l^{2}-\frac{1}{4}\right)\frac{1}{\rho^{2}}\right]\
,$ (5a) $K_{2r}=\frac{i}{4}\left[2\rho\frac{\partial}{\partial\rho}+1\right]\
,$ (5b)

where $\omega_{r}/2\pi$ is the trap radial frequency and
$p_{\rho}=-i\hbar\frac{\partial}{\partial\rho}$. Because the angular momentum
operator commutes both with the generators of the symplectic groups and the
quantum Hamilton function $H$, we can restrict the study of the quantum system
to a subspace of the Hilbert space [33] with axial angular momentum
$L_{z}=\hbar lI$, where $l$ represents the orbital quantum number [32]. The
commutation relations are

$\left[K_{0r},K_{1r}\right]=iK_{2r}\ ,\;\left[K_{0r},K_{2r}\right]=-iK_{1r}\
,\;\left[K_{1r},K_{2r}\right]=iK_{0r}\ $ (6)

Using eq. (5) we infer

$K_{0r}^{2}-K_{1r}^{2}=\frac{1}{4}\left(-\rho^{2}\frac{\partial^{2}}{\partial\rho^{2}}+l^{2}-\frac{1}{4}\right)+iK_{2r}.$
(7)

The Casimir operator for the radial symplectic group states is [28, 16, 19,
32]

$C_{2r}=K_{0r}^{2}-K_{1r}^{2}-K_{2r}^{2}=\frac{l^{2}-1}{4}=k\left(k-1\right)\
,$ (8)

therefore

$k_{r12}=\frac{1\pm l}{2}\Rightarrow
k_{r}=\frac{l+1}{2},\;k_{r}>0\;\mbox{and}\;l\geq 0\ .$ (9)

We choose a potential with axial symmetry [28, 32]

$\Phi\left(\vec{r},t\right)=A(t)g\left(\rho,z\right)\ ,$ (10)

where $A\left(t\right)$ is a time periodic function, of period
$T=2\pi/\Omega$, where $\Omega$ is the frequency of the RF trapping voltage,
and $g$ is a function of $\rho^{2}$ and $z^{2}$

$g\left(\rho,z\right)=\sum_{k\geq 1}c_{k}H_{2k}\left(\rho,z\right)\ .$ (11)

The $H_{2k}$ polynomials are harmonic, of degree $k$ in $\rho^{2}$ and
$z^{2}$. The electric potential described by eq. (10) is characteristic to a
trap that exhibits both axial and radial (with respect to the $xOy$ plan)
symmetry. In case of a Penning trap $A$ is constant in time, but for most RF
(Paul) traps $A\left(t\right)=U_{0}+V_{0}\cos{\Omega t}$, where $U_{0}$ is the
d.c. voltage and $V_{0}$ is the RF trapping voltage. In case of harmonic
potentials $c_{k}=0\ ,\mbox{for}\;k>1$

$g\left(\rho,z\right)=\frac{1}{r_{0}^{2}+2z_{0}^{2}}\left(\rho^{2}-2z^{2}\right)\
,\;c_{2}=-\frac{1}{r_{0}^{2}+2z_{o}^{2}}\ ,$ (12)

where $r_{0}$ and $z_{0}$ represent the semiaxes of the combined quadrupole
trap and $\rho^{2}=x^{2}+y^{2}$. The particular case of an ideal Paul trap is
obtained for zero magnetic field. In case of a RF ion trap with octupole
anharmonicity, the electric potential given by eq. (11) becomes

$g\left(\rho,z\right)=c_{1}H_{2}\left(\rho,z\right)+c_{2}H_{4}\left(\rho,z\right)+\ldots\\\
H_{2}\left(\rho,z\right)=\frac{1}{2}\left(2z^{2}-\rho^{2}\right)\
,\;H_{4}\left(\rho,z\right)=\frac{1}{8}\left(8z^{4}-24z^{2}\rho^{2}+3\rho^{4}\right)\
,$ (13)

where $H_{2}k$ represent the harmonic polynomials [28]. We denote the elastic
constants as

$K_{r}=\frac{M\omega_{c}^{2}}{4}-2Qc_{2}A(t)\ ,\;K_{a}=4Qc_{2}A(t)\ .$ (14)

The quantum Hamilton function that describes an ion (particle) of electric
charge $Q$, mass $M$, and orbital angular momentum $\hbar l$, confined within
an octupole ion trap characterized by a RF anharmonic electric potential, can
be expressed as [32]

$H_{l}=-\frac{\hbar^{2}}{2M}\left(\frac{\partial^{2}}{\partial\rho^{2}}-\frac{l^{2}}{\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}\right)-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
z^{2}}+\frac{K_{r}}{2}\rho^{2}+\frac{K_{a}}{2}z^{2}-\\\
\frac{\omega_{c}}{2}\hbar l+QA\left(t\right)P\left(\rho^{2},z^{2}\right)\ ,$
(15)

where $\omega_{c}=QB_{0}/M$ stands for the cyclotronic frequency, and $B_{0}$
is the axial magnetic field. The anharmonic part can be expressed as

$P\left(\rho^{2},z^{2}\right)=\sum_{k\geq 2}c_{k}H_{2k}\left(\rho,z\right)\;,$
(16)

where $c_{k}$ are constants [32], and $H_{2k}$ are harmonic polynomials
(spherical harmonics functions) of $k$ degree in $\rho^{2}$ and $z^{2}$,
defined as

$H_{2k}\left(\rho,z\right)=\sum_{j=0}^{k}\frac{\left(2k\right)!\rho^{2j}z^{2(k-j)}}{4^{j}\left(2k-2j\right)!\left(j!\right)^{2}}\
.$ (17)

In order to characterize the quantum system, we must build the Hamilton
function. We investigate the Schrödinger equation

$i\hbar\frac{\partial\chi}{\partial t}=H_{l}\chi\ ,$ (18)

where the quantum Hamilton function $H_{l}$ results from eq. (15),
particularized for combined quadrupole or octupole traps. Thus, the quantum
Hamilton function describes an algebraic model when
$P\left(\rho^{2},z^{2}\right)$ is a polynomial function. Such model is linear
for quadrupole traps ($P$ is a linear combination of $\rho^{2}$ and $z^{2}$).
We choose $\mu=-l^{2}$ and we reintroduce the Lie algebra generators for the
radial symplectic group

$K_{0,1r}^{\prime}=\frac{1}{2\hbar\omega_{r}}\left[\pm\frac{\hbar^{2}}{2M}\left(\frac{\partial^{2}}{\partial\rho^{2}}+\frac{\mu}{\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}\right)+\frac{M\omega_{r}^{2}}{2}\rho^{2}\right]\
,$ (19a)
$K_{2r}^{\prime}=\frac{i}{2}\left(\rho\frac{\partial}{\partial\rho}+1\right)\
,$ (19b)

and

$\left[K_{0r}^{\prime},K_{1r}^{\prime}\right]=iK_{2r}^{\prime}\
,\;\left[K_{2r}^{\prime},K_{1r}^{\prime}\right]=iK_{0r}^{\prime}\ ,\\\
\left[K_{2r}^{\prime},K_{0r}^{\prime}\right]=iK_{1r}^{\prime}\ .$ (20)

Using eqs. (19) we infer

$\rho^{2}=\frac{2\hbar}{M\omega_{r}}\left(K_{0r}^{\prime}+K_{1r}^{\prime}\right)\,$
(21a)
$-\frac{\hbar^{2}}{2M}\left(\frac{\partial^{2}}{\partial\rho^{2}}+\frac{\mu}{\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}\right)=\hbar\omega_{r}\left(K_{0r}^{\prime}-K_{1r}^{\prime}\right)\
.$ (21b)

Using eqs. (3), the quantum Hamilton function described by eq. (15) changes
accordingly

$H_{l}=\hbar\omega_{r}\left(K_{0r}^{\prime}-K_{1r}^{\prime}\right)+2\hbar\omega_{a}\left(K_{0a}-K_{1a}\right)+\frac{K_{r}}{2}\frac{2\hbar}{M\omega_{r}}\left(K_{0r}^{\prime}+K_{1r}^{\prime}\right)+\\\
\frac{K_{a}}{2}\frac{2\hbar}{M\omega_{a}}\left(K_{0a}+K_{1a}\right)-\frac{\omega_{c}}{2}\hbar
l+QA\left(t\right)P\left(\rho^{2},z^{2}\right)\ .$ (22)

We consider our trap to exhibit axial symmetry (in the coordinates $\rho$ and
$z$) and rotation symmetry with respect to the $xOy$ plan. From eqs. (20) and
(21) we infer

$\left(K_{0r}^{\prime}+K_{1r}^{\prime}\right)\left(K_{0r}^{\prime}-K_{1r}^{\prime}\right)=C_{2}+K_{2r}^{\prime
2}-iK_{2r}^{\prime}\ ,$ (23)

and the Casimir operator is

$C_{2r}=-\frac{1}{4}\left(\mu+1\right)\Rightarrow k_{r}=\frac{l+1}{2}\
,\;\mbox{with}\;l\geq 0\ .$ (24)

Hence, the study of the Hamilton function for an ion trapped within the
combined (or Paul) trap we have considered is reduced to the study of a linear
Hamilton system for the real symplectic group $Sp\left(2,{\mathbb{R}}\right)$.
Then

$P\left(\rho^{2},z^{2}\right)=DH_{4}\left(\rho,z\right)=D\left(8z^{4}-24z^{2}\rho^{2}+3\rho^{4}\right)$
(25)

where $D$ stands for a coefficient that depends on the trap geometry. We
perform an average on the generators of the symplectic group, and on
$\rho^{4},\rho^{2}z^{2}$ and $z^{4}$. The anharmonic term is

$\Phi_{anh}=A\left(t\right)DH_{4}\left(\rho,z\right)\ ,$ (26)

while in the quasipotential approximation case (ideal Paul trap), the
anharmonic term can be expressed as

$\Phi_{anheff}=C_{4}H_{4}\left(\rho,z\right)+C_{6}H_{6}\left(\rho,z\right)\
,\\\
H_{6}\left(\rho,z\right)=\frac{1}{16}\left(-5\rho^{6}+90\rho^{4}z^{2}-120\rho^{2}z^{4}+16z^{6}\right)$
(27)

where the $C_{4}$ and $C_{6}$ coefficients depend on the specific trap
geometry. Hence, the quantum Hamilton function can be expressed as

$H_{l}=H_{anh}+\left\\{\begin{array}[]{ll}Q\Phi_{anh}&\mbox{Paul and Penning
trap,}\\\ Q\Phi_{anheff}&\,\mbox{ideal Paul trap}\ .\end{array}\right\\}$ (28)

Eq. (28) gives the quantum Hamilton function for the particular traps we have
chosen, including the pseudopotential approximation case that characterizes an
ideal Paul trap.

### Schrödinger equation solutions and energy spectrum

The solutions of the Schrödinger equation for the Hamilton function described
by eq. (22) are [28, 23]

$\Psi_{k_{a}m_{a}k_{r}m_{r}l}=\frac{1}{\sqrt{\rho}}\exp\left[il\left(\theta+\frac{\omega_{c}}{2}t\right)-i\varphi\right]\psi_{k_{a}m_{a}}\left(z_{a}\right)\psi_{k_{r}m_{r}}\left(z_{r}\right)\;,$
(29)

where $m_{a},m_{r}\in{\mathbb{N}}$ are natural numbers,
$z_{a},z_{r}\in{\mathbb{Z}}$ are complex coordinates within the unit disk
$\left(|z_{a}|<1,\;|z_{r}|<1\right)$, and $\varphi\in{\mathbb{R}}$ represents
a real phase such as

$\varphi=\left(k_{a}+m_{a}\right)\varphi_{a}+\left(k_{r},m_{r}\right)\varphi_{r}\;.$
(30)

The $z_{a},z_{r},\varphi_{a}$ and $\varphi_{r}$ variables are solutions of the
differential equations

$i\frac{dz_{j}}{dt}=\alpha_{j}+\frac{\beta_{j}}{2}\left(z_{j}^{2}+1\right)\;,\;\frac{d\varphi_{j}}{dt}=\alpha_{j}+\frac{\beta_{j}}{2}\left(z_{j}+z_{j}^{*}\right)\;,\;j=a,r.$
(31)

The functions $\psi_{k_{a}m_{a}}\left(z_{a}\right)$ and
$\psi_{k_{r}m_{r}}\left(z_{r}\right)$ represent the coherent symplectic
vectors, for symplectic axial and radial dynamic groups. In order to
explicitly define these vectors, we recall the orthonormal system of coherent
states associated to the group SU(1,1) [34]. We use the raising and lowering
operators $K_{\pm j}$ given by eq. (36), and the canonical base [16, 23]

$\phi_{k_{j}m_{j}}=\left[\frac{\Gamma\left(2k_{j}\right)}{m_{j}!\Gamma\left(2k_{j}+m_{j}\right)}\right]^{1/2}\left(K_{+j}\right)^{m_{j}}\phi_{k_{j}0}\;,$
(32)

where the normalized vector satisfies $K_{+j}\phi_{k_{j}0}=k_{j}\phi_{k_{c}0}$
and $K_{-j}\phi_{k_{j}0}=0$. We introduce the unitary evolution operator
defined as

$U\left(z_{j}\right)=\exp\left(z_{j}K_{+j}\right)\exp\left(\lambda_{j}K_{0j}\right)\exp\left(-z_{j}^{*}K_{-j}\right)\;,$
(33)

where $\lambda_{j}=\ln(1-z_{j}z_{j}^{*})$. The coherent symplectic vectors are
defined as

$\psi_{k_{j}m_{j}}\left(z_{j}\right)=U\left(z_{j}\right)\phi_{k_{j}m_{j}}\;.$
(34)

The quasienergy spectrum is thus determined by the
$\Psi_{k_{a}m_{a}k_{r}m_{r}l}$ wavefunctions and the quasienergies [23]

$E_{k_{a}m_{a}k_{r}m_{r}l}=2\hbar\left[\mu_{a}\left(k_{a}+m_{a}\right)+\mu_{r}\left(k_{r}+m_{r}\right)-\frac{\omega_{c}l}{4}\right]\;,$
(35)

where $\mu_{a}$ and $\mu_{r}$ represent the Floquet exponents for the stable
solutions of the Hill equations.

The quantum Hamilton function for a system of $N$ ions confined within a
cylindrically symmetric RF quadrupole ion trap is found in [23]. The center of
mass (CM) Hamilton function is similar to the single particle Hamiltonian
given in [32], which means that the quasienergy and coherent states are almost
the same (the single ion mass is substituted with the product of masses in
case of the multi-ion system). When the interaction potential is invariant to
translations and homogeneous of rank $-2$, soluble models result (e.g., for
Calogero type potentials). Such models enable obtaining explicit bases and
systems of coherent states in order to investigate trapped ion systems [18].

##  3 Quasiclassical dynamics in nonlinear combined and Paul traps

### Classical Hamilton function for an anharmonic octupole and quadrupole
trap. Equations of motion

We recall some definitions with respect to coherent states built over weight
states of discrete positive series of the symplectic group
$Sp\left(2,{\mathbb{R}}\right)$ [35]. The Lie algebra associated to this group
is spanned by three generators, namely $K_{0},K_{1}$, and $K_{2}$, where the
raising and lowering operators $K_{\pm}=K_{1}\pm iK_{2}$, obey the following
commutation relations [28, 17, 19]

$\left[K_{0},K_{\pm}\right]=\pm K_{\pm}\ ,\;\ \left[K_{-},K_{+}\right]=2K_{0}\
,$ (36)

The eigenvalues of the Casimir operator are denoted by $k\left(k-1\right)$,
where $k$ stands for the Bargmann index for unitary irreducible
representations (UIR) [35]. In such case we denote the basis vector as
$\left|m,k\right\rangle;m=0,1,\ldots$, where $K_{0}$ is diagonal. The action
of the operators $K_{+}$ and $K_{-}$ is defined in [19, 35]. We introduce the
following squeezed coherent states for the group
$Sp\left(2,{\mathbb{R}}\right)$:

$\left|z,m,k\right\rangle=U\left(z\right)\left|m,k\right\rangle\ ,$ (37)

where the unitary evolution operator for the quantum Hamiltonian (which
satisfies the Schrödinger equation) is [17]

$U\left(z\right)=\exp\left(zK_{+}\right)\exp\left(\lambda
K_{0}\right)\exp\left(-z^{*}K_{-}\right)\ ,$ (38)

are the operators corresponding to the group representation with
$\left|z\right|<1$ and $\lambda=\ln\left(1-zz^{*}\right)$. For $m=0$, the
geometric construction of Perelomov results [16, 19]. The matrix elements for
the generators $K_{+},K_{-}$, and $K_{0}$ in the squeezed coherent states
$\left|z,k,m\right\rangle$ can be expressed as [19, 32]

$\begin{array}[]{lll}\left\langle
z,k,m\right|K_{+}\left|z,k,m\right\rangle=2\left(k+m\right)z^{*}\left(1-zz^{*}\right)^{-1}&&\\\
\left\langle
z,k,m\right|K_{-}\left|z,k,m\right\rangle=2\left(k+m\right)z\left(1-zz^{*}\right)^{-1}&&\\\
\left\langle
z,k,m\right|K_{0}\left|z,k,m\right\rangle=\left(k+m\right)\left(1+zz^{*}\right)\left(1-zz^{*}\right)^{-1}&&\
.\end{array}$ (39)

We denote $\Omega=K_{0}+K_{1}=K_{0}+\frac{1}{2}\left(K_{+}+K_{-}\right)$. Then

$U^{\dagger}\left(z\right)\Omega
U\left(z\right)=\frac{\left(1+z\right)\left(1+z^{*}\right)}{\left(1-zz^{*}\right)}\left(2K_{0}+e^{-i\varphi}K_{-}+e^{i\varphi}K_{+}\right)\
,$ (40)

where $U^{\dagger}$ represents the Hermitian adjoint of the unitary operator
$U\left(z\right)$, and

$e^{i\varphi}=\frac{1+z}{1+z^{*}}\ ,\;e^{-i\varphi}=\frac{1+z^{*}}{1+z}\
,\;\varphi\in{\mathbb{R}}\ .$ (41)

Moreover

$\left[U^{\dagger}\left(z\right)\Omega
U\left(z\right)\right]^{n}=U^{\dagger}\left(z\right)\Omega^{n}U\left(z\right)=\frac{\left(1+z\right)^{n}\left(1+z^{*}\right)^{n}}{\left(1-zz^{*}\right)^{n}}E^{n}\
,\\\ E=2K_{0}+e^{-i\varphi}K_{-}+e^{i\varphi}K_{+}\ .$ (42)

From eqs. (39 \- 42) we infer

$\left\langle
z,k,m\right|\Omega^{n}\left|z,k,m\right\rangle=\frac{\left(1+z\right)^{n}\left(1+z^{*}\right)^{n}}{\left(1-zz^{*}\right)^{n}}\left\langle
k,m\right|E^{n}\left|k,m\right\rangle\ .$ (43)

Eq. 43 allows the explicit calculus of the energy function for any dynamical
group. By calculus we infer

$E^{2}=4K_{0}^{2}+e^{-2i\varphi}K_{-}^{2}+e^{2i\varphi}K_{+}^{2}+K_{-}K_{+}+K_{+}K_{-}+\\\
2e^{-i\varphi}\left(K_{0}K_{+}+K_{-}K_{0}\right)+2e^{i\varphi}\left(K_{0}K_{+}+K_{+}K_{0}\right)\
,$ (44)
$E^{3}=8K_{0}^{3}+e^{-3i\varphi}K_{-}^{3}+e^{3i\varphi}K_{+}^{3}+2e^{-2i\varphi}\left(K_{-}^{2}K_{0}+K_{0}K_{-}^{2}+K_{-}K_{0}K_{-}\right)\\\
+2e^{2i\varphi}\left(K_{+}^{2}K_{0}+K_{0}K_{+}^{2}+K_{+}K_{0}K_{+}\right)+4e^{-i\varphi}\left(K_{0}K_{-}K_{0}+K_{-}K_{0}^{2}+K_{0}^{2}K_{-}\right)+\\\
4e^{i\varphi}\left(K_{0}K_{+}K_{0}+K_{+}K_{0}^{2}+K_{0}^{2}K_{+}\right)+e^{i\varphi}\left(K_{+}^{2}K_{-}+K_{-}K_{+}^{2}+K_{+}K_{-}K_{+}\right)+\\\
e^{-i\varphi}\left(K_{-}K_{+}K_{-}+K_{+}K_{-}^{2}+K_{-}^{2}K_{+}\right)+\\\
2\left(K_{0}K_{+}K_{-}+K_{+}K_{0}K_{-}+K_{-}K_{+}K_{0}+K_{+}K_{-}K_{0}+K_{0}K_{-}K_{+}+K_{-}K_{0}K_{+}\right)\
.$ (45)

We weight on the vacuum state vector $\left|0\right\rangle$ and obtain

$K_{-}\left|0\right\rangle=0\ ,\left\langle 0|0\right\rangle=0\
,K_{0}\left|0\right\rangle=k\left|0\right\rangle\ ,\left\langle
0\right|K_{+}=0\ .$ (46)

Eq. (36) leads to

$K_{-}K_{+}=K_{+}K_{-}+2K_{0}\ ,K_{-}K_{0}=K_{0}K_{-}+K_{-}\ ,\\\
K_{0}K_{+}=K_{+}K_{0}+K_{+}\ .$ (47)

Then, the Casimir operator can be expressed as

$C_{2}=K_{0}^{2}-\frac{1}{4}\left(K_{+}+K_{-}\right)^{2}-\frac{1}{4i^{2}}\left(K_{+}-K_{-}\right)^{2}=k\left(k-1\right)\
,$

and

$K_{+}K_{-}+K_{-}K_{+}=2K_{0}^{2}-2k\left(k-1\right)\ .$ (48)

Using eqs. (46) we infer

$\left\langle E^{2}\right\rangle=\left\langle
0|E^{2}|0\right\rangle=4k^{2}+2k\ ,$ (49a) $\left\langle
E^{3}\right\rangle=\left\langle 0|E^{3}|0\right\rangle=8k^{3}+12k^{2}+4k\ .$
(49b)

If instead of the vacuum state vector $\left|0\right\rangle$ we weight on a
vector of an orthonormal system of state vectors $\left|k,m\right\rangle$,
namely $\left\langle k,m|k,m^{\prime}\right\rangle=\delta_{mm^{\prime}}$, we
obtain

$\left\langle k,m|E^{2}|k,m\right\rangle=2k\left(2k+1\right)+12km+6m^{2}\ ,$
(50a) $\left\langle
k,m|E^{3}|k,m\right\rangle=4k\left(k+1\right)\left(2k+1\right)+4mk\left(5+12k\right)+4m^{2}\left(15k+1\right)+20m^{3}\
,$ (50b)

where $K_{0}\left|k,m\right\rangle=\left(k+m\right)\left|k,m\right\rangle$
[19] and we make use of eqs. (46) and (49). Using eq. (47), we are able to
evaluate the final term from eq. (45) as

$F=2K_{0}\left(K_{+}K_{-}+K_{-}K_{+}\right)+2\left(K_{+}K_{-}+K_{-}K_{+}\right)K_{0}+\\\
2K_{0}\left(K_{+}K_{-}+K_{-}K_{+}\right)+2\left(K_{+}K_{-}+K_{-}K_{+}\right)\
.$ (51)

We denote

$\xi=\frac{\left(1+z\right)\left(1+z^{*}\right)}{1-zz^{*}}=S_{n}\
,\;n=1,2,\ldots\ ,$ (52)

where

$S_{1}=2\xi\left(k+m\right)=S_{1}\left(z,k,m\right)$ (53a)
$S_{2}=\xi^{2}\left[2k\left(2k+1\right)+12km+6m^{2}\right]$ (53b)
$S_{3}=\xi^{3}\left[4k\left(k+1\right)\left(2k+1\right)+4mk\left(5+12k\right)+4m^{2}\left(15k+1\right)+20m^{3}\right]\
.$ (53c)

By taking into account eq. (39) [28, 19] and

$\left\langle z,k,m|A|z,k,m\right\rangle=\left\langle
k,m|U^{-1}\left(z\right)AU\left(z\right)|k,m\right\rangle\ ,$ (54)

we infer

$\left\langle
z,k,m|K_{0}-K_{1}|z,k,m\right\rangle=\left(k+m\right)\frac{\left(1-z\right)\left(1-z^{*}\right)}{1-zz^{*}}\
,$ (55)

where $z$ stands for the squeezed state parametre, $k$ and $m$ represent pure
harmonic oscillator states, while $K_{0}-K_{1}$ is the kinetic energy. The
classical energy (Hamilton) function for an anharmonic octupole ion trap
described by eq. (22), can be expressed as

${\mathcal{H}}_{cl}=\hbar\omega_{r}\left(k_{r}+m_{r}\right)\frac{\left(1-z_{r}\right)\left(1-z^{*}_{r}\right)}{1-z_{r}z^{*}_{r}}+2\hbar\omega_{a}\left(k_{a}+m_{a}\right)\frac{\left(1-z_{a}\right)\left(1-z^{*}_{a}\right)}{1-z_{a}z^{*}_{a}}+\\\
\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)\frac{\left(1-z_{r}\right)\left(1-z^{*}_{r}\right)}{1-z_{r}z^{*}_{r}}+\frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)\frac{\left(1-z_{a}\right)\left(1-z^{*}_{a}\right)}{1-z_{a}z^{*}_{a}}+H_{anh}-\frac{\omega_{c}}{2}\hbar
l\ ,$ (56)

where

$H_{anh}=QA\left(t\right)\langle
P\left(\rho^{2},z^{2}\right)\rangle=QA\left(t\right)\sum_{k\geq 1}c_{k}\langle
H_{2k}\rangle.$ (57)

We have denoted by $\langle X\rangle$ the expectation value of the $X$
operator in the coherent state $\Phi_{m_{a}m_{r}}\left(z_{a},z_{r}\right)$. In
particular we obtain

$H_{anh}=\left\\{\begin{array}[]{ll}QA\left(t\right)D\left\langle
H_{4}\right\rangle\ ,\mbox{combined (Paul and Penning) trap}\ ,&\\\
QC_{4}\left\langle H_{4}\right\rangle+QC_{6}\left\langle
H_{6}\right\rangle,\mbox{ideal Paul trap}&\ .\end{array}\right.$ (58)

The second term in eq. (58) characterizes the pseudopotential aproximation
case (the ideal radiofrequency or Paul trap) [28]. Then

$\left\langle H_{4}\right\rangle=8S_{2a}-24S_{1r}S_{1a}+3S_{2r}\ ,$ (59)
$\left\langle
H_{6}\right\rangle=16S_{3a}-120S_{2a}S_{1r}+90S_{1a}S_{2r}-5S_{3r}\ ,$ (60)
$H_{6}=16z^{6}-120z^{4}\rho^{2}+90z^{2}\rho^{4}-5\rho^{6}\ .$ (61)

We denote

$S_{jr}=S_{j}\left(z_{r},k_{r},m_{r}\right)\ \
,\;\;S_{ja}=S_{j}\left(z_{a},k_{a},m_{a}\right)\ \ ,\;\;1\leq j\leq 3\ ,$

and then introduce the Husimi $Q$ representation

$S_{j}\left(z,k,m\right)=\xi^{j}Q_{j}\left(k,m\right)\ .$ (62)

Moreover

$Q_{1}\left(k,m\right)=2\left(k+m\right)$ (63a)
$Q_{2}\left(k,m\right)=2k\left(2k+1\right)+12km+6m^{2}$ (63b)
$Q_{3}\left(k,m\right)=4k\left(k+1\right)\left(2k+1\right)+4mk\left(5+12k\right)+4m^{2}\left(15k+1\right)+20m^{3}\
.$ (63c)

The parametre values are known

$k_{a}=\frac{1}{4},\frac{3}{4}\
;\;\;m_{a},m_{r}=0,1,\ldots\;;\;\;k_{r}=\frac{l+1}{2}\ ,$

where $l$ stands for the orbital quantum number for the trapped particle
(ion), and the orbital angular momentum is $L^{2}=l\left(l+1\right)I$. By
minimizing ${\mathcal{H}}$ after $z$ we infer the approximations for the
quantum energies. Using cylindrical coordinates $\left(z=\rho e^{i\theta}\
,\;z^{*}=\rho e^{-i\theta}\right)$, we introduce the complex variables

$\xi_{n}=\frac{\left(1+z\right)\left(1+z^{*}\right)}{1-zz^{*}}\
,\;\eta_{n}=\frac{\left(1-z\right)\left(1-z^{*}\right)}{1-zz^{*}}\ ,\;n=a,r\
.$ (64)

The classical Hamilton function for an anharmonic octupole trap can be
expressed as

${\mathcal{H}}_{cl}=\hbar\omega_{r}\left(k_{r}+m_{r}\right)\eta_{r}+2\hbar\omega_{a}\left(k_{a}+m_{a}\right)\eta_{a}+\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)\xi_{r}+\\\ \frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)\xi_{a}-\frac{\omega_{c}}{2}\hbar
l\\\
+\left\\{\begin{array}[]{ll}QA\left(t\right)D\left[8S_{2a}-24S_{1r}S_{1a}+3S_{2r}\right]\longmapsto\mbox{combined
trap,}&\\\ QC_{4}\left(8S_{2a}-24S_{1r}S_{1a}+3S_{2r}\right)+\\\
QC_{6}\left(16S_{3a}-120S_{2a}S_{1r}+90S_{1a}S_{2r}-5S_{3r}\right)\longmapsto\mbox{pseudopot.
approx.}\ ,\end{array}\right.$ (65)

where the first expression refers to the case of an octupole combined Paul and
Penning trap, while the second case describes the pseudopotential
approximation in case of an ideal Paul trap. The energy function associated to
the quantum Hamilton function $H_{l}$ is a classical type Hamiltonian
$H_{cl}$, whose values are exactly the expectation values of $H_{l}$ on the
coherent states $\psi_{k_{a}0}(z_{a})$ and $\psi_{k_{r}0}(z_{r})$. Therefore

$\frac{\partial}{\partial z}=\frac{\partial\xi}{\partial
z}\frac{\partial}{\partial\xi}+\frac{\partial\eta}{\partial
z}\frac{\partial}{\partial\eta}\;\ ,\;\;\frac{\partial}{\partial
z^{*}}=\frac{\partial\xi}{\partial
z^{*}}\frac{\partial}{\partial\xi}+\frac{\partial\eta}{\partial
z^{*}}\frac{\partial}{\partial\eta}\ .$ (66)

Quantum dynamics for ions confined in quadrupole and octupole combined and
ideal Paul traps can be characterized by solutions of the Schrödinger
equation, where the $H=H_{r}$ Hamilton function corresponds to radial motion,
while $H=H_{a}$ describes the axial motion. We apply the TDVP and obtain the
classical Hamilton function for the particular traps we considered

${\mathcal{H}_{cl}}=\left\langle z,k,m|H_{l}|z,k,m\right\rangle\ ,$

which determines an equation of motion in the classical Lobacevski phase space
$\left|z\right|<1$ :

$\dot{z}=\left\\{z,{\mathcal{H}_{cl}}\right\\}\ ,$ (67)

where $\left\\{,\right\\}$ stands for the generalized Poisson bracket. We
finally obtain

$\left\\{z,{\mathcal{H}_{cl}}\right\\}=\frac{\left(1-zz^{*}\right)^{2}}{2i\left(k+m\right)}\frac{\partial{\mathcal{H}_{cl}}}{\partial
z^{*}}\ ,$ (68)

which represents the equation of motion for an ion confined within an octupole
combined (or Paul) trap. Thus, we show that the coherent state approach leads
to quantum solutions that are explicitly constructed as functions of the
classical trajectories on the phase space.

### Hamilton function for an anharmonic combined trap. Equations of motion

According to eq. (65) the classical Hamilton function can be expressed as

${\mathcal{H}}_{cl\
anh}=\hbar\omega_{r}\left(k_{r}+m_{r}\right)\eta_{r}+2\hbar\omega_{a}\left(k_{a}+m_{a}\right)\eta_{a}+\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)\xi_{r}+\\\ \frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)\xi_{a}-\frac{\omega_{c}}{2}\hbar
l+QA\left(t\right)D\left[8S_{2a}-24S_{1r}S_{1a}+3S_{2r}\right]\ .$ (69)

The expressions for $S_{1a},S_{2a},S_{1r}$ and $S_{2r}$ result from eqs.
(53a). We turn to the expression of the Hamilton function described by eq.
(69). The points of minimum of the Hamilton function are characterized by the
following equations

$\frac{\partial{\mathcal{H}}_{lanh}}{\partial\xi_{a}}=0\Rightarrow\frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)+\\\
QA\left(t\right)D\left[16\xi_{a}\cdot
Q_{2}\left(k_{a},m_{a}\right)-96\xi_{r}\left(k_{a}+m_{a}\right)\left(k_{r}+m_{r}\right)\right]=0$
(70a)
$\frac{\partial{\mathcal{H}}_{lanh}}{\partial\xi_{r}}=0\Rightarrow\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)+\\\
QA\left(t\right)D\left[-96\xi_{a}\left(k_{a}+m_{a}\right)\left(k_{r}+m_{r}\right)+6\xi_{r}Q_{2}\left(k_{r},m_{r}\right)\right]=0$
(70b) $\frac{\partial{\mathcal{H}}_{lanh}}{\partial\eta_{a}}=0\Rightarrow
2\hbar\omega_{a}\left(k_{a}+m_{a}\right)=0$ (70c)
$\frac{\partial{\mathcal{H}}_{lanh}}{\partial\eta_{r}}=0\Rightarrow
2\hbar\omega_{r}\left(k_{r}+m_{r}\right)=0\ .$ (70d)

Using eq. (64), we infer

$1-zz^{*}=2\frac{1+zz^{*}}{\xi+\eta}\
,\;z+z^{*}=\frac{1}{2}\left(1-zz^{*}\right)\left(\xi-\eta\right)\ .$ (71)

According to eq. (66) and eqs. (64)

$\frac{\partial{\mathcal{H}}}{\partial
z^{*}}=\frac{\partial{\mathcal{H}}}{\partial\xi}\left(\frac{1+z}{1-zz^{*}}\right)^{2}-\frac{\partial{\mathcal{H}}}{\partial\eta}\left(\frac{z-1}{1-zz^{*}}\right)^{2}$
(72)

From eqs. (64) and (71), after some calculus, we obtain

$\dot{\xi}=\left(\frac{1+z^{*}}{1-zz^{*}}\right)^{2}\dot{z}+\left(\frac{1+z}{1-zz^{*}}\right)^{2}z^{*}\
.$ (73a)
$\dot{\eta}=-\left(\frac{1-z^{*}}{1-zz^{*}}\right)^{2}\dot{z}-\left(\frac{1-z}{1-zz^{*}}\right)^{2}z^{*}\
.$ (73b)

From eq. (67), we can express $z^{*}$ as:

$z^{*}=-\frac{\left(1-zz^{*}\right)^{2}}{2i\left(k+m\right)}\frac{\partial{\mathcal{H}}}{\partial
z}\ ,$ (74)

while eqs. (66) yields

$\frac{\partial{\mathcal{H}}}{\partial
z}=\frac{\partial{\mathcal{H}}}{\partial\xi}\left(\frac{1+z^{*}}{1-zz^{*}}\right)^{2}-\frac{\partial{\mathcal{H}}}{\partial\eta}\left(\frac{z^{*}-1}{1-zz^{*}}\right)^{2}\
.$ (75)

If we replace eqs. (67) and (74) in eq. (73) and do the math, we infer

$2i\left(k+m\right)\dot{\xi}=4\frac{z-z^{*}}{1-zz^{*}}\frac{\partial{\mathcal{H}}}{\partial\eta}\
.$ (76)

If the same mathematical operations are applied to eq. (73b), the result is

$2i\left(k+m\right)\dot{\eta}=-4\frac{z-z^{*}}{1-zz^{*}}\frac{\partial{\mathcal{H}}}{\partial\xi}\
.$ (77)

We return to eq. (64) and then write eq. (71) as

$\frac{2}{1-zz^{*}}=\frac{\xi+\eta}{2}+1\ .$ (78)

Then

$z+z^{*}=2\frac{\xi-\eta}{\xi+\eta+2}\
,\;zz^{*}=\frac{\xi+\eta-2}{\xi+\eta+2}\ .$ (79)

Using eqs. (79) we infer

$\left(z-z^{*}\right)^{2}=16\frac{1-\xi\eta}{\left(\xi+\eta+2\right)^{2}}\ .$
(80)

Hence

$\frac{z-z^{*}}{1-zz^{*}}=i\varepsilon\sqrt{\left|1-\xi\eta\right|}\;\;,\;\varepsilon=\pm
1\ .$ (81)

Considering eq. (81), the equations of motion (76) and (77) can be expressed
as

$2\left(k+m\right)\dot{\eta}=-4\varepsilon\sqrt{\left|1-\xi\eta\right|}\frac{\partial{\mathcal{H}}}{\partial\xi}$
(82a)
$2\left(k+m\right)\dot{\xi}=4\varepsilon\sqrt{\left|1-\xi\eta\right|}\frac{\partial{\mathcal{H}}}{\partial\eta}\
.$ (82b)

Further on, we denote

$i\sigma=\frac{z-z^{*}}{1-zz^{*}}\Rightarrow-\sigma^{2}=1-\xi\eta\;,$ (83)

while the time derivative is

$2\sigma\dot{\sigma}=\dot{\xi}\eta+\xi\dot{\eta}\ .$ (84)

Then eqs. (82a) and (82b) can be expressed as

$\left(k+m\right)\dot{\eta}=-2\sigma\frac{\partial{\mathcal{H}}}{\partial\xi}\
,\;\left(k+m\right)\dot{\xi}=2\sigma\frac{\partial{\mathcal{H}}}{\partial\eta}\;.$
(85)

By comparing eqs. (85) and (84), we infer the following equations of motion

$\left(k+m\right)\dot{\sigma}=\eta\frac{\partial{\mathcal{H}}}{\partial\eta}-\xi\frac{\partial{\mathcal{H}}}{\partial\xi}\
.$ (86)

### Quasipotential approximation case - Ideal RF (Paul) trap case

The expression of the classical Hamilton function for an octupole trap, in
case of the pseudopotential approximation is (see eq. (65))

${\mathcal{H}}_{cl\
anheff}=\hbar\omega_{r}\left(k_{r}+m_{r}\right)\eta_{r}+2\hbar\omega_{a}\left(k_{a}+m_{a}\right)\eta_{a}+\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)\xi_{r}+\\\ \frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)\xi_{a}-\frac{\omega_{c}}{2}\hbar
l+QC_{4}<H_{4}>+QC_{6}<H_{6}>$ (87)

where $H_{4}$ and $H_{6}$ are characterized by eqs. (59 \- 60). Hence

$S_{3a}=4k_{a}\left(k_{a}+1\right)\left(2k_{a}+1\right)+4k_{a}m_{a}\left(5+12k_{a}\right)+4m_{a}^{2}\left(15k_{a}+1\right)+20m_{a}^{3}\
,$ (88)
$S_{3r}=4k_{r}\left(k_{r}+1\right)\left(2k_{r}+1\right)+4k_{r}m_{r}\left(5+12k_{r}\right)+4m_{r}^{2}\left(15k_{a}+1\right)+20m_{r}^{3}\
.$ (89)

The points of minimum that characterize the equilibrium configurations for
trapped ions are given by the equations

$\frac{2\hbar
K_{a}}{M\omega_{a}}\left(k_{a}+m_{a}\right)+QC_{4}\left[16\xi_{a}Q_{2}\left(k_{a},m_{a}\right)-96\xi_{r}\left(k_{a}+m_{a}\right)\left(k_{r}+m_{r}\right)\right]+\\\
QC_{6}[48\xi_{a}^{2}Q_{3}\left(k_{a},m_{a}\right)-240\xi_{a}\xi_{r}Q_{2}\left(k_{a},m_{a}\right)Q_{1}\left(k_{r},m_{r}\right)+90\xi_{r}Q_{1}\left(k_{a},m_{a}\right)Q_{2}\left(k_{r},m_{r}\right)]=0$
(90) $\frac{2\hbar
K_{r}}{M\omega_{r}}\left(k_{r}+m_{r}\right)+QC_{4}\left[-96\xi_{a}\left(k_{a}+m_{a}\right)\left(k_{r}+m_{r}\right)+6\xi_{r}Q_{2}\left(k_{r},m_{r}\right)\right]+\\\
QC_{6}[-120\xi_{a}^{2}Q_{2}\left(k_{a},m_{a}\right)Q_{1}\left(k_{r},m_{r}\right)+180\,\xi_{a}\xi_{r}Q_{1}\left(k_{a},m_{a}\right)Q_{2}\left(k_{r},m_{r}\right)-15\xi_{r}^{2}Q_{3}\left(k_{r},m_{r}\right)]=0\
.$ (91)

In the pseudopotential approximation case, the points of minimum of the
dequantified Hamilton function define the equilibrium configurations for
trapped ions (Coulomb crystals), of interest for implementing and scaling
quantum logic (see [36]). The equilibrium configurations are found by
performing numerical simulations, which is verified for example in case of a
2D octupole Paul trap, where axially symmetric structures of Coulomb crystals
and phase transitions have been reported [37, 38]. The advantage of a linear
RF octupole trap lies in a large region within the trap centre where the
electric potential is flat, which makes it very suited for cold molecule
spectroscopy or frequency metrology [39]. Moreover, experiments performed with
higher order multipole 2D RF ion traps show the occurrence of multiple regions
of dynamical stability, some of them located near the trap electrodes where
the electrical potential is steep. Numerical simulations performed validate
experimental data and observations [40]. Recent investigations are focused on
producing planar Coulomb crystals using specially designed linear multipole
traps [41, 37, 42]. Calculations performed for time-dependent dynamics in RF
ion traps were also generalized to multispecies ion crystals in multipole
traps [43].

The possibility to perform precision numerical simulations of the associated
(semi)classical or quantum dynamics for any kind of engineered quantum system
is an issue of large importance, and the methods presented in [44] can be
directly adapted to any given Hamiltonian, including the model we propose.
Numerical simulations can be used to validate the method suggested in the
paper and find the equilibrium configurations for trapped ions or the
eigenvalues of the normal modes of oscillation [39]. The results in the paper
are also valid for different quantum simulator architectures [45].

##  4 Conclusions

The paper characterizes the evolution of squeezed coherent states for
quadrupole and octupole RF and combined traps, an issue not studied in the
literature. We propose a dequantization algorithm by which the energy function
associated to the trapped ion system is a classical Hamilton function, whose
values represent the expectation values of the quantum Hamiltonian on squeezed
coherent states. Such formalism can be applied to Hamiltonians that are
nonlinear in the infinitesimal generators of a dynamical symmetry group, such
as the case of quadrupole and octupole ion traps with a RF anharmonic electric
potential, which we demonstrate in the paper. An algorithm is proposed, that
extends the results to any electrical multipole. We build the quantum Hamilton
function associated to an ion confined within a nonlinear combined (or RF) ion
trap. Moreover, we also build the classical Hamilton functions particularized
for such traps. The classical (Hamilton) equations of motion in complex
coordinates are obtained for an ion confined in an octupole combined trap. The
results obtained are also valid in case of the center of mass (CM) motion, for
a system made of identical ions.

To summarize, the coherent state formalism has been developed, particularized
and applied in case of quadrupole and octupole RF and combined traps. We show
that the dynamic group associated to an ion confined within such particular
traps is the product of the symplectic axial group, a symplectic radial group
and the group of rotations around the symmetry axis. The study of the Hamilton
function for an ion trapped in an anharmonic combined (or RF) trap we have
considered is reduced to the study of the linear Hamilton system for the real
symplectic group $Sp\left(2,{\mathbb{R}}\right)$. In agreement with the TDVP
applied on coherent states, if a harmonic analytical electric potential is
added, the classical equations of motion turn into coupled nonlinear Hill
equations. We also show that the coherent state approach leads to quantum
solutions that are explicitly constructed as functions of the classical
trajectories on the Lobacevski phase space.

Dequantization by means of coherent states enables achieving exact connections
between quantum and classical dynamical systems, that exhibit (axial)
symmetry. For an ion confined within a quadrupole ion trap (QIT), the
requantization operation corresponds to a geometrical quanzation. The
method(s) and formalism suggested are developed and implemented for the
particular traps we considered. The points of minimum of the dequantified
Hamilton function can be identified through standard numerical programming.
The algorithm we suggest, combined with the latest experimental techniques,
can be applied to investigate entanglement, quantum logic, and quantum
simulations in a wide range of physical systems.

##  5 acknowledgements

The author would like to acknowledge support provided by the Ministery of
Education, Research and Inovation from Romania (ANCS-National Agency for
Scientific Research), contracts PN 16470103, UEFISCDI 90/06.10.2011, and ROSA
Contract Nr. 53 / 20.11.2013 - LEOPARD.

The author (B. M.) is grateful to Prof. Marian Apostol (IFIN-HH) and Prof.
Günther Werth for very useful comments and suggestions on the paper.

## References

## References

  * [1] M. Orszag, Quantum Optics: Including Noise Reduction, Trapped Ions, Quantum Trajectories, and Decoherence, Springer Intl. Publishing, 2016.
  * [2] K. Blaum, F. Herfurth, Trapped Charged Particles and Fundamental Interactions, Vol. 749 of Lecture Notes in Physics, Springer, Berlin Heidelberg, 2008.  doi:10.1007/978-3-540-77817-2.
  * [3] S. Haroche, J.-M. Raimond, Exploring the Quantum: Atoms, Cavities and Photons, Oxford University Press, 2013.
  * [4] D. Leibfried, R. Blatt, C. Monroe, D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75 (2003) 281–324.  doi:https://doi.org/10.1103/RevModPhys.75.281.
  * [5] D. J. Wineland, Nobel Lecture: Superposition, entanglement, and raising Schrödinger’s cat, Rev. Mod. Phys. 85 (3) (2013) 1103–1114.  doi:https://doi.org/10.1103/RevModPhys.85.1103.
  * [6] M. Knoop, N. Madsen, R. C. Thompson (Eds.), Physics with Trapped Charged Particles: Lectures from the Les Houches Winter School, Imperial College Press, 2014.
  * [7] A. Sinclair, An Introduction to Trapped Ions, Scalability and Quantum Metrology, in: E. Andersson, P. Öhberg (Eds.), Quantum Information and Coherence, Vol. 67 of Scottish Graduate Series, Springer, 2011, pp. 211–246.  doi:10.1007/978-3-319-04063-9.
  * [8] R. Blatt, C. F. Roos, Quantum simulations with trapped ions, Nature Phys. 8 (2012) 277–284.  doi:http://dx.doi.org/10.1038/nphys2252.
  * [9] M. Knoop, I. Marzoli, G. Morigi (Eds.), Ion Traps for Tomorrow Applications, Vol. Course 189 of Proc. Int. School of Physics Enrico Fermi, IOS Press and Società Italiana di Fisica, Bologna, Italy, 2015.
  * [10] B. C. Sawyer, J. W. Britton, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, J. J. Bollinger, Spectroscopy and Thermometry of Drumhead Modes in a Mesoscopic Trapped-Ion Crystal Using Entanglement, Phys. Rev. Lett. 108 (21) (2012) 213003.  doi:https://doi.org/10.1103/PhysRevLett.108.213003.
  * [11] H.-Y. Lo, D. Kienzler, L. de Clercq, M. Marinelli, V. Negnevitsky, B. C. Keitch, J. P. Home, Spin-motion entanglement and state diagnosis with squeezed oscillator wavepackets, Nature 521 (7552) (2015) 336–339.  doi:https://doi.org/10.1038/nature14458.
  * [12] M. Combescure, A quantum particle in a quadrupole radio-frequency trap, Ann. Inst. H. Poincaré A: Physique théorique 44 (3) (1986) 293–314.
  * [13] M. Combescure, The Quantum Stability Problem for Time-Periodic Perturbations of the Harmonic Oscillator, Ann. Inst. H. Poincaré A: Physique théorique 47 (1) (1987) 63–83.
  * [14] M. Combescure, The quantum stability problem for some class of time-dependent hamiltonians, Ann. Phys. 185 (1) (1988) 86–110.  doi:https://doi.org/10.1016/0003-4916(88)90259-X.
  * [15] A. Davydov, On Close Relationship Between Classical Time-Dependent Harmonic Oscillator and Non-relativistic Quantum Mechanics in One Dimension, Int. J. Theor. Phys. 50 (5) (2011) 1451–1467.  doi:http://dx.doi.org/10.1007/s10773-010-0654-1.
  * [16] A. Perelomov, in: W. Beiglböck, J. L. Birman, T. Regge, E. H. Lieb, W. Thirring (Eds.), Generalized Coherent States and Their Applications, Vol. Texts and Monographs in Physics, Springer-Verlag, Berlin, 1986.
  * [17] M. Combescure, A generalized coherent state approach of the quantum dynamics for suitable time-dependent hamiltonians, Journées équations aux dérivées partielles (1992) 1–13.
  * [18] M. Combescure, D. Robert, Coherent States and Applications in Mathematical Physics, Springer, 2012.
  * [19] V. N. Gheorghe, F. Vedel, Quantum dynamics of trapped ions, Phys. Rev. A 45 (1992) 4828–4831.  doi:10.1103/PhysRevA.45.4828.
  * [20] G.-Z. Li, G. Werth, The combined trap and some possible applications, Phys. Scr. 46 (6) (1992) 587.  doi:https://doi.org/10.1088/0031-8949/46/6/017.
  * [21] J. Walz, I. Siemers, M. Schubert, W. Neuhauser, R. Blatt, Motional Stability of a Nonlinear Parametric Oscillator: Ion Storage in the RF Octupole Trap, Europhys. Lett. 21 (2) (1993) 183.  doi:https://doi.org/10.1209/0295-5075/21/2/011.
  * [22] B. M. Mihalcea, C. M. Niculae, V. N. Gheorghe, On the Multipolar Electromagnetic Traps, Rom. Journ. Phys. 44 (5–6) (1999) 543–550.
  * [23] V. N. Gheorghe, G. Werth, Quasienergy states of trapped ions, Eur. Phys. J. D 10 (2000) 197–203.  doi:10.1007/s100530050541.
  * [24] T. Nakamura, S. Ohtani, M. Wada, K. Okada, I. Katayama, H. A. Schuessler, Ion dynamics and oscillation frequencies in a linear combined trap, J. Appl. Phys. 89 (5) (2001) 2922–2931.  doi:https://doi.org/10.1063/1.1345514.
  * [25] A. Contreras-Astorga, D. J. F. C., M. Velázquez, Coherent states for quadratic Hamiltonians, J. Phys. A: Math. Theor. 44 (3) (2010) 035304.  doi:https://doi.org/10.1088/1751-8113/44/3/035304.
  * [26] V. G. Ibarra-Sierra, J. C. Sandoval-Santana, J. L. Cardoso, A. Kunold, Lie algebraic approach to the time-dependent quantum general harmonic oscillator and the bi-dimensional charged particle in time-dependent electromagnetic fields, Ann. Phys. 362 (2015) 83–117.  doi:http://dx.doi.org/10.1016/j.aop.2015.07.022.
  * [27] J. R. Klauder, B. S. Skagerstam, Coherent States: Applications in Physics and Mathematical Physics, World Scientific, Singapore, 1985.
  * [28] F. G. Major, V. N. Gheorghe, G. Werth, Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement, 1st Edition, Vol. 37 of Springer Series on Atomic, Optical and Plasma Physics, Springer-Verlag, Berlin Heidelberg, 2005.
  * [29] P. Kramer, M. Saraceno (Eds.), Geometry of the Time-Dependent Variational Principle in Quantum Mechanics, Vol. 140 of Lecture Notes in Physics, Springer, Berlin, 1981.
  * [30] M. M. Nieto, D. R. Truax, Coherent states sometimes look like squeezed states and vice versa: the Paul trap, New J. Phys. 2 (2000) 18.1–18.9.  doi:https://doi.org/10.1088/1367-2630/2/1/318.
  * [31] A. A. Abrikosov, Jr., E. Gozzi, D. Mauro, Geometric dequantization, Ann. Phys. 317 (1) (2005) 24–71.  doi:https://doi.org/10.1016/j.aop.2004.12.001.
  * [32] B. M. Mihalcea, Semiclassical dynamics for an ion confined within a nonlinear electromagnetic trap, Phys. Scr. T143 (2011) 014018.  doi:http://dx.doi.org/10.1088/0031-8949/2011/T143/014018.
  * [33] P. Busch, P. Lahti, J.-P. Pellonpää, K. Ylinen, Quantum Measurement, Theoretical and Mathematical Physics, Springer International Publishing, Switzerland, 2016.  doi:10.1007/978-3-319-43389-9.
  * [34] K. Nemoto, Generalized coherent states for SU(n) systems, J. Phys. A: Math. Gen. 33 (17) (2000) 3493–3506.  doi:http://dx.doi.org/10.1088/0305-4470/33/17/307.
  * [35] A. J. Dragt, E. Forest, K. B. Wolf, Foundations of a Lie algebraic theory of geometrical optics, in: J. S. Mondragón, K. B. Wolf (Eds.), Lie Methods in Optics, Vol. 250 of Lecture Notes in Physics, Springer, Berlin Heidelberg, 1986, pp. 105–157.  doi:10.1007/3-540-16471-5_4.
  * [36] M. Marciante, C. Champenois, A. Calisti, J. Pedregosa-Gutierrez, M. Knoop, Ion dynamics in a linear radio-frequency trap with a single cooling laser, Phys. Rev. A 82 (3) (2010) 033406.  doi:https://doi.org/10.1103/PhysRevA.82.033406.
  * [37] K. Okada, T. Takayanagi, M. Wada, S. Ohtani, H. A. Schuessler, Observation of ion Coulomb crystals in a cryogenic linear octupole rf ion trap, Phys. Rev. A 80 (4) (2009) 043405.  doi:https://doi.org/10.1103/PhysRevA.80.043405.
  * [38] M. Marciante, C. Champenois, A. Calisti, M. Knoop, Structural phase transitions in multipole traps, Appl. Phys. B 107 (4) (2012) 1117–1123.  doi:https://doi.org/10.1007/s00340-012-4883-2.
  * [39] C. Champenois, M. Marciante, J. Pedregosa-Gutierrez, M. Houssin, M. Knoop, M. Kajita, Ion ring in a linear multipole trap for optical frequency metrology, Phys. Rev. A 81 (4) (2010) 043410.  doi:https://doi.org/10.1103/PhysRevA.81.043410.
  * [40] B. M. Mihalcea, L. C. Giurgiu, C. Stan, G. T. Visan, M. Ganciu, V. Filinov, D. Lapitsky, L. Deputatova, R. Syrovatka, Multipole electrodynamic ion trap geometries for microparticle confinement under standard ambient temperature and pressure conditions, J. Appl. Phys. 119 (11) (2016) 114303.  doi:https://doi.org/10.1063/1.4943933.
  * [41] D. J. Wineland, D. Leibfried, Quantum information processing and metrology with trapped ions, Laser Phys. Lett. 8 (3) (2011) 175–188.  doi:https://doi.org/10.1002/lapl.201010125.
  * [42] C. Champenois, J. Pedregosa-Gutierrez, M. Marciante, D. Guyomarc’h, M. Houssin, M. Knoop, A Double Ion Trap for Large Coulomb Crystals, Vol. 1521 of AIP Conf. Proc., AIP, 2013, pp. 210–219.  doi:https://doi.org/10.1063/1.4796077.
  * [43] H. Landa, M. Drewsen, B. Reznik, A. Retzker, Modes of oscillation in radiofrequency Paul traps, New J. Phys. 14 (9) (2012) 093023.  doi:https://doi.org/10.1088/1367-2630/14/9/093023.
  * [44] K. Singer, U. Poschinger, M. Murphy, P. Ivanov, F. Ziesel, T. Calarco, F. Schmidt-Kaler, Colloquium: Trapped ions as quantum bits: Essential numerical tools, Rev. Mod. Phys. 82 (3) (2010) 2609–2632.  doi:https://doi.org/10.1103/RevModPhys.82.2609.
  * [45] B. Yoshimura, M. Stork, D. Dadic, W. C. Campbell, J. K. Freericks, Creation of two-dimensional Coulomb crystals of ions in oblate Paul traps for quantum simulations, EPJ Quantum Technology 2 (1) (2015) 2.  doi:https://doi.org/10.1140/epjqt14.


# Laser dynamics in nonlinear transparent media with electron plasma
generation: effects of electron-hole radiative recombinations

P. Kameni Nteutse1, Alain M. Dikandé2 and S. Zekeng1 1 Laboratory of
Mechanics, Materials and Structures, Department of Physics, Faculty of
Science, University of Yaounde I P.O. Box 812 Yaounde, Cameroon.  2Laboratory
of Research on Advanced Materials and Nonlinear Science (LaRAMaNS), Department
of Physics, Faculty of Science, University of Buea P.O. Box 63 Buea, Cameroon.
<EMAIL_ADDRESS>
###### Abstract

The performance of optical devices manufactured via laser micromachining on
nonlinear transparent materials, relies usually on three main factors which
are the characteristic laser parameters (i.e. the laser power, pulse duration
and pulse repetition rate), characteristic properties of host materials (e.g.
their chromatic dispersions, optical nonlinearities or self-focusing features,
etc.) and the relative importance of physical processes such as the avalanche
impact ionization, multiphoton ionization and electron-hole radiative
recombination processes. These factors act in conjunction to impose the regime
of laser operation, in particular their competition determines the appropriate
laser operation regime. In this work a theoretical study is proposed to
explore the effects of the competition between multiphoton absorption, plasma
ionization and electron-hole radiative recombination processes, on the laser
dynamics in transparent materials with Kerr nonlinearity. The study rests on a
model consisting of a K-order nonlinear complex Ginzburg-Landau equation,
coupled to a first-order equation describing time variation of the electron
plasma density. An analysis of stability of continuous waves, following the
modulational-instability approach, reveals that the combination of multiphoton
absorption and electron-hole radiative recombination processes can be
detrimental or favorable to continuous-wave operation, depending on the group-
velocity dispersion of the host medium. Numerical simulations of the model
equations in the full nonlinear regime, reveal the eixstence of pulse trains
the amplitudes of which are enhanced by the radiative recombination processes.
Numerical results for the density of the induced electron plasma feature two
distinct regimes of time evolution, depending on the strength of the electron-
hole radiative recombination processes.

  * Resubmission date:

January 2021

Keywords: Laser-matter interactions, Continuous waves, Nonlinear transparent
media, Pulses, Plasma ionization

##  1 Introduction

Femtosecond laser micromachining nowadays offers the most reliable and
portable tool in a broad range of modern industrial material processing [1, 2,
3, 4], its applications extend from accurate manufacturing of electronic
devices, to fine drilling and machining of hard metals, ceramics and soft
plastics into various micro textures for improvement of functions and
properties of end products [2, 3, 4, 5]. In these applications laser pulses
focused on a dielectric medium are absorbed via nonlinear photo-ionization
mechanism [6, 7], leading to a permanent modification of material structure at
scales on the order of nanometers. In the specific context of transparent
materials [6, 7, 8], at low pulse powers the modification will be a smooth
refractive index change which can be exploited advantageously in the
fabrication of photonic devices [2, 9]. However, at higher pulse powers the
modification gives rise to more complex processes such as birefringences,
periodic nanoplanes aligning themselves orthogonally to the laser polarization
to form periodic nanogratings, change in the electronic structure due to
electron and hole productions from charge ionization with the generation of
electron plasma, electron-hole radiative recombination processes [2, 6, 9],
etc..

In accordance with the fineness required for the end product, optical fields
used in laser material processing can be grouped in two categories [2, 9, 10]:
continuous-wave (CW) lasers, which usually extend up to several kilowatts, and
pulsed lasers with an average power spanning well below one kilowatt thus
providing a wide range of wavelengths and pulse duration, as well as pulse
repetition rates [10]. Due to these attributes, pulse lasers can allow
micromachining with high resolution in depth and therefore offer a rich
potential for applications in drilling [2], cutting [11], welding [12, 13],
ablation [12], material surface texturing and scripting [12, 13]. We remark
that besides their short duration and high powers, femtosecond pulse lasers
have been most attractive owing to their minimal thermal drawbacks [14, 15].
Indeed femtosecond lasers are able to accumulate heat such as to minimize
defect-induced damages, as a matter of fact this heat accumulation prevents
undesired physical casualties as for instance the formation of microcracks,
material bending, etc. [16] notably during laser processing involving
transparent materials [17, 18, 19].

Theoretical investigations of femtosecond laser processing on transparent
materials have attracted a little attention in the past [7, 20, 21]. Yet
theoretical studies are useful for they provide fundamental knowledge relative
to the global picture of the laser dynamics, relevant for a good understanding
of its distinct possible operation regimes in the prospect of an optimization
of the micromachining technology. Indeed, since femtosecond lasers are optical
fields with duration far below picosecond, they belong to a specific class of
lasers known as ultrashort lasers [22]. Lasers in this specific class operate
typically in pulsed modes of relatively high powers, nevertheless in some
contexts they can be tailored to operate in the CW regime [2, 9]. This is for
instance the case when their input powers are below the typical power of a
high-intensity optical pulse, or when the input field is of low power and is
designed to grow upon propagation from CW mode to a high-intensity pulse. Such
growth can be regarded as an instability-induced dynamical transition of the
CW laser, this instability is best known as modulational-instability process
[23, 24, 25, 26] and involves CW breakup into high-power optical fields,
leading ultimately to a pulse via a regime dominated by pulse-train structures
[27, 28, 29, 30].

It follows that for a good understanding of laser processing we must have a
clear picture of possible dynamical regimes of the laser. In nonlinear optical
materials this issue can be understood in terms of laser self-starting
dynamics, where the input is assumed to be a CW field whose amplitude grows
upon propagation until a threshold amplitude. Beyond this threshold amplitude
the CW mode becomes modulational unstable, typically this instability will
first generate weakly nonlinear pulse trains which decay subsequently into
high-intensity temporal pulses.

In this work we examine the dynamics of a model of femtosecond laser
interacting with Kerr nonlinear media, taking into consideration material
modification and the generation of an electron plasma. In refs. [7, 20, 21],
mathematical models were proposed to describe femtosecond laser interaction
with Kerr nonlinear optical media. These models are generally composed of a
complex Ginzburg–Landau (CGL) equation for laser propagation, coupled to a
first-order time evolution equation for the electron plasma density.
Multiphoton ionization result in a K-order nonlinear term both in the CGL
equation and in the rate equation. In these previous works [7, 20, 21] the
authors discussed only thermal aspects of laser interactions with Kerr-type
transparent materials, in this respect they established that multiphoton
ionization processes would drastically minimize thermal drawbacks.

The model considered in the present study emanes from two previous models,
namely the model proposed in ref. [21] where account was taken of the
contribution of avalanche ionization due to K-photon absorption processes to
the electron plasma generation, and the model of ref. [20] where electron-
recombination processes were taken into consideration. In our case we extend
the model to the physical context in which electron-hole radiative
recombination processes are relevant. This translates into an additional term
in the time evolution equation for the plasma density, quadratic in the
density variable and competing with the avalanche impact ionization and the
multi-photon ionization processes.
To start we analyze the system dynamics in the CW regime. In this respect we
find steady-state solutions, and explore their stability by means of the
modulational-instability theory. A global stability picture will be proposed
in terms of a two-dimensional complex parameter space, mapped by the real and
imaginary parts of the coefficient of spatial amplification of noise over a
finite range of values of the modulation frequency. Next we discuss the system
dynamics in the full nonlinear regime, starting with fixed-point solutions and
then carrying out numerical simulations, using a sixth-order Runge-Kutta
algorithm adopted from Luther [32], to generate nonlinear time series of the
laser amplitude and the electron plasma density, for different values of the
radiative recombination coefficient and of the photon number $K$.

##  2 The model, CW solutions and modulational instability

Consider an optical field propagating along the $z$ axis of a transparent
medium with Kerr nonlinearity. We assume that the energy stored by the
propagating laser modifies the electronic structure of the bulk material,
creating an electron plasma of variable density and multiphoton ionization
processes. Instructively the problem was addressed experimentally [7, 10, 20,
21] in the contexts of femtosecond laser pulses focused on silica materials.
It was established that for a strong focusing geometry the femtosecond laser
can cause bulk damages in the optical material, followed by a narrow track
with submicron width indicating a filamentary propagation of the pulse laser.

In refs. [7, 10, 20, 21] numerical simulations were carried out to investigate
the distribution of the energy stored in the bulk material by the propagating
laser. While the laser equation was exactly similar in all these previous
studies [7, 10, 20, 21], the time-evolution equations for the electron plasma
density were different. The difference resided in the nature of physical
processes assumed to contribute to the electron plasma ionization in each of
these studies, as a matter of fact the work of ref. [7] was mainly interested
in the impact of a linear recombination process on the time evolution of the
electron plasma density, whereas in refs. [10, 20, 21] the linerar
recombination process was neglected in favor of radiative recombinations only
[10], avalanche impact only [20] or a combination of avalanche impact
ionization and linear recombination processes [21]. However in none of these
studies a detailed analysis of the system dynamics was considered.

In the present study we pay interest to the system dynamics, considering a
model that combines the effects of avalanche ionization process, the electron-
hole radiative recombination and multiphoton ionization processes. The model
can be represented by a propagation equation for the laser field given in
terms of the cubic complex Ginzburg-Landau equation with a $K$-order nonlinear
term [7, 10, 20, 21], coupled to a first-order ordinary differential equation
accounting for time evolution of the electron plasma density, i.e.:

$\displaystyle i\frac{\partial u}{\partial z}$ $\displaystyle=$
$\displaystyle\delta\frac{\partial^{2}u}{\partial
t^{2}}-\sigma|u|^{2}u-i\gamma(1-i\omega_{0}\tau_{0})\rho u-i\mu|u|^{2K-2}u,$
(1) $\displaystyle\frac{\partial\rho}{\partial t}$ $\displaystyle=$
$\displaystyle\nu|u|^{2}\rho+\alpha|u|^{2K}-a\rho^{2}.$ (2)

Instructively, we assumed a paraxial approximation for the field propagation
considering the propagation axis to be $z$, and ignored beam diffractions in
the plane transverse to the propagation axis. The first term in the right-hand
side of eq. (1) accounts for the group-velocity dispersion, the second term
accounts for the intrinsic (i.e. Kerr) nonlinearity of the host material, the
third term describes both the plasma absorption and laser defocusing, and the
last term describes multiphoton absorption processes with $\mu$ the $K$-photon
ionization rate.

Eq. (2) describes the temporal evolution of the electron plasma density
$\rho$, we remark that this equation differs from those considered in refs.
[7, 20, 21] where the same problem was discussed. Indeed eq. (2) takes into
account all together the contributions from avalanche impact ionization (first
term), multiphoton ionization (second term) and electron-hole radiative
recombination (last term). This remark also holds for the study carried out in
ref. [31], where electron-hole radiative recombinations were not taken into
consideration. Parameters in eq. (1) and eq. (2) have the following physical
meanings and characteristic units (where necessary) [7, 10, 20, 21]:

  * •

$\delta$ (in fs${}^{2}/m$) is the coefficient of group-velocity dispersion,

  * •

$\sigma$ (in $m^{2}/W$) is the Kerr nonlinearity coefficient,

  * •

$\omega_{0}$ (in MHz) and $\tau_{0}$ (in fs) are the frequency and
characteristic relaxation time of the electron plasma respectively,

  * •

$\gamma$ (in $m^{2}$) is the cross section for inverse Bremsstrahlung,

  * •

$\nu$ (in $m^{2}/eV$) is the coefficient of avalanche impact ionization,

  * •

$\alpha$ is the multiphoton ionization coefficient in the dense medium,

  * •

and $a$ (in $m^{3}/s$) is the electron-hole radiative recombination
coefficient.

Typical values for these parameters are found in most experimental works
dealing with the problem, as for instance in the study of femtosecond laser
filamentations in transpartent media [10], laser micro-modification of fused
silica [7] and so on. However, for a theoretical study such as the present
one, experimental values of these parameters are of no useful given that the
mathematical model, represented by eqs. (1) and (2), involves normalized
variables such as the propagation time $t$, the propagation distance $z$, the
field amplitude $u$ as well as some coefficients in these two equations [10,
7]. Therefore in our study we shall select arbitrary but reasonable values for
characteristic parameters of the model, but keeping track of their appropriate
signs.

Being two nonlinear equations, general solutions to the coupled set eqs.
(1)-(2) are nonlinear waves. Nevertheless, provided specific conditions,
linear solutions including harmonic waves and CWs can also exist for the same
set. Thus steady-state CW solutions to eqs. (1)-(2) can be expressed:

$u(z)=\sqrt{I_{p}}\exp\left(iP_{c}z\right),\quad\rho=\rho_{0},$ (3)

which upon substitution in (1)-(2) give:

$\displaystyle P_{c}$ $\displaystyle=$ $\displaystyle\sigma
I_{p}+\gamma\omega_{0}\tau_{0}\rho_{0},\qquad
I_{p}=\left(\frac{\sigma}{\mu\omega_{0}\tau_{0}}\right)^{\frac{1}{K-2}},$ (4)
$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\frac{\nu}{2a}I_{p}\left[1-\sqrt{1+\frac{4a\alpha\sigma}{\nu^{2}\mu\omega_{0}\tau_{0}}}\right].$
(5)

Here $P_{c}$, the CW wave number, is fixed by the input power $I_{p}=|u|^{2}$
as well as the equilibrium value $\rho_{0}$ of the electron plasma density
$\rho$. Given that these two last quantities (i.e. $I_{p}$ and $\rho_{0}$)
depend on characteristic parameters of the model, they cannot be arbitrary and
hence can be tuned by varying characteristic parameters of the model.

To look into the stability of the CW solutions eqs. (3), we carry out a
modulational-instability analysis assuming a small-amplitude perturbation
$f(z,t)$ of the CW amplitude $\sqrt{I_{p}}$. Steady-state solutions to eqs.
(1)-(2) thus become:

$u(z,t)=\left[u_{0}+f(z,t)\right]\exp(iP_{c}z),\qquad\rho(t)=\rho_{0}+\delta\rho(t).$
(6)

Replacing in eqs. (1)-(2) and linearizing, we obtain:

$\displaystyle i\frac{\partial f}{\partial z}$ $\displaystyle-$
$\displaystyle\delta\frac{\partial^{2}f}{\partial
t^{2}}+C^{(1)}_{K}\,(f+f^{*})=-i\gamma(1-i\omega_{0}\tau_{0})u_{0}\delta\rho(t),$
(7) $\displaystyle\frac{\partial\delta\rho(t)}{\partial t}$ $\displaystyle-$
$\displaystyle q_{0}\delta\rho(t)=C^{(2)}_{K}\,(f+f^{*}),$ (8)

Where,

$\displaystyle C^{(1)}_{K}$ $\displaystyle=$ $\displaystyle\sigma
u^{2}_{0}+i\mu(K-1)u^{2K-2}_{0},\quad
C^{(2)}_{K}=\frac{1}{\gamma}(\alpha\gamma K-\mu\nu)u^{2K-1}_{0},$
$\displaystyle q_{0}$ $\displaystyle=$ $\displaystyle\nu
u^{2}_{0}+\frac{2a\mu}{\gamma}u^{2K-2}_{0},$ (9)

$f^{*}$ in the linear equations (7) and (8) denotes the complex conjugate of
$f$. The first-order inhomogeneous linear equation (8) can be solved by means
of Green’s function technique [25, 31], yielding:

$\displaystyle\delta\rho(t)=C^{(2)}_{K}\int_{-\infty}^{t}(f+f^{*})e^{-q_{0}(t^{\prime}-t)}dt^{\prime}.$
(10)

Because of the presence of $f^{*}$ in eq. (7), we must consider its complex
conjugate. With this consideration, as solution to eq. (7) we pick
$f(z,t)=A_{1}\exp\left(\kappa z+i\Omega t\right)$ and
$f^{*}(z,t)=A_{2}\exp\left(\kappa z+i\Omega t\right)$, where $\kappa$ is the
coefficient of spatial amplification of the perturbation and $\Omega$ is the
modulation frequency for the perturbation. These solutions lead to the
following matrix-form eigenvalue problem:

$\displaystyle\kappa\left(\begin{array}[]{c}A_{1}\\\ A_{2}\end{array}\right)$
$\displaystyle=$ $\displaystyle\left[\left(\begin{array}[]{cc}M_{1}&M_{2}\\\
M^{*}_{2}&M^{*}_{1}\end{array}\right)-S_{0}\left(\begin{array}[]{cc}1&1\\\
1&1\end{array}\right)\right]\left(\begin{array}[]{cc}A_{1}\\\
A_{2}\end{array}\right)$ (19) $\displaystyle-$ $\displaystyle
T_{0}\left(\begin{array}[]{cc}N_{1}&N_{1}\\\
N^{*}_{1}&N^{*}_{1}\end{array}\right)\left(\begin{array}[]{c}A_{1}\\\
A_{2}\end{array}\right),$ (24)

where:

$\displaystyle M_{1}$ $\displaystyle=$ $\displaystyle
i(\delta\Omega^{2}+\sigma u_{0}^{2}),\quad M_{2}=i\sigma u_{0}^{2},$ (25)
$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle\mu(K-1)u_{0}^{2K-2},\quad N_{1}=1-i\omega_{0}\tau_{0},$ (26)
$\displaystyle T_{0}$ $\displaystyle=$
$\displaystyle(K\alpha\gamma-\mu\nu)u^{2K}_{0}\left(i\Omega-\nu
u^{2}_{0}-2a\mu u^{2K-2}_{0}/\gamma\right)^{-1}.$ (27)

The two possible eigenvalues of the above $2\times 2$ matrix equation are:

$\displaystyle\kappa_{1,2}$ $\displaystyle=$
$\displaystyle-\mu(K-1)u^{2K-2}_{0}\left[1+\frac{\gamma(K\alpha\gamma-\mu\nu)u_{0}^{2}}{\mu(K-1)(i\Omega\gamma-\nu
u_{0}^{2}\gamma-2a\mu u_{0}^{2K-2})}\right]$ (28) $\displaystyle\pm$
$\displaystyle\sqrt{(S_{0}+T_{0})^{2}-(\delta\Omega^{2}+\sigma
u_{0}^{2})^{2}+\sigma^{2}u_{0}^{4}-2T_{0}\delta\Omega^{2}\omega_{0}\tau_{0}},$

where subscripts $1,2$ refer to the plus (+) and minus (-) signs respectively.
From the standpoint of modulational-instability theory, CW stability will
depend on the sign and nature (i.e. real or complex) of the two eigenvalues
$\kappa_{1,2}$. Most generally the following possible situations are expected:

  * •

When the real part of $\kappa$ is zero, the CW solution will be always stable
irrespective of the sign of its imaginary part.

  * •

When the real part of $\kappa$ is negative, the CW solution will be
asymptotically stable (i.e. is stabilized after some roundtrips) irrespective
of the sign of its imaginary part.

  * •

When the real part of $\kappa$ is positive, the CW regime will be always
unstable.

Given that the two eigenvalues are functions of the modulation frequency
$\Omega$, we find it more appropriate to first consider the CW stability at
zero modulation frequency. In this later case the eigenvalues are:

$\kappa_{1}=0,\qquad\kappa_{2}=-2\mu(K-1)u^{2K-2}_{0}+\frac{2\gamma(K\alpha\gamma-\mu\nu)u_{0}^{2K}}{\nu
u_{0}^{2}\gamma+2a\mu u^{2K-2}_{0}}.$ (29)

It turns out that laser self-starting (i.e. CW instability) will be favored
provided $\kappa_{2}>0$, or in terms of formula (29);

$a<\frac{K\gamma(\alpha\gamma-\mu\nu)}{2\mu^{2}(K-1)I^{K-2}_{p}}.$ (30)

Quantitatively, this condition implies two possible characteristic values of
the radiative recombination coefficient $a$ above which laser self-starting
can occur: One is negative for $\alpha\gamma<\mu\nu$ and hence is nonphysical,
whereas the positive and physical one is conditioned by $\alpha\gamma>\mu\nu$
and is:

$\displaystyle
a_{th}=\frac{K\gamma(\alpha\gamma-\mu\nu)}{2\mu^{2}(K-1)I^{K-2}_{p}}.$ (31)

In concrete terms the quantity $a_{th}$ sets a threshold value of the
electron-ion recombination coefficient, above which the laser will self-start.
Due to the strong dependence of $\kappa_{1,2}$ in eq. (28) on the modulation
frequency, discussing CW stability from the analytical expressions of
$\kappa_{1,2}$ for arbitrary nonzero values of $\Omega$ is far from being an
easy task. Therefore we resort to a global analysis, by mapping the two
eigenvalues onto a plane Re($\kappa$)-Im($\kappa$) describing a two-
dimensional complex parameter space, where Re($\kappa$) and Im($\kappa$) are
real and imaginary parts respectively of the eigenvalue $\kappa$. In this
parametric representation, the modulation frequency $\Omega$ plays the role of
a parameter and so can span a broad range of values, which in our case will be
the finite interval $-5\leq\Omega\leq 5$. The first figures we consider are
parametric representations of Im($\kappa$) as a function of Re($\kappa$), for
some selected combinations of values of key characteristic parameters of the
model. To be more explicit, the four graphs in figs. 1 and 2 represent
Im($\kappa$) as a function of Re($\kappa$) in the anomalous dispersion regime
($\delta<0$: fig. 1) and normal dispersion regime ($\delta>0$: fig. 2)
respectively, for $K=2$, $3$, $4$ and $5$. Values of model parameters are
given in the captions, and different curves in each graph correspond to
different values of the radiative recombination coefficient $a$. Recall that
the sign of $\delta$ determines the dispersion regime [24, 25, 34], indeed a
positive $\delta$ corresponds to a normal group-velocity dispersion, whereas a
negative $\delta$ will correspond to an anomalous group-velocity dispersion
well known [24, 25, 34] to favor the generation of pulse structures, of course
provided the intrinsic refractive index of the host medium is of a self-
focusing Kerr nonlinearity.

Figure 1: (Color online) Imaginary versus real parts of $\kappa_{1}$ (full
curve) and $\kappa_{2}$ (dashed curve) for $K=2,3,4,5$. The radiative
recombination coefficient $a$ is varied as $a$= $0$, $0.001$, $0.002$,
$0.003$, $0.004$. $\alpha=0.6$ ,$\nu=0.5$, $\mu=0.1$, $I_{p}=2.5$,
$\omega_{0}\tau_{0}=0.2$, $\sigma=0.8$, $\delta=-0.5$, $\gamma=0.1$.

Figure 2: (Color online) Imaginary versus real parts of $\kappa_{1}$ (full
curve) and $\kappa_{2}$ (dashed curve) for $K=2,3,4,5$. The radiative
recombination coefficient $a$ is varied as $a$= $0$, $0.001$, $0.002$,
$0.003$, $0.004$. $\alpha=0.6$ ,$\nu=0.5$, $\mu=0.1$, $I_{p}=2.5$,
$\omega_{0}\tau_{0}=0.2$, $\sigma=0.8$, $\delta=0.5$, $\gamma=0.1$.

Fig. 1 suggests that small values of $K$ are expected to favor laser self-
starting in the anomalous dispersion regime. As $K$ increases, the real part
of the largest eigenvalue gradually shifts to the negative branch and
consequently, the CW regime is stabilized. On the contrary fig. 2 indicates
that in the normal dispersion regime, CW operation will be favored for small
values of the multiphoton absorption rate $K$. As $K$ increases, CW modes
become unstable thus favoring laser self-starting.

##  3 Nonlinear regime

In the full nonlinear regime, solutions to the laser equation eq. (1) are
high-intensity fields which can be represented as real-amplitude pulses,
undergoing spatio-temporal modulations i.e. [33, 34]:

$u(z,t)=g(z,t)\exp i\left[\phi(z,t)-\omega z\right],$ (32)

where $g$ is the real amplitude and $\phi$ is the modulation phase. We
introduce a reduced time as $\tau=t-vz$, in which $v$ is the pulse inverse
velocity such that $\omega$ in formula (32) emerges more explicitly as a
nonlinear shift in the propagation constant. Letting $g(z,t)\equiv g(\tau)$
and $\phi(z,t)\equiv\phi(\tau)$, inserting eq. (32) into eq. (1) and eq. (2)
and separating real parts from imaginary parts, we obtain the following set of
coupled first-order nonlinear ordinary differential equations:

$\displaystyle\left(vM+\omega+\delta
M^{2}+\gamma\omega_{0}\tau_{0}\rho\right)g-\delta y^{\prime}+\sigma g^{3}$
$\displaystyle=$ $\displaystyle 0,$ (33) $\displaystyle\left(\delta
M^{\prime}-\gamma\rho\right)g+\left(v+2\delta M\right)y-\mu g^{2K-1}$
$\displaystyle=$ $\displaystyle 0,$ (34) $\displaystyle\rho^{\prime}-\nu
g^{2}\rho-\alpha g^{2K}+a\rho^{2}$ $\displaystyle=$ $\displaystyle 0,$ (35)

with $y=g^{\prime}$, $M=\phi^{\prime}$ and the prime symbol refers to
derivative with respect to $\tau$. Let us focus on the system dynamics in the
particular case $v=0$ [26, 33, 34]. For this value of $v$ the set of coupled
first-order nonlinear ordinary differential equations (33)-(35) reduces to:

$\displaystyle M^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\gamma\rho}{\delta}-\frac{2My}{g}+\frac{\mu
g^{2K-2}}{\delta},$ $\displaystyle y^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\left(\omega+\delta
M^{2}+\gamma\omega_{0}\tau_{0}\rho\right)g}{\delta}+\frac{\sigma
g^{3}}{\delta},$ $\displaystyle g^{\prime}$ $\displaystyle=$ $\displaystyle
y,$ $\displaystyle\rho^{\prime}$ $\displaystyle=$ $\displaystyle\nu
g^{2}\rho+\alpha g^{2K}-a\rho^{2}.$ (36)

Our first interest will be on the singular solutions to this system, which are
their fixed points, with the aim to probe the effects of important
characteristic parameters of the model such as the radiative recombination
coefficient $a$ and the multiphoton absorption rate $K$, on equilibrium
solutions of the laser amplitude $g$ and instantaneous frequency $M$, as well
as of the electron plasma density $\rho$.

###  3.1 Fixed points

Singular solutions to the set of first-order nonlinear ordinary differential
equations (36), which are their fixed points, are the roots of the following
nonlinear system:

$\displaystyle\frac{\gamma\rho}{\delta}+\frac{\mu g^{2K-2}}{\delta}$
$\displaystyle=$ $\displaystyle 0,$ (37)
$\displaystyle\frac{\left(\omega+\delta
M^{2}+\gamma\omega_{0}\tau_{0}\rho\right)g}{\delta}+\frac{\sigma
g^{3}}{\delta}$ $\displaystyle=$ $\displaystyle 0,$ (38) $\displaystyle y$
$\displaystyle=$ $\displaystyle 0,$ (39) $\displaystyle\nu g^{2}\rho+\alpha
g^{2K}-a\rho^{2}$ $\displaystyle=$ $\displaystyle 0.$ (40)

From this system we derive:

$M^{2}=\frac{1}{\delta}\left[\mu\omega_{0}\tau_{0}g^{2K-2}-\omega-\sigma
g^{2}\right],\qquad\rho=g^{K}\,\sqrt{\left(\frac{\alpha\gamma-\mu\nu}{a\gamma}\right)}.$
(41)

Remarkable enough formula (41) suggests that irrespective of the value of $K$,
the electron plasma density $\rho$ will be zero when the laser amplitude $g$
is zero.

Fig. 3 represents the variations of the two fixed points of the laser
amplitude $g$ as a function of $\omega$, for $K=2,3,4,5$ and values of model
parameters given in the captions. It should be noted that by fixed points of
$g$ we understand its extrema, i.e. its maximum and minimum which are obtained
by annihilating $M$ in formula (41).

Figure 3: (Color online) Fixed points of the laser amplitude $g$ as a function
of $\omega$, for $K=2,3,4,5$, $\sigma=0.8$, $\omega_{0}\tau_{0}=0.2$,
$\mu=0.25$.

According to fig. 3, for $K=2$ the laser dynamics is dominated by weakly
nonlinear pulse trains with maximum amplitudes for $\omega=0$. An increase of
$K$ (see graphs for $K=3$, $4$ and $5$) favors strongly nonlinear pulse trains
of larger amplitudes [33, 34]. We have also plotted the instantaneous
frequency $M$ as a function of the amplitude $g$ of laser (fig. 4), and the
laser amplitude $g$ as function of the electron plasma density $\rho$ (fig.
5), for four different values of $K$. Remark that the expression of $M$ given
in formula (41) does not contains the radiative recombination coefficient $a$
but is controlled mainly by the laser propagation constant $\omega$, whereas
$g$ as a function of $\rho$ extracted from eq. (41) depends on $a$ but not on
$\omega$. Therefore, we have chosen to plot $M$ as a function of $g$ by
considering two cases i.e., the case $\omega=0$ and the case of finite nonzero
value of $\omega$ as one sees in fig. 4.

Figure 4: (Color online) Variation of the instantaneous frequency $M$ with the
amplitude $g$ of laser, for $K=2,3,4,5$ and two different values of $\omega$
(Solid curve correspond to $\omega=0$, and dashed curve to $\omega=0.5$).
Values of other parameters are $\sigma=0.8$, $\omega_{0}\tau_{0}=0.2$,
$\delta=-0.5$ , $\mu=0.25$.

The different curves in the graphs of fig. 4 show that $M$ is enhanced by an
increase of $g$, and that there is a threshold value of the amplitude beyond
which the instantaneous frequency is expected to decrease to zero. As it is
apparent, this threshold value of $g$ (and consequently of $M$) is decreased
with an increase of $K$.

Figure 5: (Color online) Variation of the laser amplitude $g$ with the
electron plasma density $\rho$, for $K=2,3,4,5$ and different values of the
radiative recombination coefficient $a$ indicated in the graphs. $\nu=0.1$
$\alpha=0.6$, $\gamma=0.1$ , $\mu=0.25$.

In fig. 5, the fixed point of the electron plasma density is always zero at
zero amplitude of the laser whatever the value of $K$, consistently with what
we learned from formula (41). However there is a drastic increase of $\rho$
with an increase of $K$, meaning that for the same laser intensity the
stronger the multiphoton absorption processes the larger will be the electron
plasma density.

###  3.2 Pulse structures

Considering the full nonlinear dynamics of the system, the set of first-order
ordinary differential equations (36) was solved numerically using a sixth-
order Runge-Kutta algorithm adapted from ref.. [32]. Because eqs. (1)-(2)
involve several parameters, all of which cannot be varied in this study, we
fixed most parameters except two ones i.e. the multiphoton absorption rate
$K$, which was given the four different values $K=2,3,4,5$, and the radiative
recombination coefficient $a$ which was varied in three distinct ranges of
values where three distinct behaviors were noticed.

Figs. 6 and 7 are time variations of the laser amplitude $g(t)$ and of the
electron plasma density $\rho(t)$, for four distinct values of $K$ and values
of model parameters listed in the figure captions.

Figure 6: (Color online) Time variation of the laser amplitude $g$, for
$K=2,3,4,5$ and different values of the radiative recombination coefficient
$a$. $\nu=0.1$ $\alpha=0.6$, $\gamma=0.1$ , $\mu=0.25$. Values of $a$ are,
from the smallest to the largest amplitudes: 0, 3.9, 5.6, 6.

Figure 7: (Color online) Time variation of the electron plasma density $\rho$,
for $K=2,3,4,5$ and values of $a$ indicated in the graphs. $\sigma=0.8$,
$\omega_{0}\tau_{0}=0.8$, $\delta=-0.5$ , $\mu=0.5$, $\nu=0.5$, $\gamma=0.18$,
$\alpha=0.6$.

Note that for each graph, we plotted dynamical quantities for different values
of $a$ starting with $a=0$ (i.e. when there is no radiative recombination).
Our objective in so doing was to highlight the qualitative and quantitative
influences of $a$, on the system dynamics.

In fig. 6, the laser amplitude $g(t)$ is manifestly a pulse train with a
maximum varying only weakly with $a$ for $K=2$. However, as $K$ is increased
the maximum of $g$ gets more and more large. Note the fall-off of $g(t)$ with
time, which is more and more pronounced as $K$ gets larger and for $a=0$. We
attribute this fall-off to a damping effect induced by the density of electron
plasma i.e. $\rho(t)$, which in this context acts like a laser gain/loss.

The electron plasma density $\rho$, plotted in fig. 7, is increasing with time
for $a=0$ irrespective of the value of $K$. However, when $a$ is increased in
a relatively large range of values, $\rho(t)$ decreases exponentially in time
tending to its equilibrium value with an increasingly sharp slope. Evidently,
the nonzero values of $a$ chosen for the numerical results just discussed, are
large enough and therefore do not enable one appreciate how the electron
plasma density $\rho$ changes from its exponentially increasing feature for
$a=0$, to an exponentially decreasing feature when $a\neq 0$. In order to earn
more insight onto the effects of $a$ on the time variation of $\rho$, in fig.
8 and fig. 9 we plot $\rho$ versus time keeping the same values of parameters
used in fig. 7, but choosing smaller values of $a$.

Figure 8: (Color online) Time variation of the electron plasma density $\rho$,
for $K=2,3,4,5$ and values of $a$ indicated in the graphs. $\sigma=0.8$,
$\omega_{0}\tau_{0}=0.8$, $\delta=-0.5$ , $\mu=0.5$, $\nu=0.5$, $\gamma=0.18$,
$\alpha=0.6$ .

Figure 9: (Color online) Time variation of the electron plasma density $\rho$,
for $K=2,3,4,5$ and values of $a$ indicated in the graphs. $\sigma=0.8$,
$\omega_{0}\tau_{0}=0.8$, $\delta=-0.5$, $\mu=0.5$, $\nu=0.5$, $\gamma=0.18$,
$\alpha=0.6$ .

To be explicit, in fig. 8 values of $a$ were chosen small but large enough
such that $\rho(t)$ is already decreasing with time. On the contrary the
values chosen for $a$ in fig. 9, are closer to zero. The physically relevant
insight is that in the range of values of $a$ for which $\rho$ decreases, an
increase of $a$ sharpens the slope of the decrease of $\rho$. When $a$ is very
small such that $\rho$ is an exponentially increasing function of time, an
increase of $a$ will soften the variation of $\rho$. It is quite apparent that
in the very small range of values of $a$, $\rho$ is periodically oscillating
as it increases exponentially with time (fig. 9). The quantity $a_{th}$ in the
graphs of fig. 8, is the characteristic value of the radiative recombination
coefficient for which the exponential variation of $\rho$ with time is
suppressed. For this value of $a$, the electron plasma density $\rho(t)$
constantly oscillates between two positive extrema and hence is always
nonzero, in average.

##  4 Conclusion

We have investigated the dynamics of a model for femtosecond laser inscription
in a transparent medium with Kerr nonlinearity. This model is an improvement
of two recent ones [20, 21], in which the contribution of electron-hole
radiative recombination processes to the plasma generation was neglected.
Mathematically the model is represented by a complex Ginzburg-Landau equation
with cubic nonlinearity plus a $K$-order nonlinearity accounting for $K$
photon absorption processes, coupled to a time first-order ordinary
differential equation with a term quadratic in the electron plasma density,
accounting for radiative recombination processes.

In view of the fact that the laser dynamics in the host material determines
the specific regime of operation of the laser during laser micromachining on
the material, we considered two distinct regimes of laser dynamics namely CWs
and pulses. First addressing CWs and their stability, a linear-stability
analysis was carried out following the modulational-instability theory. From
this analysis we constructed a global stability map, which enabled us explore
parameter regimes in which CW operations could be stable.

In the nonlinear regime, an ansatz was introduced which aimed at representing
the laser as a pulse field with a real amplitude and real phase. With the help
of this ansatz the system dynamics was transformed into a set of four first-
order ordinary differential equations. Numerical solutions to these equations
established the pulse-train structure of the laser amplitude, the power of
which was increased by an increase of the radiative recombination coefficient
$a$. The electron plasma density was found to increase or decrease
exponentially in time with an oscillating amplitude, depending on the order of
magnitude of $a$.

Though the model treated in this work is rich, it remains limited for it
concerns only few out of a wide range of possible physical contexts. Indeed
not all transparent materials oppose a nonlinear optical response of the Kerr
type to field propagation, in fact the Kerr response is usually considered as
a very weak nonlinearity. There exists transparent materials with stronger
optical nonlinear responses, in general their refractive index is a series of
powers of the optical field intensity extending beyond the quadratic term.
Usually this is also represented by saturable-nonlinearity functions [25]. Few
examples of such materials include doped silica fiber glasses, where the
doping with rare-earth ions provides means of controlling the strength of
nonlinearity of the transparent host. A study similar to the present one, but
taking into consideration the saturable nonlinearity of the host material,
will surely enrich current knowledge on the fundamental physics of CW and
pulse operations in femtosecond laser micromachining involving transparent
materials with non-Kerr optical nonlinearity [25].

A. M. Dikandé thanks the Ministry of Higher Education of Cameroon (MINESUP)
for financial assistance, within the framework of the presidential grant
designated ”Research Modernization Allowances”.

## References

## References

  * [1] F. Dausinger, F. Lichtner and H. Lubatschowski, Femtosecond Technology for Technical and Medical Applications (Springer, Berlin, 2004).
  * [2] G. Chryssolouris, Laser machining: theory and practice ($1^{rst}$ edt., Springer, New York, 1991).
  * [3] N. B. Dahotre and A. Samant, Laser Machining of Advanced Materials (CRC Press, Taylor and Francis, 2014).
  * [4] J. Krüger and W. Kautek, Laser Phys. 9, 30 (1999).
  * [5] M. Stolze, T. Herrmann and J. A. L’huillier, Laser Applications in Microelectronic and Optoelectronic Manufacturing (LAMOM) XXI, Proceedings of SPIE 9735, 97350O, March 14th (2016).
  * [6] R. R. Gattass and E. Mazur, Nature Photon. 2, April (2008).
  * [7] A. V. Dostovalov, A. A. Wolf, V. K. Mezentsev, A. G. Okhrimchuk and S. A. Babin, Opt. Expr. 23, 32541 (2015).
  * [8] A. G. Okhrimchuk, V. K. Mezentsev, H. Schmitz, M. Dubov and I. Bennion, Laser Phys. 19, 1415 (2009).
  * [9] S. M. Eaton, G. Cerullo and R. Osellame, Fundamentals of Femtosecond Laser Modification of Bulk Dielectrics, in: ”Femtosecond Laser Micromachining”, Topics in Applied Physics (R. Osellame, G. Cerullo and R. Ramponi, eds.) vol 123. (Springer, Berlin, 2012).
  * [10] A. Couairon and A. Mysyrowicz, Phys. Rep. 441, 47 (2007).
  * [11] H. Shin and D. Kim, Opt. Laser Technol. 102, 1 (2018).
  * [12] P. K. Shukla, J. Lawrence and Y. Zhang, Opt. Laser technol. 75, 40 (2015).
  * [13] Zh. Wang,Y. B. Li, F. Bai, C. W. Wang and Q. Z. Zhao, Opt. Laser Technol. 81, 60 (2016).
  * [14] J. W. Jung and C. M. Lee, Cutting temperature and laser beam temperature effects on cutting tool deformation in laser-assisted machining, in: Proceedings of IMECS II, pp. 1817–1822 (2009).
  * [15] H. B. Xu, J. Hu, H. Sheng and Z. C. Du, J. Phys. Conf. Ser. 633, 0120901 (2015).
  * [16] S. M. Eaton, H. Zhang and P. R. Herman, Opt. Expr. 13, 4708 (2005).
  * [17] D. Savastru, R. Savastru, S. Miclos and I. I. Lancranjan, Opt. Laser Eng. 110, 288 (2018).
  * [18] K. Minoshima, A. W. Kowalevicz, I. Hartl, E. P. Ippen and J. G. Fujimoto, Opt. Lett. 26, 1516 (2001).
  * [19] S. Nolte, S., M. Will, J. Burghoff and A. Tünnermann, J. Mod. Opt. 51, 2533 (2004).
  * [20] L. Sudrie, A. Couairon, M. Franco, B. Lamouroux, B. Prade, S. Tzortzakis and A. Mysyrowicz, Phys. Rev. Lett. 89, 186601 (2002).
  * [21] J. S. Petrovic, V. Mezentsev, H. Schmitz and I. Bennion, Opt. Quant. Electron. 39, 939 (2007).
  * [22] M. Wollenhaupt, A. Assion and T. Baumert, Short and Ultrashort Laser Pulses, in: ”Springer Handbook of Lasers and Optics” (Springer, Berlin, 2012).
  * [23] T. B. Benjamin and J. E. Feir, J. Fluid Mech. 27, 417 (1967).
  * [24] C. J. Chen, P. K. A. Wai and C. R. Menyuk, Opt. Lett. 20, 350 (1995).
  * [25] A. M. Dikandé, J. Voma Titafan and B. Z. Essimbi, J. Opt. 19, 105504 (2017).
  * [26] E. M. Ntongwe and A. M. Dikandé, Opt. Quant. Electron. 51, 361 (2019).
  * [27] A. M. Dikandé, Phys. Rev. A81, 013821 (2010).
  * [28] D. Jr. Jubgang Fandio and A. M. Dikandé, J. Opt. Soc. Am. B34, 2721 (2017).
  * [29] R. D. Dikandé Bitha and A. M. Dikandé, Phys. Rev. A97, 033813 (2018).
  * [30] R. D. Dikandé Bitha and A. M. Dikandé, Eur. Phys. Jour. D73, 152 (2019).
  * [31] P. Kameni Nteutse, A. M. Dikandé and S. Zekeng, Opt. Quantum Elec. 52, 313 (2019).
  * [32] H. Luther, Math. Comp. 22 434 (1968).
  * [33] J. M. Soto-Crespo, N. Akhmediev and G. Town, J. Opt. Soc. Am. 19, 234 (2002).
  * [34] F. G. Mbieda Ngomegni, A. M .Dikandé and B. Z. Essimbi, Phys. Scrip. 95, 025502 (2020).
  * [35] J. Lawrence (ed.), Advances in Laser Materials Processing: Technology, Research and Application (second ed., Elsevier, Amsterdam, 2018).