Dynamical analysis of mass-spring models using Lie algebraic methods

TL;DR

This paper applies Lie algebraic methods to analyze the dynamics of mass-spring vibrational systems, providing solutions for various configurations including finite arrays, bridging classical and quantum-inspired techniques.

Contribution

It introduces Lie algebraic techniques to solve complex vibrational systems, extending classical methods with quantum-inspired mathematical tools.

Findings

Solutions for finite circular and linear arrays using classical methods.

Application of Lie algebraic methods to more complex arrays.

Unified framework connecting classical and quantum approaches.

Abstract

The dynamical analysis of vibrational systems of masses interconnected by restitution elements each with a single degree of freedom, and different configurations between masses and spring constants, is presented. Finite circular and linear arrays are studied using classical arguments, and their proper solution is given using methods often found in quantum optical systems. We further study some more complicated arrays where the solutions are given by using Lie algebras.

Equations

$${\ddot{q}}_0+k(2q_0-q_N-q_1)$$

$${\ddot{q}}_j+k(2q_j-q_{j+1}-q_{j-1})$$

$${\ddot{q}}_N+k(2q_N-q_0-q_{N-1})$$

$$\ddot{\mathbf{Q}}=\mathbb{M}\mathbf{Q},$$

$${\mathbb{M}}_{\mathrm{c}}=k\left(\begin{array}{cccccc}-2& 1& 0& 0& \cdots & 1\\ 1& -2& 1& 0& \cdots & 0\\ 0& 1& -2& 1& \cdots & 0\\ 0& 0& 1& -2& \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & \ddots & 1\\ 1& 0& 0& \cdots & 1& -2\end{array}\right),$$

$$\mathbf{Q}={\left(\begin{array}{ccccc}{q}_0& q_1& q_2& \cdots & q_N\end{array}\right)}^T,$$

$$\mathbf{Q}\left(t\right)=\cosh (t\sqrt{{\mathbb{M}}_{\mathrm{c}}})\mathbf{Q}(0),$$

$$\mathbb{V}=\left(\begin{array}{cccccc}0& 1& 0& 0& \cdots & 0\\ 0& 0& 1& 0& \cdots & 0\\ 0& 0& 0& 1& \cdots & 0\\ 0& 0& 0& 0& \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & \ddots & 1\\ 1& 0& 0& \cdots & 0& 0\end{array}\right).$$

$$\mathbb{F}=\frac{1}{\sqrt{N+1}}\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ {\lambda }_0& {\lambda }_1& {\lambda }_2& \cdots & {\lambda }_N\\ {\lambda }_0^2& {\lambda }_1^2& {\lambda }_2^2& \cdots & {\lambda }_N^2\\ {\lambda }_0^3& {\lambda }_1^3& {\lambda }_2^3& \cdots & {\lambda }_N^3\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\lambda }_0^N& {\lambda }_1^N& {\lambda }_2^N& \cdots & {\lambda }_N^N\end{array}\right)$$

$${\lambda }_j=\exp \left(\mathrm{i}\frac{2\pi }{N+1}j\right),j=0,1,2,...,N$$

$$\Lambda =\left(\begin{array}{ccccc}{\lambda }_0& 0& 0& \cdots & 0\\ 0& {\lambda }_1& 0& \cdots & 0\\ 0& 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & {\lambda }_N\end{array}\right)$$

$$\mathbb{Q}\left(t\right)$$

$$=\mathbb{F}\cosh \left[t\sqrt{k(\Lambda +{\Lambda }^{\ast })-2k\mathbb{I}}\right]{\mathbb{F}}^{†}\mathbf{Q}(0),$$

$$\mathbb{F}:=\frac{1}{\sqrt{N+1}}\sum _{n,m=0}^N\exp \left[\mathrm{i}\frac{2\pi mn}{N+1}\right]|m⟩⟨n|,$$

$$\cosh$$

$$=\sum _{m=0}^N\cos \left[2t\sqrt{k}\sin \left(\frac{\pi m}{N+1}\right)\right]|m⟩⟨m|.$$

$$\mathbf{Q}(t)=$$

$$\times \cos \left[2t\sqrt{k}\sin \left(\frac{\pi m}{N+1}\right)\right]|n⟩,$$

$$q_n(t)=$$

$$\times \cos \left[2t\sqrt{k}\sin \left(\frac{\pi m}{N+1}\right)\right].$$

$$m_0{\ddot{q}}_0+\left(k_0+k_1\right)q_0-k_1{q}_1$$

$$m_j{\ddot{q}}_j+\left(k_j+k_{j+1}\right)q_j-k_{j+1}{q}_{j+1}-k_j{q}_{j-1}$$

$$m_N{\ddot{q}}_N+\left(k_N+k_{N+1}\right)q_N-k_N{q}_{N-1}$$

$$\ddot{\mathbf{Q}}(t)={\mathbb{M}}_{\mathrm{L}}\mathbf{Q}(t),$$

$${\mathbb{M}}_{\mathrm{L}}$$

$$+\sum _{j=0}^{N-1}\frac{k_{j+1}}{m_{j+1}}|j+1⟩⟨j|+\sum _{j=0}^{N-1}\frac{k_{j+1}}{m_j}|j⟩⟨j+1|$$

$$\mathbf{Q}(t)=\sum _{j=0}^N{q}_j\left(t\right)|j⟩.$$

$$\mathbf{Q}(t)=\cosh \left(t\sqrt{{\mathbb{M}}_{\mathrm{L}}}\right)\mathbf{Q}\left(0\right).$$

$${\mathbb{M}}_{\mathrm{L}}=-2\sum _{j=0}^N|j⟩⟨j|+\sum _{j=0}^{N-1}(|j+1⟩⟨j|+|j⟩⟨j+1|),$$

$${\ddot{q}}_0+k(2q_0-q_1)$$

$${\ddot{q}}_j+k(2q_j-q_{j+1}-q_{j-1})$$

$${\ddot{q}}_N+k(2q_N-q_{N-1})$$

$${\mathbb{M}}_{\mathrm{L}}=\mathbb{SD}{\mathbb{S}}^{-1},$$

$$\mathbb{S}=\sum _{i,j=0}^N\frac{U_i(y_{j+1})}{\sqrt{{\sum }_{s=0}^N[U_s(y_{j+1}){]}^2}}|i⟩⟨j|,$$

$$\mathbb{D}=2k\sum _{j=0}^N\left[\cos ({\varphi }_{j+1})-1\right]|j⟩⟨j|.$$

$$\mathbf{Q}(t)=\mathbb{S}\cosh (t\sqrt{\mathbb{D}}){\mathbb{S}}^{-1}\mathbf{Q}(0).$$

Full text

modal analysis, particle kinematics, quantum optics

I Introduction

The dynamical description of systems with coupled subelements is a well known subject in classical mechanics literature Meirovitch1986 ; landauM ; Kotkin1980 . It is a typical task to take the stated system of dynamical coupled equations, that emerge from the equilibrium analysis, and solve them with some differential equations or eigenvalue techniques. A simple model that describes a chain of classical particles (atoms) harmonically coupled with their nearest neighbors and subjected to a periodic on-site (substrate) potential has become in recent years one of the fundamental and universal models of low dimensional nonlinear physics. In spite of the fact that a link with the classical model is not often stated explicitly in many applications, many kind of nonlinear problems involving the dynamics of discrete nonlinear chains are in fact based on the classical formulation introduced in the papers by Ya. Frenkel and T. Kontorova YaFrenkel1938 ; T.A.Kontorova1938 ; T.A.Kontorova1938a ; Braun1998 ; Allen1998 , who suggested to use this kind of nonlinear chain to describe, in the simplest way, the structural dynamics of a crystal lattice in the vicinity of the dislocation core Braun1998 ; Allen1998 .

Furthermore, the study of these particular systems, that are intrinsically classical due to their macroscopic mechanical nature, can be linked to the study of a propagated light field under a waveguide array, which in turn can also be described by classical means. Some work has been made on photonic lattices to study the analogy between quantum systems and classical light propagation Makris2006 ; Keil2011 ; PerezLeija2010 . Under this scheme, the system of equations that describes the time evolution of the oscillation amplitude can be written in terms of operators that fulfill certain known commutation relations from quantum mechanics. In this sense, the system of equations of a tight- binding model to first neighbors and periodic boundary conditions can be written using the discrete Fourier transform PerezLeija2016 . On the other hand, linear finite arrays in both paradigms, light or mechanical, exhibit a similitude when the same methods to obtain solutions and insights are applied. The emulation of quantum mechanical properties with classical propagated light is well known Makris2006 ; PrezLeija2012 , following here that the nearest neighbor interaction of mechanical systems mimics some of the features encountered in the classical light counterpart.

The goal of this manuscript is to take a few well-known classical physical systems and fully solve them with techniques found often in the description and solution of quantum mechanical problems. In Section II, we start with the analysis of a finite set of masses connected to each other by springs with constants of harmonic restitution following Hook’s law, constrained to move with a single degree of freedom in a circle. The set of differential equations that describe the system has a periodic boundary condition, which allows the first and last element to be coupled. In Section III, we restate the problem removing the boundary condition of the circular array and keeping only a finite chain fixed at both or one of the edges. In this class of arrays we distinguish when all the restitution coefficients and masses are equal III.1, which leads to an analytical solution given by Chebyshev polynomials of the second kind. The rise of traveling normal waves is observed. If the restitution coefficients follows some other law in function of the position of the masses, then other phenomena are present. We engineer an interaction matrix III.2 where his diagonalization is given in terms of Kravchuk functions that, in turn, are solutions of the discrete and finite harmonic oscillator of $\mathfrak{su}(2)$ Atakishiyev2008 ; ata . Here, the solution is also oscillatory but with persistence of the poles and nodes presented by the Kravchuk functions that are in the core hypergeometrical functions. We observe bouncing of amplitudes that become rapidly a complete interference patterns when the propagation time is sufficiently large, giving no recombination nor recovering of initial conditions. We conclude with the approach and solution of a problem that relates to the masses and springs through binomial coefficients III.3. Here, although the analytical solution exists and belongs too to the realm of the well-known algebra $\mathfrak{su}(2)$, but it is somewhat not easy to calculate explicitly, and instead numerical results are presented.

II Circular finite array

Let us consider the interaction between a finite set of masses labeled by $m_{j},\;j=0,1,2,...,N$, where the dimension of the set is $d_{N}=N+1$. Geometrically arranged in a circle, the interaction of $m_{j}$’s is mediated by springs with equal restitution constants $k$, as it is shown schematically in Fig. 1. The position of every single mass is labeled by the canonical coordinate $q_{j}(t),\;j=0,1,2,...,N$. Taking $m_{j}=1$ for every $j$ in the set, the coupled equations of motion of this system is given by Meirovitch1986

[TABLE]

where for the sake of simplicity, we drop the time dependence of the space coordinates, $q_{j}(t)\equiv q_{j}$.

It is straightforward to cast the set of differential equations (II) onto the matrix form

[TABLE]

where $\mathbb{M}$ is the tridiagonal real matrix plus bounded corners, with dimensions $d_{N}^{2}$, explicitly given by

[TABLE]

being $\mathbf{Q}$ the $d_{N}$ column vector of canonical coordinates $q_{j}$’s

[TABLE]

where clearly $q_{j}$ is time dependent, so $\mathbf{Q}\equiv\mathbf{Q}(t)$.

Because we are dealing with a set of second order differential equations in (II), we need to establish initial conditions, one set in the positions $\mathbf{Q}(0)$, and other set in the velocities $\dot{\mathbf{Q}}(0)$. A careful looking onto the matrix form (3) suggest an ansatz for the initial value problem (2); mainly, it needs to be a continuous differentiable real function with no parity associated, thus we propose it to be

[TABLE]

which can be verified as a legal solution by direct substitution in Eq. (2).

Now is time to restructure the problem (2) in order to use some of the methods encountered in quantum optics. First, we notice that the coefficient matrix (3) has a tridiagonal form and the corners occupied with ones, which suggests to use the well known London operators $\{\mathbb{V},\mathbb{V}^{\dagger}\}$, whose matricial representation is given by London1926 ; PerezLeija2016

[TABLE]

It can be shown that the matrix $\mathbb{V}$ obeys the spectral decomposition $\mathbb{V}={\mathbb{F}}\Lambda{\mathbb{F}}^{\dagger}$ moya2011differential , where $\mathbb{F}$ is the discrete Fourier transform given by the Vandermonde confluent matrix PerezLeija2016 ; MoyaCessa2018

[TABLE]

with

[TABLE]

the $j$th root of the unity, and

[TABLE]

as the eigenvalue diagonal matrix.

Using the former equations (6), (7) and (8), we can substitute in (5) and simplify the proposed solution as

[TABLE]

where it is worth to notice that the terms inside the square root of (II) are a sum of pure diagonal matrices.

Taking the basis vectors to be a Kronecker basis of dimension $d_{N}$ of the form $\ket{j}:=\delta_{n,j}$ for $n=0,1,2,\dots,N$; that is, defined $j$, the ket $\ket{j}$ is a vector that has zeros everywhere, except when $n=j$; when the Hilbert space dimension is infinite, these states are known as a Fock states in the quantum mechanical realm. These kets can be used to define the Fourier kernel (7) as an operator $\mathbb{F}$ in terms of the outer product of their elements as PerezLeija2016 ,

[TABLE]

for which can be easily verified the inverse property $\mathbb{F}^{\dagger}=\mathbb{F}^{-1}$. Also the hyperbolic cosine in (9) can be expressed, with the help of Eqs. (8) and (7), as

[TABLE]

With this restructured problem, we need only to establish the initial conditions on the positions and leave the initial velocities unknown. The initial conditions read $\mathbf{Q}(0)=w\ket{l}$ in terms of the orthonormal basis $\ket{l},\;l\in\{0,1,2,\cdots,N\}$, and it is a vector in a finite Hilbert space of dimension $N+1$. This means that the $l$-th mass in the circle is excited with an excitation weight $w$; in other words, an amplitude of initial perturbation. We are now in place to forward the proposed solution (5) onto a closed analytical form, with the help of (II), (11), (II) and the initial condition. Remembering that $\braket{r}{s}=\delta_{r,s}$, we arrive to

[TABLE]

which means that for every single solution in the vector $\mathbf{Q}(t)$, we obtain the solely dynamic dictated for every mass in the array by

[TABLE]

In Figure 2, we plot the time evolution (II) for $N=30$ and an initial condition $\mathbf{Q}(0)=\ket{1}$. As can be seen, the general mechanical behavior is oscillatory. It is important to notice, that because of the periodic boundaries in the array, part of the initial amplitude is transmitted from the $q_{j}$ mass to the neighbor $q_{j+N}$; in this case some of the initial amplitude in $q_{0}$ at time $t_{0}$ is transmitted to $q_{30}$ at a time $t+\Delta t$, explaining that the coupled evolution plot have interference at some time $t_{i}$. It is equivalent to say, because the circular arrange of the masses, that the initial condition on $\ket{0}$ when the system evolves, looks like two initial conditions, one at $t_{0}$ and the other, retarded, at time $t+\Delta t$.

III Linear finite array of $N+1$ masses

We now consider the task to determine the dynamic evolution of a linear finite array of masses $m_{j},\;j=0,1,2,...,N$, with restitution elements $k_{j},\;j=0,1,2,3,...,N$, as seen above in Section II, and which schematic representation is shown in Fig. 3. In this case, the system of ordinary differential equations that governs the evolution is given by Meirovitch1986

[TABLE]

This system can be packed into the matrix representation

[TABLE]

where $\mathbb{M}_{L}$ is represented in the finite basis $\ket{l}$ as

[TABLE]

and

[TABLE]

Using the ansatz previously given, the proper solution of (16) can be stated as

[TABLE]

It is clear that the real difficulty lies on two facts of the matrix function $\cosh\left(t\sqrt{\mathbb{M}_{\mathrm{L}}}\right)$; first, we need to evaluate the matrix $\mathbb{M}_{L}$ at time $t$; second, this evaluation need to operate onto the initial condition on the right. The last issue can be simple tackle in the case when $\mathbb{M}_{L}$ can be diagonalized, this lets us obtain a simple form of the solution. The task now is to examine some cases where the matrix could be diagonalized in function of the nature of $m_{j}$ and $k_{j}$.

III.1 $m_{j}$ and $k_{j}$ equal to one for all $j$’s

It is straightforward to take all masses and restitution elements as one; that is $m_{j}=1$ and $k_{j}=1$ for all $j\in{0,1,2,\cdots N}$ as they exhibit no dependence across the array positions nor the time parameter. Rewriting (III) with the new conditions, we obtain

[TABLE]

where, in contrast with the matrix in (3), the matrix in (20) lacks of the ones at the corners, which maintain joined the array into a circle. The linear finite array implies that the chain of masses could have not material or physical constrictions on the borders; that is, just contour conditions to keep the system well defined. The set of coupled differential equations governing the evolution of this system is Meirovitch1986

[TABLE]

To solve this problem, we follow the presentation in rmf57.2 proposing the spectral decomposition of the interaction matrix as

[TABLE]

where the matrix operator $\mathbb{S}$ is defined as

[TABLE]

being $U_{i}(x)$ the Chebyshev polynomials of second kind Abramowitz1964 ; Olver2010 ; $y_{j}=\cos(\phi_{j})$ with $\phi_{j}=\frac{j\pi}{N+2},\;j=0,1,2,...,N$ as the $N+1$ roots of the polynomials. The matrix $\mathbb{D}$ is diagonal with their elements given by Abramowitz1964 ; Olver2010

[TABLE]

Then, in analogy with Sec. II, the solution of the matrix equation $\ddot{\mathbf{Q}}\left(t\right)=\mathbb{M}_{\mathrm{L}}$ is find in the factorization

[TABLE]

Given the specific initial condition $\mathbf{Q}(0)=w\ket{l}$, with $l$ the number of the mass that is displaced with an amplitude $w$ respect to the equilibrium position, we arrive to the explicit form of the solution for every single mass $n$

[TABLE]

where it is worth to notice that the denominator is a sum of $N+1$ squared Chebyshev polynomials, that can be identified as a normalization dependence on the position of every single mass.

In Figure 4, we plot the temporal evolution (26) for an initial condition $\mathbf{Q}(0)=\ket{15}$ with $N=30$; that is, $N+1$ elements in the array. A centralized initial condition is elected, because in the former system of equations (III) the boundary conditions forbids the interaction between the first and final elements of the array, so no amplitude is transmitted from $q_{0}$ to $q_{N}$. Setting the initial condition in the center of the array permits us observe that the propagation follows a symmetric evolution until it arrives to the edges, rebounds and recombine with the amplitudes coming from the center of the array, giving place to interference at times $t>20$. Of course, the points where the propagation meets the edges is a function of the pair values $\{m_{j},k_{j}\}$, because of our unit value election, the presented behavior follows.

III.2 Kravchuk interaction

We can do a step forward and engineer an iteration matrix $\mathbb{M}_{L}$ which diagonalization is given in terms of hypergeometrical functions, more specific, discrete orthogonal polynomials. Following Regniers Regniers2009 , we can propose the interaction matrix to be of the form

[TABLE]

where the coefficients follow the laws

[TABLE]

It can be proved that $p=1/2$ gives the diagonalization

[TABLE]

where $\mathbb{D}=-\mathrm{diag}(0,1,\cdots,N)$ and the matrix elements of $\mathbb{U}$ are defined by

[TABLE]

with $w(j)=\binom{N}{j}p^{j}(1-p)^{N-j}$, and $h_{i}=\left(\frac{1-p}{p}\right)^{i}/\binom{N}{i}$. The functions $k_{i}(j)$ are the symmetric Kravchuk polynomials, whose use is extensive in the description of the discrete and finite harmonic oscillator ata ; Weimann2016 . These polynomial are defined in terms of the Gaussian hypergeometric function as

[TABLE]

thus, the election of $p=1/2$ is justified for the obtainment of the symmetric functions that fulfill the requirements dictated by the spectral theorem (29).

In order to solve the problem

[TABLE]

we can write the matrices $\mathbb{U}$ and $\mathbb{D}$, in the discrete basis $\ket{l}$, as

[TABLE]

As $\mathbb{D}$ is diagonal, we can obtain the proper factorization of $\mathbf{Q}(t)$ as in the previous cases in the form

[TABLE]

and then, with the definitions of the matrices (33) and the initial condition $\mathbf{Q}(0)=w\ket{l}$, we obtain the vector solution as

[TABLE]

Finally projecting over the $m$th element in the array, we obtain the single dynamical function for the position

[TABLE]

for $m\in\{0,1,\cdots,N\}$ and $l$ the index of the initial condition, that is, which mass is excited at time $t=0$.

This last exercise was somehow easy because we can diagonalize $\mathbb{M}_{L}$ in terms of a pure diagonal matrix $\mathbb{D}$ and the Kravchuk matrices $\mathbb{U}$. It is important to remark that the functions (30) are the solutions of the discretization of the quantum harmonic oscillator embedded in the compact algebra $\mathfrak{su}(2)$. This functions are a feasible approximation of the Hermite-Gauss functions in the discrete and finite space of the algebra.

In Figure 5, we plot the temporal evolution of (36), here we set a dimension $N=16$ and an initial condition with unit amplitude at $\ket{0}$. We see that the slope of the propagation is nearly acute, giving that the initial excitation transmitted to the adjacent masses reach the edge at $q_{N}$ very rapidly, where a bouncing is observed; after that, the propagation follows another acute slope but with some recombination of amplitudes occurring. After the second bounce we observe not a regular pattern, just recombination and interference of amplitudes. This last behavior let us assume that the nearest neighbor interaction (27) presents some complex features after some time in the evolution (36).

III.3 Case $m_{j}=\binom{N}{N-j}$, $k_{j}=jm_{j}$

In this section, the masses $m_{j}$ and the restitution elements $k_{j}$ will be defined in terms of the binomial coefficients. We propose to use

[TABLE]

in such way that both quantities follow some sort of binomial distribution. It is important to notice that $m_{j}$ is well defined for all the values of $j$, but that the restitution elements has a null value when $j=0$. This last issue is not a problem nor a slip in the statements, it just says that the first restitution element is present with a intrinsic value of restitution equal to zero, following that the mathematical constraint is well defined. Physically this is the same picture of Figure 3, but the left and right material attachments are leave free.

Doing the calculations for the coefficients in (III), we arrive to the definition of $\mathbb{M}_{L}$ as

[TABLE]

The diagonal is constant with value $-N$ and the off-diagonals has up and downward values, from $1$ to $N$ and vice versa, respectively.

The coefficients previously showed can be cast to the operational representation as

[TABLE]

in such way that (III.3) can be rewritten as

[TABLE]

The operator given in (III.3) obeys the commutation relations

[TABLE]

thus, they are a representation of the Lie algebra $\mathfrak{su}(2)$. Using the commutation relations, we transforms the operators (III.3) in such way that $\mathbb{M}_{L}$ is diagonal. It is long, but straightforward, to arrive to the set of transformations

[TABLE]

Hence, the interactions matrix can be written as

[TABLE]

This last equation can be carry to a full diagonal form if we choose the parameters to be $\alpha=-1$ and $\beta=1/2$,

[TABLE]

where $\Lambda=\sum_{j=0}^{N}j\ket{j}\bra{j}$. Inverting (46), we finally obtain that

[TABLE]

The diagonal representation of $\mathbb{M}_{L}$ given above, let us now look for the solution of the initial value problem (19) as

[TABLE]

Expression (48) is an analytic closed expression; however, their explicit calculation is cumbersome. Therefore, instead of giving the long-complete expression, we accelerate the process using the matrix representation (39) for a fixed $N$ and make a numeric evaluation for the initial condition $\mathbf{Q}(0)=w\ket{l}$. Also due to the definition of the masses and restitution elements, when we excite masses near $j=0$ the dynamics is not interesting, because around this specific place the restitution constant vanish. So, to obtain relevant results in the dynamics, we propose to use initial conditions around the middle mass. Figure 6 shows a numerical evaluation of (48) when the initial condition is a coherent state around the center of the array.

IV Conclusions

In this work we present a set of mechanical systems that are solved using mathematical methods and arguments often encountered in the analysis and solution of quantum optical phenomena. We found that the time evolution of the position amplitude of a chain of masses has some resemblance with light propagation in graded indexed waveguide arrays; we may conclude that in some particular cases (specific values of spring constants and masses) there is an isomorphism between both systems. This conclusion may be relevant because it motivates the search of relations between mechanical systems of coupled oscillators and systems that obey Schrödinger and Helmholtz-like equations Unpublished . Finally, we may say that the use of these quantum optical methods are fully equivalent to those that follow Hamiltonian or Lagrangian developments for many-body interactions. The simplicity of the solution arise from the fact that the interaction matrix is fully diagonalizable.

V Acknowledgments

Alejandro R.U. acknowledge CONACyT for their financial support in the development of this work through Ph.D. grant #449192. I.R.-P. thanks Prof. J. Récamier for his hospitality at ICF-UNAM and acknowledge partial support from DGAPA UNAM project PAPIIT IN111119, and Beca de Colaboración INAOE.