1 Introduction↩︎

The Jaynes–Cummings model (JCM), first introduced in 1963 by Edwin Jaynes and Fred Cummings [1], is one of the paradigmatic models of quantum light–matter interaction. It describes the coherent exchange of excitations between a two-level atom and a single quantized field mode and has played a central role in quantum optics, cavity quantum electrodynamics, and quantum information science [2]–[9]. Its interaction has no direct classical counterpart in the sense that the energy exchange occurs in discrete quanta, while the eigenstates of the coupled Hamiltonian form the well-known dressed-state basis, which is generally entangled across the atom and field degrees of freedom [10]. These features make the JCM a particularly suitable setting in which to investigate how mathematical structures such as symmetries and conservation laws reflect the underlying quantum dynamics.

The conventional analytical treatment of the JCM is based on the decomposition of the dynamics into invariant excitation subspaces and on the diagonalization of the Hamiltonian within each of them. Although this construction provides the standard dressed-state solution, the dynamics can also be formulated as a coupled system of partial differential equations for the atomic components and the field mode in the characteristic-function representation [11], [12]. This formulation retains the full atom–field structure while providing a natural framework for studying the model from the perspective of differential equations.

In this work, we explore what Lie symmetry theory reveals about this PDE representation of the JCM. Symmetry analysis identifies continuous transformation groups admitted by the dynamical equations and uses them to construct invariant solutions, reduce the number of independent variables, and expose structural relations that are not evident from the Hamiltonian diagonalization alone [13]–[15]. In this sense, the objective is not to replace the standard solution of the JCM, but to characterize it from a symmetry-based perspective and to determine whether the admitted symmetries produce additional classes of physically admissible solutions.

The group-analysis procedure is naturally semi-inverse. One first derives the solution forms imposed by a selected symmetry and subsequently determines the initial and boundary conditions compatible with them. In the present problem, physical admissibility is fixed by the quantum-mechanical conditions imposed on the characteristic function and on the reduced atomic density matrix. This distinction is essential: the differential equations may admit mathematically valid invariant solutions that do not represent physical quantum states. The symmetry analysis therefore provides both a classification of invariant solutions and a criterion for determining which of them survive the physical constraints of the model.

The application of this method recovers the familiar dressed-state dynamics as one of the invariant solutions of the PDE system. At the same time, another subalgebra of group generators leads to a distinct class of invariant solutions whose radial dependence is expressed in terms of Heun polynomials. These solutions extend beyond the conventional dressed-state construction and are compatible with generic coupled atom–field configurations, although their complete physical interpretation remains to be established. Their appearance illustrates one of the main outcomes of the symmetry approach: it can reveal admissible analytical structures that are not naturally suggested by the usual Hamiltonian treatment.

The differential formulation also provides a natural setting for the analysis of non-trivial conservation laws. For this purpose, we employ the methodology developed by Vinogradov’s school, which relates conserved currents of a system of differential equations to their generating functions  [16]–[18]. Within this framework, the standard conservation of the total number of excitations is recovered directly from the PDE system. In addition, further conserved currents are obtained whose densities and fluxes involve atomic populations, coherence, purity, and characteristic-function moments associated with the field dynamics.

Of particular interest is a conservation-law balance equation for the quantity \(P_a-f_a^2\), where \(P_a\) is the purity of the reduced atomic state and \(f_a\) is its coherence function. Its time variation is determined by the atom–field coupling strength and by combinations of atomic expectation values with derivatives of the field characteristic functions. This result does not imply that entanglement itself is directly conserved. Rather, it identifies a conserved current involving quantities whose evolution is linked to the redistribution of purity, coherence, and correlations between the atom and the field. For globally pure atom–field states, changes in local purity are directly associated with the development of bipartite entanglement, which gives these conservation laws a natural connection with the entanglement dynamics of the JCM. A related association between conserved quantities and entanglement dynamics has also been suggested in Ref. [19].

The aim of this work is therefore to examine the mathematical and physical outcomes of applying symmetry analysis and conservation-law theory to the JCM. The method provides a different characterization of known solutions, reveals an additional class of invariant solutions, recovers the standard excitation-number constant of motion, and produces non-trivial conserved currents associated with fundamental quantum properties of the atom–field dynamics. To the best of our knowledge, this is the first application of Vinogradov’s conservation-law framework to a composite quantum system, suggesting a possible route for the analysis of more complex quantum models formulated as systems of differential equations.

This paper is organized as follows. In Sec. 2, we derive the dynamical equations of the JCM in the characteristic-function representation and introduce the transformations used throughout the analysis. A brief review of the conventional Hamiltonian and dressed-state construction is included in Appendix 7. In Sec. 3, we determine the admitted point and generalized symmetries and derive the invariant solutions associated with the relevant subalgebras. In Sec. 4, we construct the corresponding conservation laws, including the excitation-number conservation law and the currents related to atomic purity, coherence, and atom–field correlations. Finally, in Sec. 5, we summarize the main results and their possible implications.

2 Differential equations of the JCM↩︎

In this section, we derive the system of differential equations governing the dynamics of the JCM in the characteristic-function representation. We start from the von Neumann equation \[\label{vNeu} \mathrm{i}{\mathrm{d}\over \mathrm{d}t}\varrho = [\hat{H},\varrho]\,,\tag{1}\] where \(\varrho\) is the density operator of the composite atom-field system and \(\hat{H}\) is the JCM Hamiltonian given by: \[\label{ham} \hat{H} = \frac{\Delta}{2}\sigma_z\otimes \openone + \openone\otimes\hat{n} + g\left( \sigma_{+}\otimes\hat{a} + \sigma_{-}\otimes\hat{a}^{\dagger} \right).\tag{2}\] Here, \(g=\lambda/(2\omega_o)\) is the dimensionless atom-field coupling strength, \(\sigma_+\) and \(\sigma_-\) are the atomic raising and lowering operators, and \(\hat{a}^{\dagger}\) and \(\hat{a}\) are the creation and annihilation operators of the field mode. The interaction term describes the coherent exchange of one excitation between the atom and the field: absorption of a field quantum raises the atom, whereas emission by the atom increases the field excitation number. This formulation allows one to describe both pure and mixed quantum states and is therefore more general than a state-vector description. Moreover, rather than working in the eigenbasis of the total Hamiltonian, we project Eq. (1 ) onto the bare atomic basis, \(\{|e\rangle\otimes\openone_f, |g\rangle\otimes\openone_f\}\). In this basis, \[\sigma_z = |e\rangle\langle e|-|g\rangle\langle g|,\quad \sigma_+ = |e\rangle\langle g|,\quad \sigma_- = |g\rangle\langle e|=\sigma_+^\dagger .\]

Writing the density operator in atomic blocks, \(\varrho_{ij}=\langle i|\varrho|j\rangle\), with \(i,j=e,g\), yields the following system of operator equations: \[\begin{align} \mathrm{i}\dot{\varrho}_{ee} &=& [\hat{n}+1/2,\varrho_{ee}] +g\left(\hat{a}\varrho_{ge}-\varrho_{eg}\hat{a}^{\dagger}\right), \nonumber\\[0.5em] \mathrm{i}\dot{\varrho}_{eg} &=& \Delta\varrho_{eg} +[\hat{n}+1/2,\varrho_{eg}] +g\left(\hat{a}\varrho_{gg}-\varrho_{ee}\hat{a}\right), \nonumber\\[0.5em] \mathrm{i}\dot{\varrho}_{ge} &=& -\Delta\varrho_{ge} +[\hat{n}+1/2,\varrho_{ge}] +g\left(\hat{a}^{\dagger}\varrho_{ee}-\varrho_{gg}\hat{a}^{\dagger}\right), \nonumber\\[0.5em] \mathrm{i}\dot{\varrho}_{gg} &=& [\hat{n}+1/2,\varrho_{gg}] +g\left(\hat{a}^{\dagger}\varrho_{eg}-\varrho_{ge}\hat{a}\right). \label{meq4} \end{align}\tag{3}\]

To apply Lie symmetry theory, we move from this operator formulation to the characteristic-function representation of the field degrees of freedom  [11]. This representation is obtained as the double Fourier transform of the Wigner function [20]. For each atomic block, we define \[\label{chf} w_{ij}(k,s,t) = \int \mathrm{d}q\,\mathrm{e}^{\mathrm{i}q k} \langle q+s/2|\varrho_{ij}(t)|q-s/2\rangle, \;i,j=e,g .\tag{4}\] Equivalently, the characteristic function of the composite system can be written as a \(2\times 2\) matrix in the atomic subspace, \[\label{dmatblocks} \mathrm{w}(\vec{r},t) = \left(\begin{array}{cc} w_{ee}(\vec{r},t) & w_{eg}(\vec{r},t)\\ w_{ge}(\vec{r},t) & w_{gg}(\vec{r},t) \end{array}\right),\tag{5}\] where \(\vec{r}=(k,s)\) (depending on the context, we will use a column vector or a row vector).

The characteristic-function description is particularly suitable for systems with continuous field variables, since it represents the dynamics in Fourier phase space. In this representation, the operator equations above become the following system of PDEs: \[\begin{align} \nonumber \mathrm{i}\hat{L} w_{ee} &=& \varkappa(s/2+\partial_s)(w_{ge}+w_{eg}) -\mathrm{i}\varkappa(k/2+\partial_k)(w_{ge}-w_{eg}),\\\nonumber &&\\\nonumber \mathrm{i}\hat{L} w_{eg} &=& \Delta w_{eg}\\\nonumber &&+{\varkappa\over 2}(s-\mathrm{i}k)(w_{gg}+w_{ee}) -\varkappa(\partial_s-\mathrm{i}\partial_k)(w_{ee}-w_{gg}),\\\nonumber &&\\ \mathrm{i}\hat{L} w_{ge} &=& -\Delta w_{ge} \label{wPDEs} \\ &&+{\varkappa\over 2}(s+\mathrm{i}k)(w_{gg}+w_{ee}) -\varkappa(\partial_s+\mathrm{i}\partial_k)(w_{ee}-w_{gg}),\nonumber \\ \nonumber &&\\\nonumber \mathrm{i}\hat{L} w_{gg} &=& \varkappa(s/2-\partial_s)(w_{ge}+w_{eg}) -\mathrm{i}\varkappa(k/2-\partial_k)(w_{ge}-w_{eg}),\\ \nonumber \end{align}\tag{6}\] where \(\partial_s := \partial/\partial s\), \(\partial_k := \partial/\partial k\), \(\varkappa=g/\sqrt{2}\) and \[\label{L-op} \hat{L}=\partial_t+s\partial_k-k\partial_s.\tag{7}\] The part \(s\partial_k-k\partial_s\) of operator \(\hat{L}\) generates rotations in the \((k,s)\) plane; consequently, the radius \(r=|\vec{r}|=\sqrt{k^2+s^2}\) is preserved along the corresponding flow. This rotational structure reflects, in Fourier phase space, the conservation of the total number of excitations in the JCM.

It is useful to introduce the combinations \[\begin{align} \nonumber w_f := w_{ee}+w_{gg}, &\qquad& w_x := w_{eg}+w_{ge},\\ w_y := {w_{ge}-w_{eg}\over \mathrm{i}}, &\qquad& w_z := w_{ee}-w_{gg}. \label{wfx-yz} \end{align}\tag{8}\] In terms of the vector \[\vec{w}:=(w_f,w_x,w_y,w_z) ,\] the system takes the compact form \[\label{afceqs} \hat{L}\vec{w}(\vec{r},t) = \mathrm{R}(\vec{r})\,\vec{w}(\vec{r},t),\tag{9}\] with \[\mathrm{R}(\vec{r}) = -\mathrm{i}\varkappa \left( \begin{array}{cccc} 0 & s & k & 0\\ s & 0 & -\mathrm{i}{\Delta/\varkappa} & -2\partial_s \\ k & \mathrm{i}{\Delta/\varkappa} & 0 & -2\partial_k \\ 0 & 2\partial_s & 2\partial_k & 0 \end{array} \right). \label{sys95w95f}\tag{10}\] The first component, \(w_f(\vec{r},t)\), corresponds to the characteristic function of the reduced field state, obtained after tracing over the atomic degrees of freedom. The remaining components are associated with the atomic observables; indeed, after tracing over the field degrees of freedom, one obtains \[\begin{align} \langle\sigma_x(t)\rangle &=& w_x(\vec{r},t)\big|_{\vec{r}=0},\\ \langle\sigma_y(t)\rangle &=& w_y(\vec{r},t)\big|_{\vec{r}=0},\\ \langle\sigma_z(t)\rangle &=& w_z(\vec{r},t)\big|_{\vec{r}=0}. \label{SIGMA95Z} \end{align}\tag{11}\]

For the subsequent analysis, we consider general atom-field initial conditions of the form \[\mathrm{w}^0(\vec{r}) := \mathrm{w}(\vec{r},0),\] subject to the normalization and Hermiticity conditions \[\label{INITIAL95W0} [\mathrm{w}^0(0)]_{11}+[\mathrm{w}^0(0)]_{22}=1,\qquad [\mathrm{w}^0(0)]_{12}=[\mathrm{w}^0(0)]_{21}^{\ast}.\tag{12}\] Moreover, \([\mathrm{w}^0(0)]_{11}\) and \([\mathrm{w}^0(0)]_{22}\) are real.

We will also consider the special case of initially factorized atom-field states, for which \[\mathrm{w}^0(\vec{r}) = \left(\begin{array}{cc} \varrho^0_{ee} & \varrho^0_{eg}\\ \varrho^0_{ge} & \varrho^0_{gg} \end{array}\right) w^0_f(\vec{r}). \label{IN95CON95W}\tag{13}\] Here, \(\varrho^0_{ij}:=\varrho_{ij}(0)=w^0_{ij}(0)\), \(w_f^0(\vec{r}):=w_f(\vec{r},0)\), while the reduced atomic density-matrix elements are given by: \[\varrho_{ij}(t):=w_{ij}(\vec{r},t)\big|_{\vec{r}=0}.\] The reduced atomic density-matrix elements must satisfy \[\begin{align} \nonumber \varrho_{ee}(t)+\varrho_{gg}(t)=1,\; \varrho_{ee}(t),\varrho_{gg}(t)\in\mathbb{R},\; \varrho_{eg}(t)=\varrho_{ge}^{\ast}(t).\\\label{atom95rho95cond} \end{align}\tag{14}\] In addition, the reduced field characteristic function \(w_f(\vec{r},t)\) must satisfy the quantum-mechanical admissibility conditions listed in Appendix 6.

Finally, because the operator \(\hat{L}\) has a natural rotational structure, it is convenient to introduce polar coordinates in Fourier phase space, \[s=r\cos\varphi,\qquad k=r\sin\varphi .\] In these variables, \[\label{L-op-polar} \hat{L}=\partial_t+\partial_\varphi .\tag{15}\] We then perform the transformation \[\label{rotpt} \vec{w}=\mathrm{T}(r,\varphi)\vec{u},\tag{16}\] where \[\mathrm{T}(r,\varphi) = \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & \mathrm{i}\cos(\varphi)/r & -\mathrm{i}\sin(\varphi)/r & 0 \\ 0 & \mathrm{i}\sin(\varphi)/r & \mathrm{i}\cos(\varphi)/r & 0 \\ 0 & 0 & 0 & 1 \end{array} \right),\] and \(\vec{u}=(u^1,u^2,u^3,u^4)\). The system for \(\vec{u}\) becomes \[\label{eqs} \hat{L}\vec{u}=\mathrm{S}(\vec{r})\vec{u},\tag{17}\] with \[\mathrm{S}(\vec{r}) = \left( \begin{array}{cccc} 0 & \varkappa& 0 & 0\\ -\varkappa r^2 & 0 & -\delta & 2\varkappa r\partial_r \\ 0 & \delta & 0 & 2\varkappa\partial_{\varphi} \\ 0 & {2\varkappa\over r}\partial_r & {2\varkappa\over r^2}\partial_{\varphi} & 0 \end{array} \right).\] This transformed system will be the starting point for the derivation of invariant solutions and conservation laws in the following sections.

3 Symmetries↩︎

In this section, we introduce the symmetry framework used to analyze the differential system obtained for the JCM. We follow the formulation in terms of generating functions of symmetries [16]–[18]. Let \[\label{PDEsystem} \vec{F}=\left(F^1,\ldots,F^m\right) =0\tag{18}\] be a determined regular system of PDEs of order \(d\), where each equation has the form \[F^\alpha = F^\alpha\left(\vec{x},\vec{u},\vec{u}^{(1)},\ldots,\vec{u}^{(d)}\right), \qquad \alpha=1,\ldots,m .\] Here, \(\vec{x}=(x^1,\ldots,x^n)\in\mathbb{R}^n\) denotes the independent variables, \(\vec{u}=(u^1,\ldots,u^m)\) denotes the dependent variables, and \(\vec{u}^{(s)}\) denotes the collection of all derivatives of \(\vec{u}\) of order \(s\), with \(1\leq s\leq d\).

The universal linearization operator [16], also known as the Fréchet derivative of the system [21], is the matrix differential operator obtained by linearizing \(\vec{F}\) with respect to the dependent variables and their derivatives, denoted by \[\label{UnivLinOp} l_{\vec{F}}=\left(\left[l_{\vec{F}}\right]_{\alpha\beta}\right), \quad \alpha,\beta=1,\dots,m.\tag{19}\] Its components are given by \[\label{l95F} \left[l_{\vec{F}}\right]_{\alpha \beta} = \sum_{0\leq |\varsigma|\leq d} {\partial F^\alpha \over \partial u_{\varsigma}^\beta} D^{(\varsigma)} .\tag{20}\] In this expression, \(\varsigma=(\varsigma^1,\ldots,\varsigma^n)\) is a multi-index, understood as an ordered set of non-negative integers, and \[|\varsigma|=\varsigma^1+\cdots+\varsigma^n,\quad D^{(\varsigma)} = D_1^{\varsigma^1}\circ\cdots\circ D_n^{\varsigma^n}.\] The notation \[u^\beta_\varsigma =\partial^{\varsigma^1}_{x^1}...\partial^{\varsigma^n}_{x^n} u^\beta\] denotes the derivative of \(u^\beta\) corresponding to the multi-index \(\varsigma\). The operator \(D_j\) is the total derivative operator with respect to \(x^j\), here expressed in expanded form as: \[\label{TOTAL95D} D_j = \partial_{x^j} + u^\alpha_j\partial_{u^\alpha} + u^\alpha_{ij}\partial_{u^\alpha_i} + u^\alpha_{ijk}\partial_{u^\alpha_{ik}} +\cdots ,\tag{21}\] with summation over repeated indices.

A function \[\vec{\phi}\left(\vec{x},\vec{u},\vec{u}^{(1)},\ldots\right) = (\phi^1,\ldots,\phi^m)\] is called the generating function, or characteristic [15], of a symmetry admitted by the system \(\vec{F}=0\), if it is a solution of the following system \[\label{l95F95phi} \left. l_{\vec{F}}\vec{\phi} \right|_{\mathcal{F}=0}=0.\tag{22}\] Here, \(\mathcal{F}=0\) denotes the infinite prolongation of the system \(\vec{F}=0\), that is, the system together with all its differential consequences \[D^{(\varsigma)}F^\alpha=0.\] The order of the generating function is determined by the highest order of derivatives of \(\vec{u}\) appearing in \(\vec{\phi}\).

In the present problem, the independent variables are \[x^1=t,\qquad x^2=r,\qquad x^3=\varphi,\] and the dependent variables are \[\vec{u}=(u^1,u^2,u^3,u^4) .\] The infinitesimal operator of symmetry \(X\) has the form \[X=\xi^i\partial_{x^i}+\eta^{\alpha}\partial_{u^{\alpha}}, \qquad i=1,2,3,\quad \alpha=1,\ldots,4. \label{X}\tag{23}\] The corresponding generating function is related to the coefficients of \(X\) by \[\label{phi95to95X} \phi^\alpha = \eta^\alpha - \xi^i u^\alpha_i.\tag{24}\] If the coefficients \(\xi^i\) and \(\eta^\alpha\) depend only on the variables \(x^i\) and \(u^\alpha\), there is a point symmetry. If they depend also on derivatives of \(\vec{u}\), the corresponding symmetry is a generalized symmetry.

For the JCM system written in the form (17 ), we now search for generating functions of first order, \[\label{phi} \vec{\phi}=\vec{\phi}\left(\vec{x},\vec{u},\vec{u}^{(1)}\right),\tag{25}\] satisfying the condition (22 ). For convenience, we introduce the radial variable \[\label{R} R=\frac{r^2}{4},\tag{26}\] which simplifies the radial part of the differential operators and any further reference to the capital letter \(R\) refers to this variable. In this variable, the universal linearization operator associated with the system (17 ) is the following \[\label{l95F95JM} l_{\vec{F}} = \begin{pmatrix} D_{t}+D_{\varphi} & -\varkappa& 0 & 0\\ 4\varkappa R & D_t+D_{\varphi} & \delta & -4\varkappa R D_R\\ 0 & - \delta & D_{t}+D_{\varphi} & -2\varkappa D_{\varphi}\\ 0 & - \varkappa D_R & -\frac{\varkappa}{2R}D_{\varphi} & D_t+D_{\varphi} \end{pmatrix}.\tag{27}\]

Solving Eq. (22 ) for first-order generating functions yields the following five independent generating functions: \[\begin{align} \tag{28} \vec{\phi}_1 &=& \left(\partial_t u^1,\partial_t u^2,\partial_t u^3,\partial_t u^4\right) , \\ \tag{29} \vec{\phi}_2 &=& \left(\partial_\varphi u^1,\partial_\varphi u^2, \partial_\varphi u^3,\partial_\varphi u^4\right) , \\ \tag{30} \vec{\phi}_3 &=& \left(u^1,u^2,u^3,u^4\right) , \\ \tag{31} \vec{\phi}_4 &=& \left(\tilde{u}^1,\tilde{u}^2,\tilde{u}^3,\tilde{u}^4\right) , \\ \tag{32} \vec{\phi}_5 &=& \left( \frac{\delta}{\varkappa}u^4+\frac{\partial_{\varphi}u^2}{2R} -\partial_R u^3, \, -2\partial_\varphi u^1, \right. \nonumber\\ && \left. 4R\left(-u^4+\partial_R u^1\right), \, \frac{\delta}{\varkappa}u^1-u^3 \right) . \end{align}\] Here, \(\vec{\tilde{u}}(t,r,\varphi)=(\tilde{u}^1,\tilde{u}^2,\tilde{u}^3,\tilde{u}^4)\) denotes an arbitrary solution of the same linear system (17 ).

The first four functions (28 )–(31 ) represent point symmetries (see the next subsection). The generating function \(\vec{\phi}_5\) does not correspond to a point symmetry. It is a first-order generalized symmetry with operator \(X^{(1)}=\phi_5^\alpha\partial_{u^\alpha}\), whose coefficients \(\eta^\alpha\) depend on the first derivatives of \(\vec{u}\).

Let us now consider two diagonal matrix differential operators: \[\mathrm{D}_t:=\operatorname{diag}(D_t,D_t,D_t,D_t), \quad \mathrm{D}_\varphi:=\operatorname{diag}(D_\varphi,D_\varphi,D_\varphi,D_\varphi).\] One can verify that these two operators commute with \(l_{\vec{F}}\): \[\mathrm{D}_t\circ l_{\vec{F}}=l_{\vec{F}}\circ \mathrm{D}_t, \quad \mathrm{D}_\varphi\circ l_{\vec{F}}=l_{\vec{F}}\circ \mathrm{D}_\varphi .\] Consequently, if \(\vec{\phi}\) is a solution to (22 ), then \(\mathrm{D}_t \vec{\phi}\) and \(\mathrm{D}_\varphi \vec{\phi}\) are also solutions. These operators are known as recursion operators since they allow us to obtain new generating functions of symmetries from known generating functions. In particular, repeated applications of \(\mathrm{D}_t\) and \(\mathrm{D}_\varphi\) to the generating functions \(\vec{\phi}_3\) (\(\vec{\phi}_1=\mathrm{D}_t\vec{\phi}_3\), \(\vec{\phi}_2=\mathrm{D}_\varphi\vec{\phi}_3\)) and \(\vec{\phi}_5\) produce generalized symmetries of arbitrarily high order. More precisely, if \(\vec{\phi}\) is a generating function of order \(q\), then \(\mathrm{D}_t^k\vec{\phi}\) and \(\mathrm{D}_\varphi^k\vec{\phi}\) are generating functions of order \(q+k\). This shows that the JCM differential system admits generalized symmetries of arbitrary order.

We also notice that the direct calculation of second-order generating functions does not yield additional independent symmetries beyond those obtained through these two recursion operators.

3.1 Point symmetries and invariant solutions↩︎

In Lie symmetry theory, point symmetries admitted by a system of differential equations can be used to construct classes of invariant solutions. These are solutions that remain unchanged under the action of a selected symmetry group. After substituting the corresponding invariant ansatz into the original system, one usually obtains a reduced system involving fewer independent variables or a simpler dependence on the original variables. This reduced system, often called a factor-system, is typically more accessible analytically and provides a systematic route for deriving exact solutions.

Using the relation between infinitesimal generators and generating functions, Eq. (24 ), applied to the first four generating functions in Eqs. (28 )–(31 ), we obtain the point symmetry generators \[X_1=\partial_t,\quad X_2=\partial_\varphi,\quad X_3=u^\alpha\partial_{u^\alpha},\quad X_4=\tilde{u}^\alpha\partial_{u^\alpha}.\] The generators \(X_1\) and \(X_2\) correspond to translations in \(t\) and \(\varphi\), respectively, while \(X_3\) generates a scaling transformation of the dependent variables \(u^\alpha\). The generator \(X_4\) represents the symmetry associated with the superposition of solutions of linear system of PDEs.

The Lie algebra of point symmetries with non-zero commutators \([X_i,X_4] = X_4\), \(i=1,2,3\) consists of the Abelian subalgebra \(L_3 = \langle X_1, X_2, X_3 \rangle\) and the ideal \(X_4\). The set of non-equivalent one-dimensional subalgebras of \(L_3\) is given by[22]: \[\label{algebra} \Theta_1=\left\langle X_2 \right\rangle,\; \Theta_2=\left\langle X_1+\alpha_0 X_2 \right\rangle,\; \Theta_3=\left\langle X_3+\alpha_0 X_1+\beta_0 X_2 \right\rangle,\tag{33}\] where \(\alpha_0\) and \(\beta_0\) are arbitrary constants.

3.1.1 Invariant solution for \(\Theta_1\)↩︎

To determine the invariant solutions associated with each subalgebra, we first find the invariants \(J\) of the corresponding basis operators. For the subalgebra \(\Theta_1=\langle X_2\rangle\), the invariance condition \(X_2J=0\) gives \(J_\alpha=u^\alpha,\) \(J_5=t\), \(J_6=r\). Therefore, the invariant solution has the form \[u^\alpha(\vec{r},t)=f^\alpha(r,t),\] that is, it has no explicit dependence on the angular variable \(\varphi\). Substitution of this ansatz into the dynamical system (17 ) gives the factor-system \[\begin{align} \partial_t f^1 &=& \varkappa f^2, \nonumber\\ \partial_t f^2 &=& -\delta f^3-\varkappa r\left(r f^1-2\partial_r f^4\right), \nonumber\\ \partial_t f^3 &=& \delta f^2, \nonumber\\ \partial_t f^4 &=& {2\varkappa\over r}\partial_r f^2. \label{eff1} \end{align}\tag{34}\] The equation for \(f^2\) can be decoupled by taking a second derivative with respect to \(t\). Using the remaining equations of the factor-system, one obtains \[\partial_t^2 f^2+\delta^2 f^2 = 4\varkappa^2 R\left(\partial_R^2 f^2-f^2\right).\] We now use separation of variables, \(f^2=T(t)Y(R)\). Substitution into the previous equation leads to \[\begin{align} T^{\prime\prime}+(\delta^2+4\varkappa^2\ell)T&=&0, \tag{35}\\ Y^{\prime\prime}-\left(1-\frac{\ell}{R}\right)Y&=&0, \tag{36} \end{align}\] where \(-4\varkappa^2\ell\), with \(\ell\neq 0\), has been chosen as the separation constant.

By introducing \[Y(R)=R\mathrm{e}^{-R}y(\tilde{x}),\qquad \tilde{x}=2R,\] Eq. (36 ) becomes Kummer’s equation, \[\tilde{x}y^{\prime\prime}+(2-\tilde{x})y^\prime +\left(\frac{\ell}{2}-1\right)y=0.\] Its general solution is [23] \[y(\tilde{x}) = c_1\,{}_1F_1(1-\ell/2;2;\tilde{x}) + c_2\,U(1-\ell/2;2;\tilde{x}),\] where \({}_1F_1\) is the confluent hypergeometric function of the first kind and \(U\) is the confluent hypergeometric function of the second kind. The asymptotic behavior of these functions as \(R\to\infty\) is given by (see Eqs. 13.2.23 and 13.2.6 of Ref. [24]) \[\begin{align} {}_1F_1(a;b;z) &\sim& \frac{\mathrm{e}^z z^{a-b}\Gamma(b)}{\Gamma(a)}, \qquad a\neq 0,-1,\ldots,\\ U(a,b,z) &\sim& z^{-a}. \end{align}\] Thus, one must either set \(c_1=0\) or impose \[1-\frac{\ell}{2}=-n,\qquad n=0,1,\ldots .\] In the first case, the function \(U(1-\ell/2;2;2R)\) is one-valued [23] only if either (i) \(1-\ell/2=-n\), in which case it is a polynomial in \(R\) and is proportional to \({}_1F_1\), or (ii) \(1-\ell/2=1\), corresponding to \(\ell=0\). Therefore, the only one-valued solution that is compatible with finiteness and integrability of the chord function is proportional to \[R\mathrm{e}^{-R}\,{}_1F_1(-n;2;2R), \qquad \ell=2(n+1).\]

Using the relation between confluent hypergeometric functions and associated Laguerre polynomials, \[\label{hyptoLag} L^{(m)}_n(x) = \frac{(m+n)!}{m!n!}\,{}_1F_1(-n;m+1;x),\tag{37}\] the radial solution can be written, up to an overall constant, as \[Y_n(r):=r^2\mathrm{e}^{-r^2/4}L^{(1)}_n(r^2/2).\] Furthermore, using the identity \[L^{(1)}_n(x) = {n+1\over x}\left[L_n(x)-L_{n+1}(x)\right],\] we can equivalently write \[Y_n(r) = 2(n+1)\mathrm{e}^{-r^2/4} \left[ L_n(r^2/2)-L_{n+1}(r^2/2) \right].\]

From Eq. (35 ), together with the quantization condition \(\ell=2(n+1)\), the time-dependent part becomes \[T^+_n(t) := \tilde{c}_1\mathrm{e}^{\mathrm{i}\lambda_n t} + \tilde{c}_2\mathrm{e}^{-\mathrm{i}\lambda_n t},\] where \(\lambda_n\) is the generalized Rabi frequency defined in Eq. (80 ). The time-dependent functions \(T^+_n(t)\) and \(T^-_n(t)\), defined below, reproduce the same oscillatory structure as the dressed-state solution in Eqs. (83 )–(84 ), encoded in the amplitudes \(a_n(t)\) and \(b_n(t)\) in Eq. (85 ), for appropriate initial conditions (see Appendix 7).

Thus, the invariant solutions \(\vec{u}_n\) obtained from \(\Theta_1\) and associated with the frequency \(\lambda_n\) are \[\begin{align} u^1_n(\vec{r},t) &=& -{\mathrm{i}g\over \lambda_n\sqrt{2}}\,Y_n(r)T^-_n(t)+U^1(r), \tag{38}\\ u^2_n(\vec{r},t) &=& Y_n(r)T^+_n(t),\\ u^3_n(\vec{r},t) &=& -{\mathrm{i}\delta\over\lambda_n}\,Y_n(r)T^-_n(t)+U^3(r),\\ u^4_n(\vec{r},t) &=& -{\mathrm{i}g\sqrt{2}\over\lambda_n}\, {Y_n^\prime(r)\over r}T^-_n(t)+U^4(r), \tag{39} \end{align}\] where \[\begin{align} T^-_n(t) &:=& \tilde{c}_1\mathrm{e}^{\mathrm{i}\lambda_n t} - \tilde{c}_2\mathrm{e}^{-\mathrm{i}\lambda_n t}, \tag{40}\\ U^1(r) &:=& {2\over r}U^{4\,\prime}(r) - {\delta\sqrt{2}\over gr^2}U^3(r), \tag{41} \end{align}\] and \(U^3(r)\) and \(U^4(r)\) are arbitrary functions.

Using the inverse transformation between the variables \(u^\alpha\) and the characteristic-function components, we have \[\begin{align} w_{ee}&=\!\!&\!\! {u^1_n(\vec{r},t)+u^4_n(\vec{r},t)\over 2}, \!\!\!\quad\!\!\! w_{gg} \!=\! {u^1_n(\vec{r},t)-u^4_n(\vec{r},t)\over 2}, \\ w_{eg} &=& (k+\mathrm{i}s){u^2_n(\vec{r},t)\over 2r^2} + (s-\mathrm{i}k){u^3_n(\vec{r},t)\over 2r^2}, \tag{42}\\ w_{ge} &=& (-k+\mathrm{i}s){u^2_n(\vec{r},t)\over 2r^2} - (s+\mathrm{i}k){u^3_n(\vec{r},t)\over 2r^2}. \tag{43} \end{align}\] One can then verify that the initial conditions (12 ) and the conditions a), b), and c) of Appendix 6 for the characteristic function \(w_f(\vec{r},t)\) are satisfied if \(a_1\) and \(b_1\) are real constants, \(U^1(r),U^3(r),U^4(r)\in\mathbb{R}\), \(U^1(r)\) is bounded, \(U^1(0)=1\), and \[\begin{align} T^+_n(t) &=& a_1\cos\lambda_n t-b_1\sin\lambda_n t,\\ T^-_n(t) &=& \mathrm{i}\left(a_1\sin\lambda_n t+b_1\cos\lambda_n t\right). \end{align}\]

The condition of initially decoupled atom–field states, Eq. (13 ), together with the boundary conditions and the conditions a), b), and c) of Appendix 6, is fulfilled by choosing \[\begin{align} u^1_n(\vec{r},t) &=& {cg\over \lambda_n\sqrt{2}}\,Y_n(r)\cos\lambda_n t +U^1_n(r), \tag{44}\\ u^2_n(\vec{r},t) &=& -cY_n(r)\sin\lambda_n t,\\ u^3_n(\vec{r},t) &=& {c\delta\over\lambda_n}\,Y_n(r)(\cos\lambda_n t-1),\\ u^4_n(\vec{r},t) &=& {cg\sqrt{2}\over\lambda_n} {Y_n^\prime(r)\over r}(\cos\lambda_n t-1), \tag{45} \end{align}\] where \[\begin{align} T^+_n(t)&=&-c\sin\lambda_n t, \qquad T^-_n(t)=\mathrm{i}c\cos\lambda_n t,\\ U^3_n(r)&=& -{c\delta\over\lambda_n}Y_n(r), \qquad U^4_n(r)= -{cg\sqrt{2}\over\lambda_n}{Y'_n(r)\over r}, \end{align}\] with \[c:={\varkappa\over (n+1)\lambda_n}.\] In this case, the reduced field characteristic function \(w_f(\vec{r},t)=u^1_n(\vec{r},t)\) takes the explicit form \[w^{(n)}_f(\vec{r},t) = {cg(\cos\lambda_n t-1)\over\lambda_n\sqrt{2}}Y_n(r) + {c\lambda_n\sqrt{2}\over g}{Y_n(r)\over r^2}.\] Moreover, \(\rho_{eg}(t)=\rho_{ge}(t)\equiv 0\), while \(w_z(\vec{r},t)=u^4_n(\vec{r},t)\). To simplify the expression for \(U^1_n(r)\), we use the relation \[{Y^\prime_n(r)\over r} = {Y_n(r)\over 2} + 2(n+1)\mathrm{e}^{-r^2/4}L_{n+1}(r^2/2).\]

3.1.2 Invariant solution for \(\Theta_2\)↩︎

We now consider the invariant solutions associated with the subalgebra \(\Theta_2\) in Eq. (33 ). Since this subalgebra depends on the parameter \(\alpha_0\), different choices of \(\alpha_0\) lead to different classes of invariant solutions. The invariants are obtained from \[\left(\partial_t+\alpha_0\partial_{\varphi}\right)J=0.\] Thus, \[J_\alpha=u^\alpha,\qquad J_5=\varphi-\alpha_0 t=:\theta,\qquad J_6=r,\] and the corresponding invariant solution has the form \[u^\alpha(\vec{r},t)=f^\alpha(\theta,r).\]

For the particular value \(\alpha_0=1\), one has \(\hat{L} f^\alpha=0\), and the system (17 ) leads to the stationary part of the solution (38 )–(39 ), obtained by setting \(T_n^\pm\equiv 0\). We therefore focus on the nontrivial case \[\alpha_0=1+\frac{\varkappa}{a},\qquad a\in\mathbb{R}\setminus\{0\}.\] In this case, the invariant ansatz yields the factor-system \[\begin{align} \partial_\theta f^1 &=& -a f^2, \tag{46}\\ \partial_\theta f^2 &=& {a\delta\over\varkappa}f^3+4aR\left(f^1-\partial_R f^4\right), \tag{47}\\ \partial_\theta f^3 &=& -{a\delta\over\varkappa}f^2-2a\partial_\theta f^4, \tag{48}\\ \partial_\theta f^4 &=& -a\partial_R f^2-{a\over 2R}\partial_\theta f^3, \tag{49} \end{align}\] where, as before, \(R=r^2/4\).

After eliminating \(f^1\), \(f^3\), and \(f^4\) from Eqs. (46 )–(49 ), one obtains a closed equation for \(f^2\): \[\label{ECV22} \partial_\theta^2 f^2 = B(R)\partial_R^2 f^2 + C(R)\partial_R f^2 + D(R)f^2,\tag{50}\] where \[\begin{align} B(R) &=& {4(aR)^2\over R-a^2}, \nonumber\\ C(R) &=& -{4a^4R\over (R-a^2)^2}, \nonumber\\ D(R) &=& -{4a^2R\over (R-a^2)^2} \left[ (R-a^2)^2 + {\delta^2\over 4\varkappa^2}(R-a^2) - {a\delta\over 2\varkappa} \right].\nonumber\\ \end{align}\]

We now use separation of variables, \(f^2=T(\theta)Y(R)\). Substitution into Eq. (50 ) gives \[\begin{align} {T^{\prime\prime}\over T} &=& B(R){Y^{\prime\prime}\over Y} + C(R){Y^\prime\over Y} + D(R) = \ell, \label{ECV22b} \end{align}\tag{51}\] where \(\ell\) is the separation constant. Hence, the radial function satisfies \[\label{Ydeq} B(R)Y^{\prime\prime}+C(R)Y^\prime+\left[D(R)-\ell\right]Y=0.\tag{52}\] With the change of variables \[Y(R)=\mathrm{e}^{-R}R^{\sqrt{-\ell}/2}h(z), \qquad z={R\over a^2},\] Eq. (52 ) is transformed into the confluent Heun equation \[\label{C95HEUN} h^{\prime\prime} + \left( \alpha_h + {\beta_h+1\over z} + {\gamma_h+1\over z-1} \right)h^\prime + \left( {\mu_h\over 2z} + {\nu_h\over 2(z-1)} \right)h =0,\tag{53}\] where \[\begin{align} \mu_h &=& \alpha_h\beta_h-\beta_h\gamma_h+\alpha_h-\beta_h-\gamma_h-2\eta_h, \nonumber \\ \nu_h &=& \alpha_h\gamma_h+\beta_h\gamma_h+\alpha_h+\beta_h+\gamma_h +2\delta_h+2\eta_h, \nonumber \end{align}\] and \[\begin{align} \alpha_h&=&-2a^2, \qquad \beta_h=\sqrt{-\ell}, \qquad \gamma_h=-2, \nonumber\\ \delta_h &=& {(4a^4-\ell)\varkappa^2-a^2\delta^2\over 4\varkappa^2}, \quad \eta_h = 1-\delta_h+{a\delta\over 2\varkappa}. \label{HH} \end{align}\tag{54}\]

Equation (53 ) has a unique particular solution, which is regular around the regular singular point \(z = 0\), given by \[\label{H95SERIES} h(z) = \mathrm{HeunC}\left(\alpha_h, \beta_h, \gamma_h, \delta_h, \eta_h; z\right) = \sum\limits_{k=0}^\infty v_k z^k,\tag{55}\] This series is convergent in the disk \(|z|<1\), and the coefficients \(v_k\) are determined by a three-term recurrence relation (see Ref. [25]). To obtain a convergent solution on the whole interval \(z\in[0,\infty)\), one can consider a polynomial solution, the necessary condition of which is \[-{\delta_h\over\alpha_h} = N_h+1+{\beta_h+\gamma_h\over 2}, \qquad N_h\in\mathbb{N}.\] Using the parameter values in Eq. (54 ), this condition gives the quantization of the separation constant, \[\ell_n = -2{a^2\over g^2} \left(\lambda_n\pm ag\sqrt{2}\right)^2,\] where \(\lambda_n\) is defined in Eq. (80 ) and \(N_h=n+1\). Moreover, for \(\ell=\ell_n\), the coefficients in the series (55 ) vanish from \(v_{n+2}\) onward. Therefore, \[\begin{align} h_{n+1}(R) &:=& \mathrm{HeunC}\left( -2a^2,\sqrt{-\ell_n},-2, a^4-{\ell_n\over 4}-{a^2\delta^2\over 2g^2}, \right. \nonumber\\ && \left. 1-a^4+{\ell_n\over 4}+{a^2\delta^2\over 2g^2} +{a\delta\over g\sqrt{2}}, {R\over a^2} \right) \end{align}\] is a polynomial of degree \(n+1\) in \(R\), with \(h_{n+1}(0)=1\). Consequently, the radial part of \(f^2(\theta,r)\) is \[\label{Rpart2invsol} Y_n(r)=\mathrm{e}^{-r^2/4}r^{\sqrt{-\ell_n}}h_{n+1}(r^2/4).\tag{56}\]

The temporal part is obtained directly from Eq. (51 ) as \[\label{Tpart2invsol} T_n^\pm(\theta) = \tilde{c}_1\mathrm{e}^{\mathrm{i}\sqrt{|\ell_n|}\theta} \pm \tilde{c}_2\mathrm{e}^{-\mathrm{i}\sqrt{|\ell_n|}\theta}.\tag{57}\] The corresponding invariant solution is \[\begin{align} \tag{58} u^1_n(\vec{r}, t ) &=&\frac{\mathrm{i}a}{\sqrt{|\ell_n|}}\, Y_n(r) T^{-}_n(\theta) + U^1(r),\\ u^2_n(\vec{r}, t ) &=& Y_n(r) T^+_n(\theta),\\ u^3_n(\vec{r}, t ) &=& \frac{\mathrm{i}a}{\sqrt{|\ell_n|}} \frac{r}{r^2 - 4a^2}\left[{\delta\sqrt{2}}r Y_n(r)/{g} \right. \nonumber \\ &-& \left. 4aY_n^\prime(r) \right] T^{-}_n(\theta) + U^3(r),\\\tag{59} u^4_n(\vec{r}, t ) &=& \frac{2\mathrm{i}a}{\sqrt{|\ell_n|}} \frac{1}{r^2 - 4a^2}\left[-{a\delta\sqrt{2}} Y_n(r)/{g} \right. \nonumber \\ &+& \left. rY_n^\prime(r) \right] T^{-}_n(\theta) + U^4(r), \end{align}\] where \(U^1(r)\), \(U^3(r)\), and \(U^4(r)\) are the same arbitrary functions as in Eq. (41 ).

It is useful to rewrite \(T_n^\pm(\theta)\) in terms of the original variables. Since \(\theta=\varphi-\alpha_0 t\) and \(s \pm \mathrm{i}k = r\mathrm{e}^{\pm \mathrm{i}\varphi}\), we have \[\begin{align} T_n^\pm(\theta) &=& \left[ \tilde{c}_1 (s+\mathrm{i}k)^{\sqrt{|\ell_n|}} \mathrm{e}^{-\mathrm{i}\alpha_0\sqrt{|\ell_n|}t} \right. \nonumber\\ && \left. \pm \tilde{c}_2 (s-\mathrm{i}k)^{\sqrt{|\ell_n|}} \mathrm{e}^{\mathrm{i}\alpha_0\sqrt{|\ell_n|}t} \right] r^{-\sqrt{|\ell_n|}} . \end{align}\] The normalization condition for the reduced field characteristic function, \(w_f(0,t)=u^1_n(0,t)=1\), is fulfilled if \(U^1(0)=1\).

The Hermiticity condition \[w_f(-\vec{r},t)=w_f^\ast(\vec{r},t)\] leads to two possible cases. In the even case, with \(m\in\mathbb{N}\), \[\tilde{c}_2=\tilde{c}_1^\ast, \; \ell_n=-4m^2, \; a=a_{2m} = {\sqrt{\lambda_n^2+8g^2m}-\lambda_n\over 2\sqrt{2}g}.\] In the odd case, \[\begin{align} \tilde{c}_2&=&-\tilde{c}_1^\ast, \quad \ell_n=-(2m+1)^2, \nonumber\\ a&=&a_{2m+1} = {\sqrt{\lambda_n^2+4g^2(2m+1)}-\lambda_n\over 2\sqrt{2}g}. \end{align}\] In both cases, \(U^1(r)\) is taken to be real. Thus, \[\begin{align} w_f^{\rm even}(\vec{r},t) &:=& \mathrm{i}{a_{2m}\over 2m} \mathrm{e}^{-r^2/4} h_{n+1}(r^2/4) \nonumber\\ &&\times \left[ \tilde{c}_1 (s+\mathrm{i}k)^{2m} \mathrm{e}^{-\mathrm{i}2m\alpha_0 t} - \tilde{c}_1^\ast (s-\mathrm{i}k)^{2m} \mathrm{e}^{\mathrm{i}2m\alpha_0 t} \right] \nonumber\\ && +U^1(r) \end{align}\] is real-valued, whereas \[\begin{align} w_f^{\rm odd}(\vec{r},t) &:=& \mathrm{i}{a_{2m+1}\over 2m+1} \mathrm{e}^{-r^2/4} h_{n+1}(r^2/4) \nonumber\\ &&\times \left[ \tilde{c}_1 (s+\mathrm{i}k)^{2m+1} \mathrm{e}^{-\mathrm{i}(2m+1)\alpha_0 t} \right. \nonumber\\ && \left. + \tilde{c}_1^\ast (s-\mathrm{i}k)^{2m+1} \mathrm{e}^{\mathrm{i}(2m+1)\alpha_0 t} \right] +U^1(r)\nonumber\\ \end{align}\] is, in general, complex-valued.

To satisfy the conditions (12 ), we take \(U^4(r)\in\mathbb{R}\) and restrict attention to the even case. In the odd case, the requirement that \([\mathrm{w}^0(\vec{r})]_{11}\) and \([\mathrm{w}^0(\vec{r})]_{22}\) be real forces \(\tilde{c}_1=0\), which reduces \(w_f^{\rm odd}\) to a static characteristic function. Furthermore, the Hermiticity condition at the origin, \[[\mathrm{w}^0(0)]_{12}=[\mathrm{w}^0(0)]_{21}^\ast,\] implies, in the limit \(s,k\to0\), that \[[\mathrm{w}^0(0)]_{12}=[\mathrm{w}^0(0)]_{21}^\ast=0.\]

Finally, let us consider the more restrictive condition of initially decoupled atom–field states, Eq. (13 ). In the limit \(\vec{r}\to0\), one obtains \(\rho_{eg}=\rho_{ge}=0\). Hence, for a factorized initial state, \([\mathrm{w}^0(\vec{r})]_{12}\) and \([\mathrm{w}^0(\vec{r})]_{21}\) must vanish identically. Using Eqs. (42 ) and (43 ), this requires \[u^2_n(\vec{r},0)\equiv 0, \qquad u^3_n(\vec{r},0)\equiv 0,\] which is possible only if \(\tilde{c}_1=0\). Therefore, in the factorized case the solution becomes static. Consequently, the nontrivial even invariant solution obtained from \(\Theta_2\) describes generic coupled atom–field initial states, rather than initially decoupled ones.

3.1.3 Invariant solution for \(\Theta_3\)↩︎

We finally consider the subalgebra \(\Theta_3\). For a nontrivial invariant solution associated with this subalgebra, it is necessary that \[\alpha_0^2+\beta_0^2\neq 0.\] We analyze the possible cases separately and show that none of them satisfies the boundary condition required for the characteristic function.

First, let \(\beta_0=0\). In this case, the invariant solution has the form \[u^\alpha(\vec{r},t)=\mathrm{e}^{t/\alpha_0}f^\alpha(\varphi,r).\] The normalization condition for the reduced field characteristic function, namely condition a) in Appendix 6, \(w_f(0,t)=1\), becomes \[f^1(\varphi,r)\big|_{\vec{r}=0}=\mathrm{e}^{-t/\alpha_0}.\] This condition cannot be satisfied for arbitrary \(t\), since the left-hand side does not depend on time.

Second, let \(\alpha_0=0\). Then the invariant solution takes the form \[u^\alpha(\vec{r},t)=\mathrm{e}^{\varphi/\beta_0}f^\alpha(t,r).\] Using \(s \pm \mathrm{i}k = r\mathrm{e}^{ \pm \mathrm{i}\varphi}\), the first component may be written as \[u^1 = \left( {s+\mathrm{i}k\over s-\mathrm{i}k} \right)^{1/(2\mathrm{i}\beta_0)} f^1(t,k^2+s^2).\] However, as \(\vec{r}\to0\), the angular factor has no unique limit. One might try to compensate this behavior by introducing a radial power through \(f^1(t,r)\), for instance using \[(s^2+k^2)^a=(s+\mathrm{i}k)^a(s-\mathrm{i}k)^a,\] with some exponent \(a\). Nevertheless, this cannot remove the angular dependence in a way that yields the required limit for the normalization condition; therefore, this case does not lead to an admissible characteristic function.

Finally, consider the case \(\alpha_0\beta_0\neq0\). The invariant solution has the form \[u^\alpha(\vec{r},t)=\mathrm{e}^{t/\alpha_0}f^\alpha(z,r),\] where \[z= {k-s\tan\left(\frac{\beta_0}{\alpha_0}t\right) \over s+k\tan\left(\frac{\beta_0}{\alpha_0}t\right)}.\] The variable \(z\) has no well-defined limit also at \(\vec{r}\to0\). Thus, either \(f^1\) must be independent of \(z\), so that \(f^1=f^1(r)\), or it must depend on a regularized combination such as \(zr^a\), with \(a>0\), so that \(zr^a\to0\) as \(r\to0\). In both cases, however, the value of \(f^1\) at the origin cannot compensate the prefactor \(\mathrm{e}^{t/\alpha_0}\) for arbitrary \(t\). Consequently, the normalization condition cannot again be satisfied.

We conclude that the invariant solutions associated with \(\Theta_3\) are not compatible with the boundary condition imposed by the characteristic-function description. Therefore, \(\Theta_3\) does not yield physically admissible invariant solutions for the JCM in the present framework.

4 Conservation laws↩︎

A conservation law associated with a system of differential equations (18 ), is represented by a conserved current \[\label{Acons} \vec{A}=(A^1,\ldots,A^n),\tag{60}\] whose divergence vanishes on the solution manifold of the system. In terms of total derivatives, this condition is written as \[\left.\nabla\cdot\vec{A}\right|_{\mathcal{F}=0} = \left. \left( D_1A^1+\cdots+D_nA^n \right) \right|_{\mathcal{F}=0} =0,\] where \(D_j\) is the total derivative with respect to \(x^j\), as defined in Eq. (21 ). In general, the components \(A^i\) may depend on the independent variables, the dependent variables, and a finite number of their derivatives, \[A^i=A^i(\vec{x},\vec{u},\vec{u}^{(1)},\ldots,\vec{u}^{(s)}).\]

Finding conserved currents directly is, in general, a nontrivial problem. However, the methodology introduced by Vinogradov  [16]–[18] provides a systematic way to associate local conservation laws with their generating functions \[\label{genfuncons} \vec{\psi}=(\psi^1,\ldots,\psi^m) .\tag{61}\] For a system of PDEs called \(l\)-normal system (see Chapter 5 of Ref. [18]), the generating functions are uniquely defined by conservation laws. Conservation laws that differ by a null divergence correspond to the same generating function. The generating function \(\vec{\psi}\) is the solution of the following equation (see also Ref. [26]) \[\label{sist95func95gen} \left. l^\star_{\vec{F}}\vec{\psi} \right|_{\mathcal{F}=0} =0.\tag{62}\] The formal adjoint \(^{\star}\) of a scalar differential operator is defined by \[\left( \sum_{\varsigma}\alpha_\varsigma D^{(\varsigma)} \right)^\star = \sum_{\varsigma} (-1)^{|\varsigma|} D^{(\varsigma)}\circ \alpha_\varsigma .\] For the matrix differential operator \(l_{\vec{F}}\), the adjoint is obtained by taking the formal adjoint of each entry and transposing the matrix: \[l^\star_{\vec{F}} = \left(\left[l_{\vec{F}}\right]^\star_{\beta\alpha}\right), \quad {\alpha,\beta=1,\dots,m}.\]

A solution \(\vec{\psi}\) of Eq. (62 ), sometimes called an adjoint symmetry or co-symmetry, corresponds to a conservation law if and only if there exists a differential operator \(\mathcal{T}\) such that \[\label{criterio} l_{\vec{\psi}} + \left.\square\right|_{\mathcal{F}}^{\star} = \left.\mathcal{T}\right|_{\mathcal{F}}\circ l_{\vec{F}}, \qquad \mathcal{T}^\star=\mathcal{T},\tag{63}\] where \(\square\vec{F}=l^\star_{\vec{F}}\vec{\psi}\), and \(l_{\vec{\psi}}\) denotes the linearization of \(\vec{\psi}\). We recall that a system \(\vec{F}=0\) is an Euler–Lagrange system if and only if \[l_{\vec{F}}=l^\star_{\vec{F}}.\] If \(\vec{\psi}\) satisfies the criterion (63 ), then the associated conserved current can be obtained from the relation \[\label{CLaw} \nabla\cdot\vec{A} = S_\alpha F^\alpha,\tag{64}\] where summation over \(\alpha\) is understood. The operators \(S_\alpha\) are defined by the condition \[\label{S} \psi^\alpha = \left. S_\alpha^\star 1 \right|_{\mathcal{F}=0}.\tag{65}\]

We now apply this construction to the JCM differential system \(\vec{F}=(F^1,F^2,F^3,F^4) =0\) given in Eq. (17 ). The system of equations under consideration is a determined regular system of equations of evolutionary type (i.e., the derivatives with respect to time are expressed explicitly through the derivatives with respect to the remaining variables), therefore, it is \(l\)-normal. Since the linearization operator (27 ) is written in the variable \(R=r^2/4\), we first work with the current \(\vec{A}=(A^t,A^R,A^\varphi)\). Thus, \[\label{div} D_tA^t+D_RA^R+D_\varphi A^\varphi = S_\alpha F^\alpha .\tag{66}\]

The formal adjoint of the linearization operator (27 ) is \[l^\star_{\vec{F}} = \begin{pmatrix} -D_t-D_\varphi & 4\varkappa R & 0 & 0\\ -\varkappa& -D_t-D_\varphi & -\delta & \varkappa D_R\\ 0 & \delta & -D_t-D_\varphi & \frac{\varkappa}{2R}D_\varphi\\ 0 & 4\varkappa(1+RD_R) & 2\varkappa D_\varphi & -D_t-D_\varphi \end{pmatrix}.\] Since \[l_{\vec{F}}\neq l^\star_{\vec{F}},\] the system (17 ) is not an Euler–Lagrange system, and therefore the conservation laws cannot be obtained directly from Noether’s theorem. Nevertheless, the operators \(l_{\vec{F}}\) and \(l^\star_{\vec{F}}\) are related by \[l^\star_{\vec{F}}\circ E_R=l_{\vec{F}}, \qquad E_R= \operatorname{diag} \left( 1,(4R)^{-1},(4R)^{-1},1 \right).\] This relation provides a direct connection between symmetry generating functions and co-symmetries. In particular, if \(\vec{\phi}\) is a generating function of a symmetry satisfying Eq. (22 ), then \[\vec{\psi}=E_R\vec{\phi}\] satisfies the adjoint equation (62 ). Equivalently, \[\vec{\phi}=E_R^{-1}\vec{\psi}.\]

4.1 Zero order generating functions↩︎

The zero-order symmetry generating functions are \[\vec{\phi}_3=\vec{u}(\vec{r},t), \qquad \vec{\phi}_4=\vec{\tilde{u}}(\vec{r},t),\] where \(\vec{\tilde{u}}\) is an arbitrary solution of the linear system (17 ). From the relation between symmetries and co-symmetries discussed above, the corresponding co-symmetries are \[\vec{\psi}_3=E_R\vec{u}, \qquad \vec{\psi}_4=E_R\vec{\tilde{u}}.\] In particular, \[\vec{\psi}_3 = \left( u^1,\frac{u^2}{4R},\frac{u^3}{4R},u^4 \right) .\] For this co-symmetry, the operator \(\square\) is given by \(\square=-E_R\), and therefore it is self-adjoint, \(\square^\star=\square\). Moreover, \(l_{\vec{\psi}_3}=E_R\), so that \[l_{\vec{\psi}_3}+\square^\star=\mathcal{O},\] where \(\mathcal{O}\) denotes the zero operator. Hence the criterion (63 ) is satisfied with \(\mathcal{T}=\mathcal{O}\), and \(\vec{\psi}_3\) is the generating function of a conservation law.

From Eq. (65 ), we have \[\vec{S}= \left( u^1,\frac{u^2}{4R},\frac{u^3}{4R},u^4 \right).\] Integrating Eq. (64 ), we obtain the conserved current \[\begin{align} A^t &=& \frac{(u^1)^2+(u^4)^2}{2} + \frac{(u^2)^2+(u^3)^2}{8R}, \\ A^R &=& -\varkappa u^2u^4, \\ A^\varphi &=& A^t-\frac{\varkappa}{2R}u^3u^4. \end{align}\] Transforming back to the Cartesian Fourier variables \((s,k)\) and using the inverse transformations (16 ) and (8 ), the current can be written as \[\begin{align} A^t &=& w_{ee}^2+w_{gg}^2-2w_{eg}w_{ge}, \\ A^s &=& 2\mathrm{i}\varkappa\, w_z w_x-kA^t, \\ A^k &=& 2\mathrm{i}\varkappa\, w_z w_y+sA^t . \end{align}\]

Although the full current is defined in the atom-field Fourier phase space, its evaluation at the origin gives a useful interpretation in terms of the reduced atomic state. At \(\vec{r}=0\), one obtains \[\begin{align} A^t\big|_{\vec{r}=0} &=& P_a(t)-f_a^2(t), \\ A^s\big|_{\vec{r}=0} &=& 2\mathrm{i}\varkappa\langle\sigma_x\rangle\langle\sigma_z\rangle, \\ A^k\big|_{\vec{r}=0} &=& 2\mathrm{i}\varkappa\langle\sigma_y\rangle\langle\sigma_z\rangle, \end{align}\] where \[P_a(t)=\varrho_{ee}^2(t)+\varrho_{gg}^2(t)+2|\varrho_{eg}(t)|^2\] is the purity of the reduced atomic density matrix, and \[f_a(t)=2|\varrho_{eg}(t)|\] is the normalized atomic coherence function. Thus, the density component of the conserved current is not the atomic purity alone, but the combination \[P_a(t)-f_a^2(t) = \varrho_{ee}^2(t)+\varrho_{gg}^2(t)-2|\varrho_{eg}(t)|^2 .\] This quantity measures a balance between the atomic population contribution to the reduced-state purity and the contribution associated with atomic coherence.

The conservation equation \[D_tA^t+D_sA^s+D_kA^k=0\] then gives, at the origin, \[\begin{align} {\mathrm{d}\over\mathrm{d}t} \left(P_a(t)-f_a^2(t)\right) &=& -2\mathrm{i}\varkappa \left. \left[ \partial_s(w_z w_x)+\partial_k(w_z w_y) \right] \right|_{\vec{r}=0} \nonumber\\ && - \left. \left(\partial_\varphi A^t\right) \right|_{\vec{r}=0}. \label{Purity95flux1} \end{align}\tag{67}\] This equation shows that, in the absence of atom–field coupling (\(\varkappa=0\)), and for regular solutions for which the angular contribution at the origin vanishes, the quantity \(P_a-f_a^2\) is conserved. When the coupling is present, its rate of change is controlled by the coupling strength and by derivatives of products of characteristic-function components. These terms encode atom–field correlations generated by the interaction and therefore describe the exchange of information between the atomic reduced state and the field degrees of freedom.

Using the relations in Eq. (76 ), Eq. (67 ) can be written as \[\begin{align} {\mathrm{d}\over\mathrm{d}t} \left(P_a(t)-f_a^2(t)\right) &=& 2\varkappa \left( \langle\sigma_z\rangle\langle \hat{p}_x^1\rangle + \langle\sigma_z\rangle\langle \hat{x}_y^1\rangle \right. \nonumber\\ && \left. + \langle\sigma_x\rangle\langle \hat{p}_z^1\rangle + \langle\sigma_y\rangle\langle \hat{x}_z^1\rangle \right) - \left. \left(\partial_\varphi A^t\right) \right|_{\vec{r}=0}.\nonumber\\\label{Purity95flux2} \end{align}\tag{68}\] In this form, the balance equation makes explicit that the variation of \(P_a-f_a^2\) is driven by atom-field coupling and by field moments weighted by atomic expectation values. Since, for a closed bipartite pure state, changes in the local purity are directly related to the generation of atom-field entanglement, Eqs. (67 ) and (68 ) provide a conservation-law expression for quantities whose evolution is linked to atom-field correlation and entanglement dynamics in the JCM.

For the invariant solutions derived in Secs. 3.1.1 and 3.1.2, the angular contribution satisfies \[\left. \left(\partial_\varphi A^t\right) \right|_{\vec{r}=0}=0.\] In particular, for the decoupled-state solutions (44 )–(45 ), one obtains \[\partial_t \left[ (u^1_n)^2+(u^4_n)^2+\left({u^2_n\over r}\right)^2 +\left({u^3_n\over r}\right)^2 \right] = 4\varkappa\,{\partial_r(u^2_nu^4_n)\over r}, \label{CL950}\tag{69}\] and \(f_a(t)\equiv0\). Therefore, \[2P_a(t) = \left[(u^1_n)^2+(u^4_n)^2\right]_{\vec{r}=0} = \frac{4g^4}{\lambda_n^4} \left(\cos\lambda_n t-1\right)^2+1.\] This illustrates explicitly how the local atomic purity oscillates under the atom–field interaction.

The second zero-order co-symmetry, \[\vec{\psi}_4 = \left( \tilde{u}^1, \frac{\tilde{u}^2}{4R}, \frac{\tilde{u}^3}{4R}, \tilde{u}^4 \right) ,\] also generates conservation laws. In this case \(l_{\vec{\psi}_4}=\mathcal{O}=\square^\star|_{\mathcal{F}}\), and the associated current is \[\begin{align} A^t &=& \tilde{u}^1u^1 + {\tilde{u}^2\over 4R}u^2 + {\tilde{u}^3\over 4R}u^3 + \tilde{u}^4u^4, \nonumber\\ A^R &=& -\varkappa\left(\tilde{u}^2u^4+\tilde{u}^4u^2\right), \\ A^\varphi &=& A^t - {\varkappa\over 2R} \left(\tilde{u}^3u^4+\tilde{u}^4u^3\right). \nonumber \end{align}\] This type of conservation law is common to linear systems of differential equations, where the superposition principle naturally produces conservation laws involving arbitrary solutions of the adjoint system (see Example 3.3, Chapter 5 of Ref. [18]).

A particularly important consequence is that the usual constant of motion of the JCM can be recovered from this construction. Consider Eq. (64 ) with \[S_1 = -D_R-RD_{RR}-{1\over 4R}D_{\varphi\varphi}, \qquad S_2=S_3=0, \qquad S_4=1.\] Using Eq. (65 ), this choice corresponds to the generating function \[\vec{\psi}_4=(0,0,0,1) .\] Therefore, there exists a current \(\vec{A}\) such that \[S_1F^1+F^4=\operatorname{div}\vec{A}.\] In the Cartesian variables \((k,s)\), the operator \(S_1\) is the negative Laplacian in Fourier phase space, \[S_1=-\Delta_{ks}=-(\partial_{ss}+\partial_{kk}).\] Thus, \[\label{OO} S_1F^1+F^4 = \partial_tu^4-\Delta_{ks}(\partial_tu^1) + \varkappa R\partial_{RR}u^2 + O(\varphi) = 0,\tag{70}\] where \(O(\varphi)\) denotes terms containing derivatives with respect to \(\varphi\). For the invariant solutions of Sec. 3.1.1, the \(\varphi\)-dependent terms vanish. Moreover, from Eq. (36 ), \(\partial_{RR}u^2=\left(1-{\ell\over R}\right)u^2\), \(u^2|_{R=0} = 0\), and the radial behavior of the admissible solutions implies \(R\partial_{RR}u^2\to0\), as \(R\to0\). Therefore, evaluating Eq. (70 ) at \(\vec{r}=0\) gives \[\left. \partial_t \left( w_z-\partial_{ss}w_f-\partial_{kk}w_f \right) \right|_{\vec{r}=0} =0.\] Using the characteristic-function identities \[\label{xp} \langle \hat{x}^2(t)\rangle = -\left.\partial_{kk}w_f(\vec{r},t)\right|_{\vec{r}=0}, \; \langle \hat{p}^2(t)\rangle = -\left.\partial_{ss}w_f(\vec{r},t)\right|_{\vec{r}=0},\tag{71}\] together with \[\langle\sigma_z(t)\rangle=w_z(\vec{r},t)\big|_{\vec{r}=0},\] we recover the conserved quantity \[\langle\hat{N}\rangle = \langle\sigma_z(t)\rangle + \langle\hat{x}^2(t)\rangle + \langle\hat{p}^2(t)\rangle = \mathrm{const.} \label{CONSTANT}\tag{72}\] which corresponds, up to the normalization conventions used for the quadrature variables, to the expectation value of the total number of excitations in the JCM.

4.2 Generating functions of non-zero order↩︎

We now consider conservation laws generated by co-symmetries depending on first and higher derivatives of the dependent variables. Starting from the symmetry generating functions \(\vec{\phi}_1,\;\vec{\phi}_2,\; \vec{\phi}_5\), given in Eqs. (28 ), (29 ), and (32 ), we obtain the corresponding adjoint symmetries through the relation \[\vec{\psi}_{1,2,5}=E_R\vec{\phi}_{1,2,5}.\] In particular, \[\begin{align} \nonumber \vec{\psi}_5 \!\!\!\!&=&\!\!\!\! \left( {\delta\over\varkappa}u^4 \!+\! {\partial_\varphi u^2\over 2R} \!-\! \partial_R u^3, \, -{\partial_\varphi u^1\over 2R}, \, -u^4\!+\!\partial_R u^1, \, {\delta\over\varkappa}u^1\!-\!u^3 \right) .\\ \end{align}\] Among these first-order adjoint symmetries, only \(\vec{\psi}_5\) satisfies the criterion (63 ) and therefore defines a conservation law. In this case, \[\begin{align} l_{\vec{\psi}_5} &=& -\square = -\square^\star \nonumber\\ &=& \begin{pmatrix} 0& (2R)^{-1} D_\varphi & -D_R & \delta/\varkappa\\ -(2R)^{-1}D_\varphi &0 & 0 & 0\\ D_R & 0 & 0 & -1\\ \delta/\varkappa& 0& -1 & 0 \end{pmatrix}. \end{align}\] The associated conserved current, written in terms of \(\vec{w}=(w_f,w_x,w_y,w_z)\) and the Cartesian Fourier variables \((s,k)\), has components \[\begin{align} A^t &=& w_f \left[ {\delta\over\varkappa}w_z + 2\mathrm{i}\left(\partial_s w_y-\partial_k w_x\right) \right] + \mathrm{i}w_z(sw_y-kw_x), \nonumber\\ A^s &=& -\varkappa \left[ k\left(w_f^2+w_z^2-w_x^2+w_y^2\right) + 2s w_xw_y \right] -kA^t, \\ A^k &=& \varkappa \left[ s\left(w_f^2+w_z^2+w_x^2-w_y^2\right) + 2k w_xw_y \right] +sA^t. \end{align}\] Thus, the conservation law is \[D_tA^t+D_sA^s+D_kA^k=0.\] Equivalently, in polar variables, this equation can be written as \[\begin{align} D_t A^t + D_\varphi \left(A^t + \varkappa(w_f^2 + w_z^2)\right) &=& \nonumber \\ &&-\varkappa(r \sin2\varphi D_r + \cos 2\varphi D_\varphi)(w_x^2 - w_y^2) \nonumber \\ &&+ 2\varkappa(r\cos2\varphi D_r - \sin 2\varphi D_\varphi)(w_x w_y). \end{align}\]

For the invariant solutions (44 )–(45 ) associated with \(\Theta_1\), one obtains \(A^t=-u_n^3u_n^4\), and the conservation law reduces to \[\partial_t(u_n^3u_n^4) = g\sqrt{2}\,{1\over r}\partial_r(u_n^2u_n^3).\]

Higher-order conservation laws can be constructed by applying the recursion operators \(\mathrm{D}_t\) and \(\mathrm{D}_\varphi\) to the co-symmetries. At second order, the following co-symmetries define conservation laws: \[\begin{align} \vec{\psi}_1^{(t)} &:=& \mathrm{D}_t\vec{\psi}_1 = \left( \partial_{tt}^2u^1, {\partial_{tt}^2u^2\over 4R}, {\partial_{tt}^2u^3\over 4R}, \partial_{tt}^2u^4 \right) , \\ \vec{\psi}_1^{(\varphi)} &:=& \mathrm{D}_\varphi\vec{\psi}_1 = \left( \partial_{t\varphi}^2u^1, {\partial_{t\varphi}^2u^2\over 4R}, {\partial_{t\varphi}^2u^3\over 4R}, \partial_{t\varphi}^2u^4 \right) , \\\nonumber \vec{\psi}_2^{(\varphi)} &:=& \mathrm{D}_\varphi\vec{\psi}_2 = \left( \partial_{\varphi\varphi}^2u^1, {\partial_{\varphi\varphi}^2u^2\over 4R}, {\partial_{\varphi\varphi}^2u^3\over 4R}, \partial_{\varphi\varphi}^2u^4 \right) .\\ \end{align}\] The co-symmetries obtained from \(\mathrm{D}\vec{\psi}_5\) (\(\mathrm{D}\) denotes either \(\mathrm{D}_t\) or \(\mathrm{D}_\varphi\)) do not satisfy the criterion (63 ) at this order.

The corresponding conserved currents have the same quadratic structure as the zero-order current. For \(\vec{\psi}_1^{(t)}\), one obtains \[\begin{align} A^t &=& {(\partial_tu^1)^2+(\partial_tu^4)^2\over 2} + {(\partial_tu^2)^2+(\partial_tu^3)^2\over 8R}, \\ A^R &=& -\varkappa\,\partial_tu^2\,\partial_tu^4, \\ A^\varphi &=& A^t - {\varkappa\over 2R}\partial_tu^3\,\partial_tu^4. \end{align}\] For \(\vec{\psi}_2^{(\varphi)}\), one obtains \[\begin{align} A^t &=& {(\partial_\varphi u^1)^2+(\partial_\varphi u^4)^2\over 2} + {(\partial_\varphi u^2)^2+(\partial_\varphi u^3)^2\over 8R}, \\ A^R &=& -\varkappa\,\partial_\varphi u^2\,\partial_\varphi u^4, \\ A^\varphi &=& A^t - {\varkappa\over 2R}\partial_\varphi u^3\,\partial_\varphi u^4. \end{align}\] Finally, for \(-2\vec{\psi}_1^{(\varphi)}\), the conserved current is \[\begin{align} \nonumber A^t &=& \partial_tu^1\,\partial_\varphi u^1 + \partial_tu^4\,\partial_\varphi u^4 + {\partial_tu^2\,\partial_\varphi u^2 + \partial_tu^3\,\partial_\varphi u^3 \over 4R},\\ && \\ A^R &=& -\varkappa \left( \partial_tu^2\,\partial_\varphi u^4 + \partial_\varphi u^2\,\partial_tu^4 \right), \\ A^\varphi &=& A^t - {\varkappa\over 2R} \left( \partial_tu^3\,\partial_\varphi u^4 + \partial_\varphi u^3\,\partial_tu^4 \right). \end{align}\]

For the invariant solution (44 )–(45 ), only the first of these second-order fluxes is nontrivial. It yields a conservation law analogous to Eq. (69 ): \[\begin{align} \partial_t \left[ (\partial_tu_n^1)^2 + (\partial_tu_n^4)^2 + \left({\partial_tu_n^2\over r}\right)^2 + \left({\partial_tu_n^3\over r}\right)^2 \right] =\nonumber \\ 4\varkappa\, {\partial_r\left(\partial_tu_n^2\,\partial_tu_n^4\right)\over r}. \end{align}\]

At third order, applying \(\mathrm{D}^2\) to \(\vec{\psi}_1\) and \(\vec{\psi}_2\) does not generate new conservation laws, whereas applying \(\mathrm{D}^2\) to \(\vec{\psi}_5\) does. For instance, the co-symmetry \(\mathrm{D}_t^2\vec{\psi}_5\) generates the current \[\begin{align} A^t &=& {\delta\over\varkappa}\partial_tu^1\,\partial_tu^4 - \partial_tu^3\,\partial_tu^4 + \partial_{tR}u^1\,\partial_tu^3 + {\partial_tu^1\,\partial_{t\varphi}u^2\over 2R}, \\ A^R &=& -\partial_tu^1\,\partial_{t\varphi}u^3 - \partial_tu^1\,\partial_{tt}u^3, \\ A^\varphi &=& A^t + \varkappa(\partial_tu^1)^2 + \varkappa(\partial_tu^4)^2 + {\varkappa\over 4R} \left[ (\partial_tu^2)^2+(\partial_tu^3)^2 \right] \nonumber\\ && - {\partial_{tt}u^1\,\partial_tu^2 + \partial_{t\varphi}u^1\,\partial_tu^2 \over 2R}. \end{align}\]

This pattern continues to higher orders. In particular, repeated applications of the recursion operators \(\mathrm{D}_t\) and \(\mathrm{D}_\varphi\) generate an infinite hierarchy of conservation laws. More precisely, \[\mathrm{D}^{2m+1}\vec{\psi}_{1,2}, \qquad \mathrm{D}^{2m}\vec{\psi}_5, \qquad m=0,1,\ldots,\] provide generating functions of conservation laws.

5 Summary↩︎

In this work, we have applied Lie symmetry theory to the system of partial differential equations describing the Jaynes–Cummings model (JCM) in the characteristic-function representation. The purpose of this analysis was to characterize the model from a symmetry-based perspective, identify its admitted invariant solutions, and investigate the conservation laws encoded in its differential formulation. The physical admissibility of the resulting solutions is determined by the quantum-mechanical conditions imposed on the field characteristic function and on the reduced atomic density matrix, which act as boundary and regularity conditions for the PDE system.

The symmetry analysis identifies two relevant classes of invariant solutions compatible with these requirements. The first class recovers the familiar dressed-state dynamics from the invariant reduction of the differential equations, thereby providing a different characterization of the standard JCM solution. The second class has a radial dependence expressed in terms of Heun polynomials and describes generic coupled atom–field configurations. This solution lies outside the conventional dressed-state construction and represents an additional invariant structure admitted by the PDE system, although its complete physical interpretation remains to be established. Other mathematically admissible invariant forms were found to be incompatible with the normalization and regularity conditions required of a physical characteristic function.

The differential system also possesses a nontrivial conservation-law structure. Using the methodology developed by Vinogradov, we recovered the well-known conservation of the total number of excitations directly from the PDE formulation. In addition, we derived conserved currents involving atomic populations, coherence, reduced-state purity, and moments of the field characteristic functions. Of particular interest is the balance equation for the quantity \(P_a-f_a^2\), where \(P_a\) is the purity of the reduced atomic state and \(f_a\) is its coherence function. Its rate of change is determined by the atom–field coupling and by combinations of atomic expectation values with field-related characteristic-function moments.

These results should not be interpreted as establishing a conservation law for entanglement itself. Rather, they show that the differential equations contain conserved currents involving quantities whose evolution is linked to the generation and redistribution of atom–field correlations. For globally pure states, variations of the reduced atomic purity are directly associated with variations of the bipartite atom–field entanglement. The derived balance laws therefore provide a possible mathematical description of how purity, coherence, and correlations are redistributed during the exchange of excitations between the two subsystems.

The JCM differential system also exhibits a rich mathematical structure. In addition to point symmetries, it admits a first-order generalized symmetry. The existence of recursion operators allows generalized symmetries of arbitrarily high order to be constructed. It is a remarkable property, since relatively few systems admit higher-order symmetries. Although the system is not of Euler–Lagrange type and Noether’s theorem cannot be applied directly, a mapping between symmetries and co-symmetries can be established. This mapping facilitates the construction of generating functions for conservation laws and leads to an infinite hierarchy of conserved currents.

Overall, the application of Lie symmetry analysis and Vinogradov’s conservation-law framework provides a complementary perspective on the JCM. It recovers known features of the model from the structure of its differential equations, identifies an additional class of invariant solutions, and reveals conservation laws that are not immediately apparent from the usual Hamiltonian treatment. These results suggest that symmetry and conservation-law methods may also be useful for investigating more complex composite quantum systems whose dynamics can be formulated as coupled systems of differential equations.

Data Availability:

All data supporting the findings of this study are contained within the manuscript.

6 The characteristic function↩︎

The characteristic function \(w(\vec{r},t)\) (or chord function)  [11] is defined as the Fourier transform of the Wigner function [20] \[\begin{align} \label{Wtransf} w(\vec{r},t)&=& \iint\mathrm{d}p\, \mathrm{d}q\; \mathrm{e}^{\mathrm{i}q k+\mathrm{i}s p}\, W(q,p,t) \\\nonumber &=& \int\mathrm{d}q\; \mathrm{e}^{\mathrm{i}q k}\, \langle q+s/2 |\, \varrho(t)\, | q-s/2\rangle\; . \end{align}\tag{73}\] This type of transformation maps the quantum mechanical phase-space representation \(\{q,p\}\), where the quasi-probability distribution function, or Wigner function, \({W}(q,p, t)\) is defined, into its Fourier image \(\{k,s\}\). Additionally, the application of position and momentum operators to the density operator can be readily translated into the characteristic function representation as (for convenience, higher powers of \(\hat{x}\) and \(\hat{p}\) are included): \[\begin{align} \tag{74} \hat{x}^n \hat{p}^m\, \varrho &\mapsto& \Big( \frac{s}{2} - \mathrm{i}\partial_k\Big)^{n} \Big( \frac{-k}{2} - \mathrm{i}\partial_s\Big)^{m}w\\ \varrho\, \hat{x}^{n}\hat{p}^{m} &\mapsto& \Big( \frac{-s}{2} - \mathrm{i}\partial_k \Big)^n \Big( \frac{k}{2} - \mathrm{i}\partial_s \Big)^m\; w\\ \hat{x}^{n}\,\varrho\, \hat{p}^{m} &\mapsto& \Big( \frac{s}{2} - \mathrm{i}\partial_k\Big)^{n} \Big( \frac{k}{2} - \mathrm{i}\partial_s\Big)^{m}w\\\tag{75} \hat{p}^{m}\,\varrho\, \hat{x}^{n} &\mapsto& \Big( \frac{-s}{2} - \mathrm{i}\partial_k\Big)^n \Big( \frac{-k}{2} - \mathrm{i}\partial_s \Big)^m\! w\, . \end{align}\] This also allows to obtain explicit expressions for the \(n\)-th order moments of products of position and momentum operators, by taking the appropriate partial derivatives at the origin of the coordinate system: \[\begin{align} \label{xpn} \frac{\langle\hat{x}^n\,\hat{p}^m \rangle+ \langle\hat{p}^m\,\hat{x}^n \rangle}{2} = (-\mathrm{i})^{n+m} \, \partial^n_k\, \partial^{m}_s \, w \big|_{k,s = 0}\; , \end{align}\tag{76}\] which also include the cases: \(\langle\hat{x}^{n}\rangle= (-\mathrm{i}\, \partial_k)^n\, w \big|_{k,s = 0}\), and \(\langle \hat{p}^{n}\rangle= (-\mathrm{i}\,\partial_s)^n\, w \big|_{k,s = 0}\).

For any characteristic function to represent a valid quantum state of a unique density operator \(\varrho\), it must satisfy some important properties:

• the normalization condition \(w(\vec{r}=0,t) = 1\),

• the Hermiticity condition: \(w(-\vec{r},t) =w^\ast(\vec{r},t)\),

• it has to be continuous and bounded and

• it has to satisfy the so called twisted positive definiteness (quantum Bochner condition) [27]: \[\label{twisted95pos} \int_{\mathbb{R}^4} \mathrm{d}\vec{r}^{\,2}\mathrm{d}\vec{r}\,'^{\,2} f^*(\vec{r}) \mathrm{e}^{\mathrm{i}\sigma(\vec{r},\vec{r}\,')}\, w_f(\vec{r}-\vec{r}\,',t)f^*(\vec{r}\,') > 0.\tag{77}\] where \[\sigma(\vec{r},\vec{r}\,') = \vec{r} \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}\vec{r}\,',\] and \(f(\vec{r})\in L^2(\mathbb{R}^2)\) is any quadratic integrable test function. This condition is very rigorous and is the quantum analog of the classical Bochner condition for characteristic functions, which states that the characteristic function of a classical probability distribution must be positive definite. In the quantum case, the twisted positive definiteness condition ensures that the characteristic function corresponds to a valid quantum state, which is represented by a positive semi-definite density operator: \(\langle \psi | \varrho | \psi \rangle \geq 0\). If these properties are fulfilled, then the characteristic function represents a valid quantum state.

7 Revisiting the Jaynes-Cummings model↩︎

The standard Jaynes-Cummings model (JCM) describes the interaction between a two-level atom and a single quantized mode of the electromagnetic field through the exchange of energy quanta. Its usual analytical treatment relies on the existence of a constant of motion associated with the total number of excitations in the composite atom-field system. In a dimensionless formulation defined with respect to the time scale of the field dynamics, this constant of motion is represented by the total excitation-number operator \[\hat{N}=\frac{\Delta}{2}\sigma_z\otimes\openone+\openone\otimes\hat{n},\] where \(\Delta=\omega_q/\omega_o\), \(\omega_q\) is the atomic level spacing measured in units of frequency, \(\omega_o\) is the characteristic frequency of the field mode, also described in the same units and \(\sigma_z\) and \(\hat{n}\) denote, respectively, the atomic inversion and the field number operator.

The eigenstates of \(\hat{N}\) define invariant excitation subspaces of the atom-field Hilbert space. In particular, each subspace is spanned by the bare states \[\{|e,n\rangle, |g,n+1\rangle\}\subset \mathcal{H}_a\otimes\mathcal{H}_f,\] where \(|e\rangle\) and \(|g\rangle\) denote the excited and ground states of the two-level atom, while \(|n\rangle\) represents a number state of the field mode.

Since the total number of excitations is conserved, the operator \(\hat{N}\) commutes with the Hamiltonian, \([\hat{H},\hat{N}]=0\). Therefore, the dynamics decomposes into invariant two-dimensional subspaces, and in the basis \(\{|e,n\rangle,|g,n+1\rangle\}\) the Hamiltonian is represented as \[\label{Hn} \hat{H}\left(\begin{array} {c} |e,n\rangle \\ |g,n+1\rangle \end{array} \right) = \left(\begin{array} {cc} {\Delta\over 2} + n & g\sqrt{n+1} \\ g\sqrt{n+1} & -{\Delta\over 2} + n + 1 \end{array} \right) \left(\begin{array} {c} |e,n\rangle \\ |g,n+1\rangle \end{array} \right).\tag{78}\]

The dressed states, namely the eigenstates of the JCM Hamiltonian in each excitation subspace, are obtained by diagonalizing this matrix. Thus, \[\begin{align} \label{Heigen} \hat{H}\left(\begin{array} {c} |+\rangle \\ |-\rangle \end{array} \right) &=& \left(\begin{array} {cc} \varepsilon_+ & 0 \\ 0 & \varepsilon_- \end{array} \right) \left(\begin{array} {c} |+\rangle \\ |-\rangle \end{array} \right), \end{align}\tag{79}\] where the corresponding eigenvalues are \(\varepsilon_+ = \varepsilon_n-\lambda_n/2\) and \(\varepsilon_- = \varepsilon_n+\lambda_n/2\), with \(\varepsilon_n=n+1/2\). The quantity \[\label{lambn} \lambda_n = \sqrt{\delta^2 + 4g^2(n+1)}\tag{80}\] is the generalized Rabi frequency, where \(\delta=\Delta-1\) is the detuning. The resonant regime corresponds to \(\delta=0\), for which the atomic transition frequency coincides with the field-mode frequency.

The dressed states are related to the bare states through the transformation \[\label{vtrans} \left(\begin{array} {c} |+\rangle \\ |-\rangle \end{array} \right) = \left(\begin{array} {cc} \cos(\theta_n) & \sin(\theta_n) \\ -\sin(\theta_n) & \cos(\theta_n) \end{array} \right) \left(\begin{array} {c} |e,n\rangle \\ |g,n+1\rangle \end{array} \right),\tag{81}\] where \[\theta_n = \arctan\left(\frac{\lambda_n-\delta}{2g\sqrt{n+1}}\right)\] is the mixing angle between the bare atom-field states. Equivalently, \(\tan(2\theta_n)=(2g\sqrt{n+1})/\delta,\) while \[\begin{align} \label{thetan} \cos\theta_n=\sqrt{\frac{\lambda_n+\delta}{2\lambda_n}}, &\quad& \sin\theta_n=\sqrt{\frac{\lambda_n-\delta}{2\lambda_n}}. \end{align}\tag{82}\]

Once the dressed states and their eigenvalues are known, the evolution of any pure state within a fixed excitation subspace can be expressed in this basis. For an initial state: \(|\psi\rangle = c_+|+\rangle + c_-|-\rangle\), with \(c_+,c_-\in\mathbb{C}\) and \(|c_+|^2+|c_-|^2=1\), the time evolution is given by \[\label{soldress} |\psi(t)\rangle = c_+ \mathrm{e}^{\mathrm{i}\varepsilon_+ t}|+\rangle + c_- \mathrm{e}^{\mathrm{i}\varepsilon_- t}|-\rangle.\tag{83}\]

Removing the global phase factor \(\mathrm{e}^{\mathrm{i}\varepsilon_n t}\), this state can be written in the bare basis as \[\label{psit} |\tilde{\psi}(t)\rangle = a_n(t)|e,n\rangle + b_n(t)|g,n+1\rangle,\tag{84}\] where \(|\tilde{\psi}(t)\rangle=\mathrm{e}^{-\mathrm{i}\varepsilon_n t}|\psi(t)\rangle\), and \[\begin{align} \label{atbt} a_n(t) &=& c_+\cos\theta_n\,\mathrm{e}^{\mathrm{i}\lambda_n t/2} - c_-\sin\theta_n\,\mathrm{e}^{-\mathrm{i}\lambda_n t/2}, \nonumber \\ b_n(t) &=& c_+\sin\theta_n\,\mathrm{e}^{\mathrm{i}\lambda_n t/2} + c_-\cos\theta_n\,\mathrm{e}^{-\mathrm{i}\lambda_n t/2}. \end{align}\tag{85}\] These amplitudes satisfy the normalization condition \(|a_n(t)|^2+|b_n(t)|^2=1\) for all times.

The dressed-state representation provides the standard analytical route for solving the JCM in each excitation subspace. More general initial conditions, however, may involve superpositions over many photon-number sectors, mixed states, or field-state representations that are not most conveniently handled by working directly in the dressed-state basis. In such cases, the system evolution can be described through the density operator, \[\varrho(t)=\hat{U}(t)\varrho(0)\hat{U}^{\dagger}(t),\] where \(\hat{U}(t)=\exp(-\mathrm{i}\hat{H}t)\) is the evolution operator. Since the Hamiltonian is block-diagonal in the excitation subspaces, the evolution operator can be written as \[\hat{U}(t)= \mathrm{e}^{-\mathrm{i}\Delta t}|g,0\rangle\langle g,0| + \sum_{n=0}^{\infty}\mathrm{e}^{-\mathrm{i}H_n t}\Pi_n,\] where \[\Pi_n = |e,n\rangle\langle e,n| + |g,n+1\rangle\langle g,n+1|\] is the projector onto the \(n\)-th excitation subspace, and \(H_n\) denotes the matrix representation in Eq. (78 ). Thus, the standard construction requires the diagonalization of each excitation block.

Alternatively, the evolution operator may be obtained through direct matrix-exponentiation methods, yielding a field-number-operator-dependent representation [28], [29]: \[\hat{U}(t) = \mathrm{e}^{-\mathrm{i}\hat{n} t} \left(\begin{array} {cc} \cos (\hat{\lambda} t) - \mathrm{i}\hat{\gamma}\sin ( \hat{\lambda} t ) & -\hat{\Omega}\sin ( \hat{\lambda} t ) \\ \mathrm{i}\hat{\Omega} \sin ( \hat{\lambda} t ) & \cos (\hat{\lambda} t) + \mathrm{i}\hat{\gamma} \sin ( \hat{\lambda} t ) \end{array} \right),\] where \[\hat{\lambda}=\frac{1}{2}\sqrt{\Delta^2+4g^2(\hat{n}+1)}, \quad \hat{\gamma}=\frac{\Delta}{\hat{\lambda}}, \quad \hat{\Omega}=\frac{2g\sqrt{\hat{n}+1}}{\hat{\lambda}}.\] This representation can be applied directly to field states expanded in the number-state basis.

7.1 Dressed state solution of the JCM in the characteristic function representation↩︎

The solution of the JCM in the dressed and the bare state basis is given by (83 ,84 ). Transformation of this solution into the characteristic function representation is achieved by applying the transformation (73 ) to each of the components of the density operator. The density operator in the bare state basis is given by: \[\begin{align} \varrho(t) &=& |a_n(t)|^2 | e,n \rangle \langle e,n | + |b_n(t)|^2 | g,n+1 \rangle \langle g,n+1 | \\\nonumber &+& a_n(t) b_n^*(t) | e,n \rangle \langle g,n+1 | + b_n(t) a_n^*(t) | g,n+1 \rangle \langle e,n |\,, \end{align}\] thus the components of the density operator can be identified as: \[\begin{align} \varrho_{ee}(t) &=& |a_n(t)|^2 | n \rangle \langle n |,\\ \varrho_{gg}(t) &=& |b_n(t)|^2 | n+1 \rangle \langle n+1 |,\\ \varrho_{eg}(t) &=& a_n(t) b_n^*(t) | n \rangle \langle n+1 |,\\ \varrho_{ge}(t) &=& b_n(t) a_n^*(t) | n+1 \rangle \langle n |\,. \end{align}\] Using the definition of the characteristic function (73 ) and the position representation of the number states: \[\psi_{n}(x)=\langle x | n \rangle = \frac{1}{\sqrt{2^n n! \sqrt{\pi}}} \, H_n\left( x \right) \mathrm{e}^{- x^2/2}\,,\] where \(H_n(x)\) are the Hermite polynomials, we obtain the the dressed state density matrix solution of the JCM: \[\varrho(\vec{r},t) = \begin{pmatrix} w^{(n)}_{ee}(\vec{r},t) & w^{(n)}_{eg}(\vec{r},t) \\ w^{(n)}_{ge}(\vec{r},t) & w^{(n)}_{gg}(\vec{r},t) \end{pmatrix}\,,\] with the density matrix components described in terms of the field characteristic function as: \[\begin{align} w^{(n)}_{ee}(\vec{r},t\,) &=& |a_n(t)|^2 \, \chi_{ee}(r),\tag{86}\\ w^{(n)}_{gg}(\vec{r},t\,) &=& |b_n(t)|^2 \, \chi_{gg}(r),\\ w^{(n)}_{eg}(\vec{r},t\,) &=& a_n(t) b_n^*(t) \, \chi_{eg}(r),\\ w^{(n)}_{ge}(\vec{r},t\,) &=& b_n(t) a_n^*(t) \, \chi_{ge}(r), \tag{87} \end{align}\] where \[\begin{align} \tag{88} \chi_{ee}(r) &=& \mathrm{e}^{-r^2/4}\, L_n(r^2/2) , \\ \tag{89} \chi_{gg}(r) &=& \mathrm{e}^{-r^2/4}\, L_{n+1}(r^2/2), \\ \tag{90} \chi_{eg}(r,\varphi) &=& - \frac{ \mathrm{e}^{- \mathrm{i}\varphi}}{\sqrt{2(n+1)}} \, { r \, \mathrm{e}^{-r^2/4} } \, L_{n}^{(1)}(r^2/2), \\ \tag{91} \chi_{ge}(r,\varphi) &=& \frac{ \mathrm{e}^{\mathrm{i}\varphi}}{\sqrt{2(n+1)}} \, { r \, \mathrm{e}^{-r^2/4}}\, L^{(1)}_{n}(r^2/2). \end{align}\] where \(L_n(x)\) and \(L_n^{m}(x)\) are the Laguerre and associated Laguerre polynomials, respectively, and \(r^2 = k^2 + s^2\). Here we recall that \(\vec{r} = (k,s) = (r\sin\varphi,r \cos\varphi)\) represents a position vector in the Fourier phase-space. Within this representation, the field dynamics can be obtained by tracing out the atom degrees of freedom, yielding: \[\begin{align} \nonumber w^{(n)}_f(r,t) &=& w^{(n)}_{ee}(r,t) + w^{(n)}_{gg}(r,t)\\\nonumber & =& |a_n(t)|^2\mathrm{e}^{-r^2/4}\, L_n(r^2/2) \\\label{w95f95CHI} &&+ |b_n(t)|^2\mathrm{e}^{-r^2/4}\, L_{n+1}(r^2/2)\,, \end{align}\tag{92}\] while the characteristic function description of the atom population in this representation, i.e.: \(w^{(n)}_z(\vec{r},t) = \mathrm{tr}_a\{\sigma_z \varrho(\vec{r},t)\}\) are given by: \[\begin{align} \nonumber w^{(n)}_z(r,t) &=& w^{(n)}_{ee}(r,t) - w^{(n)}_{gg}(r,t)\\\nonumber & =& |a_n(t)|^2\mathrm{e}^{-r^2/4}\, L_n(r^2/2) \\ &&- |b_n(t)|^2\mathrm{e}^{-r^2/4}\, L_{n+1}(r^2/2). \end{align}\]

By considering \[\begin{align} \nonumber a_n(t) &=& c_+ \cos\theta_n\, \mathrm{e}^{-\mathrm{i}t \lambda_n/2} - c_- \sin\theta_n\, \mathrm{e}^{\mathrm{i}t \lambda_n/2}, \\ b_n(t) &=& c_+ \sin\theta_n\, \mathrm{e}^{-\mathrm{i}t \lambda_n/2} + c_- \cos\theta_n\, \mathrm{e}^{\mathrm{i}t \lambda_n/2}, \end{align}\] with \(\cos\theta_n\) and \(\sin\theta_n\) defined as in (82 ), and by taking \[c_{+,-} = \sqrt{\frac{1}{2} \mp \frac{\delta}{2\lambda_n}},\] we have \[\begin{align} a_n(t) &=& -\mathrm{i}\frac{2g\sqrt{n + 1}}{\lambda_n} \sin\frac{\lambda_n t}{2}, \nonumber \\ b_n(t) &=& \cos\frac{\lambda_n t}{2} + \mathrm{i}\frac{\delta}{\lambda_n} \sin\frac{\lambda_n t}{2}. \end{align}\] Thus, for example, Eq. (92 ) converts to the following expression \[w^{(n)}_f(\vec{r},t) = g^2 \frac{1 - \cos\lambda_n t}{\lambda_n^2} Y_n(r) + \mathrm{e}^{-r^2/4} L_{n+1}(r^2/2).\]

After transforming the characteristic-function components (86 )–(87 ) to the variables \(\vec{u}_\chi(\vec{r},t)\), we find that the resulting solution is related to the invariant solution (44 )–(45 ) through the superposition principle of the linear system. More precisely, the two solutions coincide after adding to \(\vec{u}_\chi\) the stationary solution \[\vec{u}_{\mathrm{st}}(\vec{r}) \!=\! \left( \mathrm{e}^{-{r^2\over 4}}L_{n+1}(r^2/2) \!+\!\frac{2Y_n(r)}{r^2}, \, 0, \, 0, \, \!-\mathrm{e}^{-{r^2\over 4}}L_{n+1}(r^2/2) \right) .\] Thus, \[\vec{u}(\vec{r},t) = \vec{u}_\chi(\vec{r},t) + \vec{u}_{\mathrm{st}}(\vec{r}),\] with the constant chosen as \[c=-\frac{g\sqrt{2}}{\lambda_n}.\] With this choice, \(\vec{u}(\vec{r},t)\) reproduces the invariant solution given in Eqs. (44 )–(45 ).

References↩︎