Some mathematical problems in a neoclassical theory of electric charges

We study here a number of mathematical problems related to our recently introduced neoclassical theory for the electromagnetic phenomena in which charges are represented by complex valued wave functions as in the Schrodinger wave mechanics. Dynamics of charges in the non-relativistic case is governed by a system of nonlinear Schrodinger equations coupled with the electromagnetic fields, and we prove for it that the centers of wave functions converge in macroscopic regimes to trajectories of points governed by the Newton's equations with the Lorentz forces provided the wave functions remain localized. Exact solutions in the form of localized accelerating solitons are found. Our studies of a class of time multiharmonic solutions of the same field equations show that they satisfy Planck-Einstein relation and that the energy levels of the nonlinear eigenvalue problem for the hydrogen atom converge to the well-known energy levels of the linear Schrodinger operator when the free charge size is much larger than the Bohr radius.


Introduction
It is well known that the concept of a point charge interacting with the electromagnetic (EM) field has fundamental problems. Indeed, in the classical electrodynamics the evolution of a point charge q of a mass m in an external electromagnetic (EM) field is governed by Newton's equation where r and v =ṙ = dr dt are respectively the charge's position and velocity, E (t, r) and B (t, r) are the electric field and the magnetic induction, and the right-hand side of the equation (1) is the Lorentz force. On other hand, if the charge's timedependent position and velocity are respectively r and v then the corresponding EM field is described by the Maxwell equations 1 c where δ is the Dirac delta-function, v (t) =ṙ (t), c is speed of light. If one would like to consider equations (1)-(3) as a closed system "charge-EM field" there is a problem. Its origin is a singularity of the EM field exactly at the position of the point charge, as, for instance, for the electrostatic field E with the Coulomb's potential q |x−r| with a singularity at x = r. If one wants to stay within the classical electromagentic theory, a possible remedy is the introduction of an extended charge which, though very small, is not a point. There are two most well known models for such an extended charge: the semi-relativistic Abraham rigid charge model (a rigid sphere with spherically symmetric charge distribution), [35,Sections 2.4,4.1,10.2,13], [32,Sections 2.2], and the Lorentz relativistically covariant model which was studied and advanced in [1], [20,Sections 16], [28], [30], [32,Sections 2,6], [33], [35,Sections 2.5,4.2,10.1], [37]. Poincaré suggested in 1905Poincaré suggested in -1906, [31] (see also [20,, [32,Sections 2.3,, [29, Section 63], [33], [37,Section 4.2] and references there in), to add to the Lorentz-Abraham model non-electromagnetic cohesive forces which balance the charge internal repulsive electromagnetic forces and remarkably restore also the covariance of the entire model.
Another well known problem of a physical nature with point charges and the classical electrodynamics is related to the Rutherford planetary model for the hydrogen atom. The planetary model of the hydrogen atom is inconsistent with physically observed phenomena such as the stability of atoms and the discreteness of the hydrogen energy levels. That inconsistency lead as we well know to introduction of non-classical models, such as the Bohr model and later on the Schrödinger's model of the hydrogen atom.
To address the above mentioned problems with the classical electrodynamics we introduced recently in [5], [6] wave-corpuscle mechanics (WCM) for charges which is conceived as one mechanics for macroscopic and atomic scales. This model is neoclassical in the sense that it is based on the classical concept of the electromagnetic field but an elementary charge is not a point but a wave-corpuscle described by a complex-valued scalar wave function ψ with the density |ψ (t, x)| 2 which is not given a probabilistic interpretation. The dynamics and shape of a wave corpuscle is governed by a nonlinear Klein-Gordon or a nonlinear Schrödinger equation in relativistic and nonrelativistic cases respectively, and a wave-corpuscle is defined as a special type of solutions to these equations. Consequently an elementary charge in the WCM does not have a prescribed geometry which differs it from the classical Abraham and Lorentz models and more recent developments (see [7], [23], [25], [18], [35] and references therein). Another important difference is that Newton's equations are not postulated as in (1) but rather are derived from the field equations in a non-relativistic regime when charges are well separated and localized. Physical aspects of the model are discussed in detail in [5], [6], where the Poincaré type cohesive forces are constructed somewhat differently then we propose here. The primary focus of this paper is on studies of mathematical properties of the field equations of the WCM.
The paper is organized as follows. In the first two subsections of this Introduction we describe the system relativistic and non-relativistic Lagrangians and the corresponding field equations, and in the third one we consider the relevant nonlinearities and their examples. Importantly, the nonlinearity depends on a size parameter a > 0 associated with a free solution (ground state). Section 2 is devoted to remote interaction regimes for many charges. We prove there that in the case of the non-relativistic field equations when a → 0 (macroscopic limit) the centers of the interacting charges defined by the formula r ℓ (t) = R 3 x ψ ℓ (t, x) 2 dx,ℓ = 1, ..., N, converge to the solutions of Newton's equations with the Lorentz forces if ψ ℓ remain localized. We also provide examples of exact solutions of the field equations in the form of accelerating solitons for which the localization assumption holds. In Section 3 we study systems of bound charges and time multiharmonic solutions to the field equations. We consider there, in particular, a connection, discovered in a different setting in [12], between the Planck-Einstein energy-frequency relation and the logarithmic nonlinearity as well as some dynamical issues related to the logarithmic nonlinearity closely related to results in [14], [16]. We continue then with a study of a system of two charges with the logarithmic nonlinearity as a model for the hydrogen atom. We prove that if κ = a 1 /a → 0, where a 1 is the Bohr radius, then the lower energy levels associated with this model converge to the well known energy levels of the linear Schrödinger operator for the hydrogen atom. The proof is based on an approach developed in [9], [10] with certain modifications.
1.1. Relativistic Lagrangian and field equations. Let us consider a system of N charges interacting directly only with the EM field described by its 4-vector potential (ϕ, A). The charges are described by their wave functions ψ ℓ with the superscript index ℓ = 1, . . . , N labeling them. The EM fields are related to the potentials ϕ, A by the following standard relations We introduce now the system Lagrangian L which is though similar to the one introduced in [5], [6] but differs from it. A different Lagrangian is introduced here in order to get the hydrogen atom frequency spectrum with desired precision. Namely, to every ℓ-th charge is assigned an adjunct potential ϕ ℓ , A ℓ describing its additional degrees of freedom and our relativistic Lagrangian is defined by where (i) the covariant derivatives are defined bỹ (ii) ψ * is complex conjugate to ψ; (iii) m ℓ > 0 is the ℓ-th charge mass, q ℓ is the value of the charge, κ 0ℓ = m ℓ c χ , and χ > 0 is a constant similar to the Planck constant = h 2π ; (iv) G ℓ ψ ℓ * ψ ℓ is a nonlinearity which will be described below.
where ρ ℓ , J ℓ are defined by (9), (10). Based on the above equations we assume that always Hence, from now on we assume that in the relativistic case the dynamics of EM fields and charges is determined by equations (17), (18), (9), (10), (13) and according to (19) and (6) the covariant derivatives∂ ℓ t and∇ ℓ in (13) involve only fields ϕ ℓ ′ , implying that the adjunct potentials completely compensate the EM self-action for every charge and effectively there is no EM self-interaction (in [5], [6] this kind of compensation was gained by using an additional nonlinear self-interaction). Let us turn now to the nonlinearities G. We introduce the nonlinearities here the same way as in [5], [6] with an intension to have a rest (ground) state for a single charge in the form of a localized wave function (form factor). Notice that for a single charge evidently N = 1, ϕ − ϕ ℓ = 0, A − A ℓ = 0 and we look for the rest solution in the form with A ℓ = 0. Substituting (21) in (13) yields Now we choose a strictly positive, monotonically decreasing radial functionψ (ground state) as a parameter for the model and determine the nonlinearity G ′ from the following charge equilibrium condition: The above equation allows to determine the nonlinearity G ′ as long asψ is a strictly positive, smooth and monotonically decreasing function of r = |x|. In Section 1.3 we consider the nonlinearity in more details providing also examples. Note that using the Lorentz invariance of the system one can easily obtain a solution which represents the charge-field moving with a constant velocity v simply by applying to the rest solution ψ ℓ , ϕ, 0 the Lorentz transformation (see [5], [6]).
1.2. Non-relativistic Lagrangian and field equations. Our non-relativistic model describes the case of charges moving with non-relativistic velocities and it is set as follows. Using the frequency-shifting substitution (21) with more general ψ = ψ ℓ ω (t, x) which depends on (t, x) we observe that the second time derivative in (13) can be written in the form We neglect the term with the factor 1 c 2 and substitute −2i m ℓ (13). Consequently, we replace the nonlinear Klein-Gordon equation (13) by the following nonlinear Schrödinger equation (where we denote the frequency shifted ψ ℓ ω by ψ ℓ ) where∇ ℓ is given by (20). Since the magnetic fields generated by moving charges also have coefficient 1 c we neglect them preserving only the external magnetic fields and replace∇ ℓ by the covariant gradient∇ ℓ ex defined bỹ So our non-relativistic Lagrangian iŝ where A ex (t, x) and ϕ ex (t, x) are potentials of external EM fields, ψ ℓ * is complex conjugate of ψ ℓ . The Euler-Lagrange equations for the electrostatic potentials have the form Assuming that that ϕ, ϕ ℓ vanish at infinity we obtain a reduced version of (19) Similarly to the relativistic case we assume now that nonrelativistic equations for dynamics of charges in the external EM field with potentials ϕ ex , A ex take the form where ϕ =ℓ is given by (23) and ϕ ℓ is determined from the equations The solution of the above equation is given by the formula The nonlinear self-interaction terms G ℓ a in (30) are determined through the charge equilibrium equation (22) and index a indicates the dependence on the size parameter a > 0 which we introduce by the formula More detailed discussion of the nonlinearity is given in the following section.
1.3. Nonlinearity, its basic properties and examples. As we have already mentioned, the nonlinear self interaction function G is determined from the charge equilibrium equation (22) based on the form factor (free ground state)ψ. Important features of our nonlinearity include: (i) the boundedness or slow subcritical growth of its derivative G ′ (s) for s → ∞ with consequent boundedness from below of the energy; (ii) slightly singular behavior about s = 0, that is for small wave amplitudes.
In this section we consider the construction of the function G, study its properties and provide examples for which the construction of G is carried out explicitly. Throughout this section we have ψ,ψ ≥ 0 and hence |ψ| = ψ.
We introduce explicitly the dependence of the free ground stateψ on the size parameter a > 0 through the following representation of the functionψ (r) whereψ 1 (r) is a function of the dimensionless variable r ≥ 0. The dependence on a is chosen so that L 2 -norm ψ a (r) does not depend on a, hence the function ψ a (r) satisfies the charge normalization condition (53) for every a > 0. Obviously, definition (34) is consistent with (22) and (33). The size parameter a naturally has the dimension of length. A properly defined spatial size ofψ a , based, for instance, on the variance, is proportional to a with a coefficient depending onψ 1 . The charge equilibrium equation (22) can be written in the following form: The functionψ a (r) is assumed to be a smooth (at least twice continuously differentiable) positive monotonically decreasing function of r ≥ 0 which is square integrable with weight r 2 , we assume that its derivativeψ ′ a (r) is negative for r > 0 and we assume it to satisfy the charge normalization condition of the form (53); such a function is usually called in literature a ground state.
Let us look first at the case a = 1,ψ a =ψ 1 ,φ a =φ 1 , for which the equation (35) yields the following representation for G ′ (ψ 2 1 ) from (35) Sinceψ 2 1 (r) is a monotonic function, we can find its inverse r = r ψ 2 , yielding Sinceψ 1 (r) is smooth and ∂ rψ1 < 0, G ′ (|ψ| 2 ) is smooth for 0 < |ψ| 2 <ψ 2 1 (0). If we do not need G ′ (s) to be smooth, we extend G ′ (s) for s ≥ψ 2 1 (0) as a constant, namely The first derivative of such an extension at s =ψ 2 1 (0) has a discontinuity point. If ψ a (r) is a smooth function of class C n , n > 2, we always can define an extension of G ′ (s) for s ≥ψ 2 1 (0) as a bounded function of class C n−2 for all r > 0 and Slowly growing (subcritical) functions G ′ (s) which are not constant for large s also can be used, see examples below.
In the case of arbitrary size parameter a > 0 we define G ′ a (s) by formula (33), and this definition is consistent with (34) and (37).
Let us take a look at general properties of G ′ (s) as they follow from defining them relations (37). In the examples below the function G ′ (s) is not differentiable at s = 0, but ifψ (r) decays exponentially or with a power law the nonlinearity g (ψ) = G ′ (|ψ| 2 )ψ as it enters the field equation (30) is differentiable for all ψ including zero, hence it satisfies the Lipschitz condition. For a Gaussianψ 1 (r) which decays superexponentially G ′ (|ψ| 2 ) is unbounded at zero and g (ψ) is not differentiable at zero. Sinceψ (|x|) > 0, the sign of G ′ 1 |ψ| 2 coincides with the sign of ∇ 2ψ 1 (|x|). At the origin x = 0 the functionψ 1 (|x|) has its maximum and, consequently, G ′ 1 (s) ≤ 0 for s close to s =ψ 2 1 (0). The Laplacian applied to the radial functionψ 1 takes the form 1 r ∂ 2 ∂r 2 rψ 1 |x| . Consequently, if rψ 1 (r) is convex at r = |x| we have ∇ 2ψ 1 (|x|) ≥ 0. Since r 2ψ 1 (r) is integrable, we naturally assume that |x|ψ 1 (|x|) → 0 as |x| → ∞. Then if the second derivative of rψ 1 (r) has a constant sign near infinity, it must be non-negative as well as G ′ 1 (s) for s ≪ 1.
In the examples we give below G ′ 1 (s) has exactly one zero on the half-axis. Example 1. Consider a form factorψ 1 (r) decaying as a power law, namelẙ where c pw is the normalization factor, c pw = 3 1/2 / (4π) 1/2 . This function evidently is positive and monotonically decreasing. Let us find now G ′ (s) based on the relations (37 The extension for s ≥ c 2 pw can be defined as a constant or the same formula (41) can be used for all s ≥ 0 since corresponding nonlinearity has a subcritical growth.
If we explicitly introduce size parameter a into the form factor using (34), we define G ′ a (s) by (33). Notice that the variance of the form factorψ 2 1 (|x|) decaying as a power law (40) is infinite.
Example 2. Now we consider an exponentially decaying form factorψ 1 of the formψ Evidentlyψ 1 (r) is positive and monotonically decreasing. The dependence r (s) defined by the relation (42) is as follows: An elementary computation shows that Combining (43) with (44) we readily obtain the following function for s ≤ c 2 e e −2 We can extend it for larger s as follows: or we can use a smooth extension as in (39). The function G ′ 1 (s) is not differentiable at s = 0. At the same time the function g (ψ) = G ′ 1 (ψ (r)) ψ if we set g (0) = 0 is continuous and g (ψ) is continuously differentiable with respect to ψ at zero and g (ψ) satisfies a Lipschitz condition. The variance of the exponential form factor ψ 1 (r) is obviously finite. To find G ′ a (s) for arbitrary a we use its representation (33).
Example 3. Now we define a Gaussian form factor by the formulå Such a ground state is called gausson in [12]. Elementary computation shows that Hence, we define the nonlinearity by the formula and refer to it as logarithmic nonlinearity. The nonlinear potential function has the form Dependence on the size parameter a > 0 is given by the formula Obviously g (ψ) = G ′ 1 |ψ| 2 ψ is continuous for all ψ ∈ C if at zero we set g (0) = 0 and is differentiable for every ψ = 0 but is not differentiable at ψ = 0 and does not satisfy the Lipschitz condition. Notice also that g (ψ) has a subcritical growth as |ψ| → ∞.

Charges in remote interaction regimes
The primary focus of this section is to show that if the size parameter a → 0 then the dynamics of centers of localized solutions is approximated by the Newton equations with the Lorentz forces. This is done in the spirit of the well known in the quantum mechanics Ehrenfest Theorem, [34,Sections 7,23]. We also provide as an example explicit wave-corpuscle solutions which have such a dynamics.
As a first step we describe basic properties of the classical solutions of (30), (31). The LagrangianL is gauge invariant and every ℓ-th charge has a 4-current ρ ℓ , J ℓ defined by which satisfies the continuity equations ∂ t ρ ℓ + ∇ · J ℓ = 0 or Note that J ℓ defined by (51) agrees with the definition of the current (10) in the Maxwell equations. Equations (52) can be obtained by multiplying (30) by ψ ℓ * and taking the imaginary part. Integrating the continuity equation we see that ψ ℓ 2 = const and we impose the following normalization condition: The motivation for this particular normalization is that this normalization allows to converge to a delta-function and, in addition, for this normalization the solution to (31) given by the formula (32) converges to the Coulomb's potential with the value q ℓ of the charge if ψ ℓ a 2 converges to a delta function. The momentum density P ℓ for the LagrangianL 0 in (26) is defined by the formula Note that so defined momentum density P ℓ is related with the current J ℓ in (51) by the formula We introduce the total individual momenta P ℓ for ℓ-th charge by and obtain the following equations for the total individual momenta where The external EM fields E ex , B ex in (56) corresponding to the potentials ϕ ex , A ex are determined by standard formulas (4). Derivation of (56) is rather elementary. Indeed, in the simplest case where A ex = 0 we multiply (30) by ∇ψ ℓ * , take the real part and integrate the result over the entire space using integration by parts. To obtain (56) in more involved general case one can similarly multiply (30) by∇ ℓ * ψ ℓ * and then integrate the result by parts using some vector algebra manipulation.
2.1. Newtonian mechanics as an approximation. Let us show now that if the size parameter a is small compared to the macroscopic scale of variation of EM fields, the charge evolution can be described approximately by the Newton equations with the Lorentz forces similar to (1). In this subsection we give a heuristic derivation. Conditions under which this kind of derivation is justified are given in the next section.
We introduce the ℓ-th charge position r ℓ (t) and velocity v ℓ (t) as the following spatial averages where the current density J ℓ is defined by (51). We show below that a combination of the continuity equation (52) with the momentum evolution equations (56) imply the following remarkable property: the positions r ℓ (t) satisfy with a high accuracy Newton's equations of motion for the system of N point charges if the size parameter a is small enough. Multiplying continuity equation (52) by x and integrating we find the following identities showing that the positions and velocities defined by formulas (58) are related exactly as in the point charge mechanics. Then integrating (54) we obtain the following kinematic representation for the total momentum which also is exactly the same as for the point charges mechanics. Relations (59) and (60) yield and we obtain from (56) the following system of equations of motion for N charges: x), E ex and B ex are defined by (4).
The derivation of the above system is analogous to the well known in quantum mechanics Ehrenfest Theorem, [34,Sections 7,23], [12]. Now we give a formal derivation of Newton's law of motion. Let us suppose that for every ℓ-th charge density ψ ℓ 2 and the corresponding current density J ℓ are localized in a-vicinity of the position r ℓ (t), and that r ℓ (t) − r ℓ ′ (t) ≥ γ > 0 with γ independent on a on time interval [0, T ]. Then if a → 0 we get where the coefficients before the Dirac delta-functions are determined by the charge normalization conditions (53) and relations (58). Using potential representations (32) we infer from (63) the convergence of the potentials ϕ ℓ to the corresponding Coulomb's potentials, namely Hence, when passing to the limit as a → 0 we can recast the equations of motion (62) as the system where and ǫ 0 → 0 as a → 0. Notice that the terms f ℓ in equations (65) coincide with the Lorentz forces and we see that the limit equations of motion obtained from (65) coincide with Newton's equations of motion for point charges interacting via the Coulomb forces and with the external EM field via corresponding Lorentz forces, namely The above formal derivation of Newton's equations of motion is very simple, but if implemented with full rigor would require to impose rather high regularity conditions on the solution and rather strong requirements on the convergence to the delta function. Note that we use essentially the fact that the nonlinearity G ′ a , that singularly depends on a as a → 0 according to (33), does not enter the system (62). In the following subsection we rewrite the system (62) in an integral form which allows to pass to the limit under less restrictive assumptions, and also study an important case of a single charge in an external field.

2.2.
Single charge in an external EM field. In this section we describe assumptions under which the above derivation of Newton's law with the Lorentz forces for localized charges can be rigorously justified. In the case of many well separated charges every single charge senses all other charges through the sum of their EM field, and it is this feature of the EM interaction makes a problem of a single charge in external EM field of special importantce.
We consider a single charge described by (30) with N = 1 in an external EM field where ϕ ex (t, x), A ex (t, x) are two times continuously differentiable functions. Now index ℓ takes only one value, but we still keep it to be consistent with our previous notation. We assume that the first spatial derivatives of E ex and B ex are bounded uniformly in (t, x). We also assume that equation (62) holds in a weaker sense which we describe below.
First let us recast (62) in a different form. Multiplying (52) by x − r ℓ (t) we obtain Note that the right-hand side of the above equation coincides with m ℓ v ℓ (t, x). Substituting this expression for v ℓ into (62) (for a single particle) and we get where Observe then that according to the charge normalization Hence (67) can be written in the form Dynamics in an external electric field. Here we consider in detail a simpler case where external magnetic field is absent, A ex = 0. In this case and the field equations take the form As we already mentioned, in this section index ℓ takes only one value. First, let us recall some properties of equations (72), (73). Note that equation (72) does not involve the potential ϕ ℓ and can be solved independently from equation (73). If the nonlinearity g(ψ) = G ′ a |ψ| 2 ψ is a continuously differentiable function of ψ with uniformly bounded derivative, and ϕ ex is a smooth bounded function, then the local existence and uniqueness of a solution of (72) with a prescribed initial data in H 1 R 3 which belongs to and is bounded in H 1 R 3 is well known, see [15,Remark 4.4.8]. This solution is defined for all t > 0 if (110) holds for G ℓ since the energy is bounded from below on the invariant set ψ ∈ H 1 , ψ = 1 . The existence and uniqueness of solutions of the initial value problem for (72) with ϕ ex = 0 and the logarithmic nonlinearity as in (48) were proven in [16], see also [15]. Now we will show that, under a localization assumption, the centers r ℓ (t) of wave function ψ ℓ defined by (58) converge to trajectories governed Newton's law of motion for point charges. To eliminate the time derivatives in equation (67) we integrate it twice in time obtaining where We consider now equation (75) on a small time interval [0, T ] with solutions in the space of continuous functions C ([0, T ]). Since the first order spatial derivatives of E ex are bounded uniformly, the operator F satisfies a global Lipschitz estimate In addition to that we have Let us take T > 0 so small that and let z 0 be a solution of the equation . Such a solution exists and is unique by the contraction principle. Note then that the corresponding r (t) = z 0 (t) + r ℓ 0 +ṙ ℓ 0 t is a solution of (75) with ǫ 4 = 0 and consequently is a solution of the equation which coincides with Newton's equation (66) for a single point charge in an external electric field. To estimate a difference between the Newtonian trajectory r (t) and the charge position r ℓ (t) defined by (58) we subtract from (77) equation for z 0 and obtain z − z 0 = F (z) − F (z 0 ) + ǫ 4 . Using (78) and smallness of T we obtain the following inequality The above argument proves the following theorem.
x) are two times continuously differentiable functions, and that the first spatial derivatives of E ex are bounded uniformly in (t, x). Suppose also that: (i) there is a family of solutions ψ ℓ = ψ a of the NLS (72) with G = G a which depend on the size parameter a > 0; (ii) r ℓ (t) as defined by (58) is a continuous function of t ∈ [0, T ], where T is a sufficiently small fixed number which satisfies (79); (iii) (75) holds where ǫ 4 satisfies the following relation Note that the function E ex (t, x) − E ex t, r ℓ (t) = 0 at x = r ℓ and it is a bounded differentiable function, implying that condition (82) is a weaker version of for every a satisfy normalization condition (53), the convergence to δ x − r ℓ (t) means a weak form of localization around r ℓ (t). In Theorem 2.3 and Remark 3 we give an example where condition (82) holds.

2.2.2.
Dynamics in a general external electromagnetic field. Now let us consider the case when the external magnetic field described by its vector potential A ex does not vanish. For localized solutions of equation (30) where N = 1 and ℓ takes one value we want to derive the convergence of r ℓ (t) defined by (58) to a solution of Newton's equation of motion Considering equation (71) as a perturbation of equation (83) let us estimate the difference between their solutions. First we derive from equation (71) estimates for ∂ t r ℓ (t) 2 similar to usual energy estimates. Multiplying (71) by ∂ t r ℓ (t) we obtain . Integration of the above equation with the respect to t yields Consequently, where constant M 1 depends on constants M and T . Let r (t) be a solution of (83) with initial data In a similar way we obtain inequality (86) for (∂ t r (t)) 2 with a constant M 1 = M 10 which does not depend on ǫ 1 + ǫ 3 . To estimate the difference r ℓ − r we subtract (83) from (71), multiply the difference by ∂ t r ℓ − r and integrate the result with respect to t obtaining Therefore, using inequality (86) for (∂ t r (t)) 2 and ∂ t r ℓ (t) 2 , we obtain Using an elementary inequality where we obtain from (89) for t ∈ [0, T ] Rewriting inequality (92) in the form we obtain the following inequality with a constant C 2 which does not depend on ǫ 5 . Hence Z (t) → 0 as ǫ 5 → 0, and using inequalities (90) and (92) we conclude that r ℓ − r → 0 and ∂ t r ℓ − r 2 → 0 as ǫ 5 → 0. The above argument proves the following theorem.
Theorem 2.2. Suppose that A ex (t, x) and ϕ ex (t, x) are two times continuously differentiable functions, and that the first spatial and time derivatives of E ex and Suppose also that (i) there is a family of solutions ψ ℓ = ψ a of the NLS equations (30) with N = 1, G = G a which depend on the size parameter a > 0; (ii) r ℓ (t) and ∂ t r ℓ (t) as defined by (58) are continuous functions of t ∈ [0, T ]; (iii) equality (88) holds, and M defined by (85) is bounded uniformly in a, 0 < a ≤ a 0 ; (iv) the following limit relation holds where ǫ 1 , ǫ 2 , ǫ 3 are defined by (69) Remark 1. Note that the derivation of Newton's law of motion for charge centers in Theorems 2.1 and 2.2 does not depend explicitly on the particular nonlinearity G and the value of χ > 0. Analysis of the hydrogen atom model in Section 5.1 shows that χ = is a natural choice.
Remark 2. Theorems similar to Theorems 2.1 and 2.2 for a system of interacting charges can be analogously formulated and proven. Among additional conditions one has to assume the condition of non-collision, that is the trajectories of Newton's system (66) have to be separated on the time interval [0, T ]: Note that one cannot assume that the electrostatic potentials ϕ ℓ ′ have uniformly bounded gradients uniformly in a since they converge to the Coulomb potentials. But looking at (82) where E ex is replaced by −∇ϕ ex − ∇ϕ =ℓ we observe that uniform boundedness of ∇ϕ ℓ ′ is required only in a vicinity of r ℓ (t) with a fixed radius d/2, whereas outside the vicinity a natural assumption is convergence to zero of the integral of ψ ℓ 2 ∇ϕ =ℓ .

Exact wave-corpuscle solutions: accelerating solitons.
Here we consider the field equations (30), (31) for a single charge, omitting the index ℓ, and present a family of their exact solutions in the form of accelerating solitons (wave-corpuscles). Such solutions provide an example for which conditions of Theorem 2.1 hold. We assume a purely electric external EM field, i.e. when A ex = 0, E ex (t, x) = −∇ϕ ex (t, x). For the purely electric external field the field equations (30), (31) take the form (72). We define the wave-corpuscle ψ, ϕ by the following formula: In the above formulaψ is the form factor satisfying (35),φ is a radial function determined from (31). We refer to the function r (t) as wave-corpuscle center. Sinceψ is center-symmetric, this definition agrees with more general definition (58). Note that in a simpler case when the external fields ϕ ex and A ex vanish, a simpler solution of (72) and (83) is provided by (94) with r (t) = r 0 + vt with a constant velocity v. In this case the wave-corpuscle solution (94) of the field equations (72)-(73) can be obtained from the rest solutionψ,φ by certain Galilean-gauge transformations. Solutions of a similar form are known in the theory of nonlinear Schrödinger equations, see [36] and references therein. For the particular case of the logarithmic nonlinearity solutions of the form (94) were found in [12] in the form of accelerating gaussons. Theorem 2.3. Suppose that ϕ ex (t, x) is a continuous function which is linear with respect to x. Then ψ defined by (94) is an exact solution to (72), provided that r (t) is determined from the equation and Proof. If the potential ϕ ex (t, x) is linear in x then for any given trajectory r (t) we can write where Observe that the representation (94) implies and by Leibnitz formula we have Substituting the expression (94) into the field equations (72) we obtain the following equation for functions v, r, s p : Using the charge equilibrium equation (22) we eliminate the nonlinearity G and the Laplacian ∇ 2 in the above equation (100) and obtain the following equation equivalent to it: Now, we equate to zero the coefficients before ∇ψ andψ in that equation, resulting in two equations: where, in view of the representation (97), the second equation in (102) can be recast as Then we equate to zero the coefficient before (x − r) and the remaining coefficient and obtain the following pair of equations: Based on the first equation (102) and the equations (104) we conclude that the wave-corpuscle defined by the formula (94) and equations (95), (96) is indeed an exact solution to the field equations (72).
Remark 3. The above construction does not depend on the nonlinearity G ′ = G ′ a as long as (22) is satisfied. It is uniform with respect to a > 0, and the dependence on a in (94) Remark 4. The form (94) of exact solutions constructed in Theorem 2.3 is the same as the WKB ansatz in the quasi-classical approach. The trajectories of the charges centers coincide with trajectories that can be found by applying well-known quasiclassical asymptotics to solutions of (30) if one neglects the nonlinearity. Note though that there are two important effects of the nonlinearity not presented in the standard quasiclassical approach. First of all, due to the nonlinearity the charge preserves its shape in the course of evolution whereas in the linear model any wavepacket disperses over time. Second of all, the quasiclassical asymptotic expansions produce infinite asymptotic series which provide for a formal solution, whereas the properly introduced nonlinearity as in (22), (36) allows one to obtain an exact solution. For a treatment of mathematical aspects of the approach to nonlinear wave mechanics based on the WKB asymptotic expansions we refer the reader to [24] and references therein.
Similarly one can consider a single charge in an external EM field which can have nonzero magnetic component. In this case the wave-corpuscle is again defined by relations (94) but equations (95) are replaced by (83), see [5] for details.

Multiharmonic solutions for a system of many charges
In this section we consider a regime of close interaction which models a system of bound charges. This regime differs significantly from the regime of remote interaction considered in Section 2. More specifically, we seek special form solutions to the field equations (30), (31) for which (i) the potentials ϕ ℓ are time independent; (ii) every wave functions ψ ℓ depends on time harmonically through gauge factor e −iω ℓ t with possibly different values of frequencies ω ℓ for different ℓ. We refer to such solutions as multiharmonic. Note that multiharmonic solutions naturally arise in the analysis of nonlinear dispersive equations and systems, see [2].
Let us consider the field equations (30) and (31) without external fields, namely with ϕ =ℓ defined by (23) The total energy of the system is given by the formula where ϕ =ℓ is defined in (106), ϕ ℓ ′ are given by (32). The energy E is a conserved quantity for the system (105), (106), it can be derived in a standard way from the Lagrangian (26) taking in account equations (105), (106). Collecting terms in the expression for E which explicitly involve ψ ℓ with a given ℓ we introduce ℓ−th charge energy in the system's field by the formula Here and below integrals with respect to spatial variables are over R 3 . We call so defined E 0ℓ the ℓ -th particle energy in system field. Note that for a single particle E 0ℓ = E ℓ = E, in the general case N ≥ 2 the total energy E does not equal the sum of E 0ℓ , see for example (144) for the case of two interacting charges; in fact the difference between the sum of E 0ℓ and the total energy E coincides with the total energy of EM fields. In accordance with (53) we assume that ψ ℓ ∈ Ξ where Theorem 3.1. Let G ℓ (s) be a function of class C 1 which is subcritical, namely there exists δ > 0 and C such that Then the functional E is bounded from below on the set Ξ N and boundedness of E implies boundedness of ψ ℓ H 1 , ℓ = 1, ..., N . Boundedness of E on Ξ N is equivalent to boundedness of all E 0ℓ , ℓ = 1, ..., N.
Proof. Note that according to (32) Splitting the domain of integration into |x − y| ≥ θ and |x − y| < θ with arbitrary small θ and using the Sobolev imbedding theorem we obtain the following inequality where ǫ can be chosen arbitrary small. From this inequality and (110) we obtain the statement of the theorem.
Let us consider now multiharmonic solutions to the system (105), (106), that is solutions of the form Note that we use lower index for ψ ℓ (x) to differ it from ψ ℓ (t, x). Then the functions ψ ℓ (x) satisfy the following nonlinear eigenvalue problem where ϕ =ℓ is defined by (23) We want the above relation (115) to hold for any regular multiharmonic solution of a system of the form (113), (106). Let us consider first the frequencies ω ℓ which can be determined from the equations as follows. Multiplying the ℓ-th equation (113) by ψ * ℓ , integrating the result with respect to the space variable and using the charge normalization condition we obtain Comparing the above with (108) we see that For two solutions{ψ ℓ } N ℓ=1 , ψ ′ ℓ N ℓ=1 we obtain the following relations linking their frequencies and energies to the nonlinearities: To find a sufficient condition for (115) to hold, we take into account the charge normalization condition ψ ℓ = 1 and observe that it is sufficient to impose the following condition for every |ψ ℓ | 2 with ψ ℓ = 1, where the constant may depend on ℓ. This condition is fulfilled if the following differential equation with independent variable s = |ψ| 2 holds: with K G that depends only on ℓ. Then (117) and the normalization condition imply and hence (115) is fulfilled. Solving the differential equation (120) we obtain the following explicit formula for the logarithmic nonlinearity G ℓ (s) = K G s ln s + Cs.
Comparing with (50) and (49) we find that if K G < 0 this is exactly the nonlinearity (50) which corresponds to the Gaussian factor. According to (49) the above formula takes the form where a ℓ is the size parameter for ℓ-th charge. If K G = 0 then G ℓ |ψ ℓ | 2 is quadratic and equation (113) turns into the linear Schrödinger equation for which fulfillment of the Planck-Einstein relation is a well-known fundamental property. It is absolutely remarkable that the logarithmic nonlinearity which is singled out by the fulfillment of the Planck-Einstein relation has a second crucial property, namely it allows a localized soliton solution, namely the Gaussian one in (47). Note that above argument can be repeated literally if equation (113) involves external timeindependent electric field as in (30), in this case ϕ =ℓ in the definition (108) of E 0ℓ has to be replaced by ϕ =ℓ + ϕ ex .
Let us study in more detail the situation with the logarithmic nonlinearity singled out by the Planck-Einstein relation. Formula (121) for the Gaussian wave function takes the form . Then for any two solutions ψ σ ℓ = ψ ℓ and ψ σ1 ℓ = ψ ′ ℓ with σ, σ 1 ∈ Σ the Planck-Einstein relation (115) is fulfilled. Proof. According to Theorem 3.2 the boundedness of E σ 0ℓ for every given σ implies that functions ψ σ ℓ are bounded in H 1 and in X. Let us consider the left-hand side of (113) and denote it by F (ψ ℓ ). Similarly to (111) we obtain that where C does not depend on f . Hence multiplication by ϕ ℓ ′ is a bounded operator from H 1 to H −1 . Using this fact to obtain continuous dependence of the term ϕ =ℓ ψ ℓ on ψ ℓ and using Lemma 9.3.3 in [15] for remaining terms in F (ψ ℓ ), we conclude that the functional F (ψ ℓ ) depends continuously in space W * defined in (138) on ψ ℓ ∈ H 1 ∩ X and it is bounded in W * . Since (117) holds for smooth rapidly decaying functions ψ ℓ , we can use in a standard way that smooth functions with compact support are dense in H 1 ∩ X and obtain (117) written in the form for functions ψ σ ℓ from H 1 ∩ X. The right-hand side does not depend on σ thanks to the normalization condition (53), and we obtain (115).
We would like to stress remarkable universality of the Planck-Einstein relation (115) proven above for multiharmonic solutions and the fact that it holds for every individual particle in a system of many interacting particles. Namely, the Planck-Einstein relation holds with the same coefficient χ for an arbitrary pair of solutions of (113), (106) and for an arbitrary component of those solutions, and the validity of this relation does not depend on the number N which equals the number of interacting particles. Observe also that the derivation of (121) and (122) is based only on properties of the nonlinear terms as in (118) and can be extended to more general systems than (105), (106). In particular, the systems may involve linear terms with potentials which explicitly depend on x as in (172).
The connection between the logarithmic nonlinearity and the Planck-Einstein frequency-energy relation was discovered in a different setting by Bialynicki-Birula and Mycielski in [12]. Note though that in [12] a system of N particles is described as in the quantum mechanics by a single wave function ψ over 3N -dimensional configuration space, whereas in our model every of N particles has its own wave function ψ ℓ depending on 3 spatial variables. This is why our approach naturally leads to the study of multi-harmonic solutions for interacting particles for which every individual wave function ψ ℓ can be associated with its individual frequency ω ℓ . Another signficant difference between our approach and the nonlinear mechanics in [12] is the way the nonlinearity enters the Lagrangian and the field equations for N ≥ 2 particles. Namely, in our approach every individual particle has its own nonlinearity whereas in [12] there is a single nonlinearity for the entire system of N particles.

3.2.
Gaussian shape as a global minimum of energy. As one can see from the previous section the logarithmic nonlinearity defined by (48) is singled out by the requirement that the Planck-Einstein frequency-energy relation holds exactly. Let us look further into the properties of the logarithmic nonlinearity.
We already observed that the equilibrium equation with the logarithmic nonlinearity admits the Gaussian function as an exact solution. This solution is a critical point of the energy, and we show below that the Gaussian function provides a unique (modulo translations) minimum of the energy functional on functions subjected to the charge normalization condition. This fact was already found in [12].
Let us consider the energy E (ψ) = E 0 (ψ) restricted to one particle with ϕ ex = 0 and G ℓ defined for a = 1 by (49), namely Of course, the case a > 0 can be considered similarly. We consider this functional on a set Ξ defined by (109). In the following theorem we show based on the logarithmic Sobolev inequality that the Gaussian function ψ g = C g e −|x| 2 /2 provides the global minimum of the energy, and that all the global minima belong to the set Ω = ψ : ψ = e iθ ψ g (· − r) , r ∈ R 3 , θ ∈ R (124) obtained from ψ g by spatial translations and gauge factor multiplication. Another proof of this statement is given in [14].
and the equality holds only for ψ ∈ Ω.
Proof. For simplicity we set here χ 2 m ℓ = 1. The value of the energy E (ψ) on ψ g can be found explicitly: Note that the Euler equation for constrained critical points of E (ψ) has the form λψ + ∇ 2 ψ = − ln |ψ| 2 + ln 1 and ψ g satisfies this equation. Now we show that ψ g provides the global minimum of E (ψ) defined by (123) for real-valued ψ. We will use the following Euclidean logarithmic Sobolev inequality for functions from H 1 R d where equality holds only for Gaussian functions (see [11]). Therefore, To find R d |∇ψ| 2 dx we consider the Euler equation (127) for the minimizer ψ. From Pohozhaev identity applied to (127) we obtain: where G |ψ| 2 = − |ψ| 2 ln |ψ| 2 + |ψ| 2 ln 1 π 3/2 − 2 .

Multiplying equation (127) by ψ and integrating the result we get
Adding (130) and (131)  Substitution of the above expression into inequality (129) implies E (ψ) ≥ 1 2 where the equality holds only for Gaussian functions. Since ψ is a minimizer, the inequality holds for any real ψ ∈ Ξ. Let us consider now a complex minimizer ψ = u + iv such that E (ψ) ≤ 1 2 . Notice that the following elementary inequality holds if u 2 + v 2 = 0 and the equality is possible only if v∇u = u∇v. The above inequality implies in particular that u 2 + v 2 1/2 ∈ H 1 . Applying once more the Euclidean logarithmic Sobolev inequality for real functions we get For such a minimizer u 2 + v 2 1/2 must coincide with a Gaussian function and equality v∇u = u∇v must hold a.e.. Since for non-zero u and v, we conclude that ln u − ln v is a constant and the complex minimizer ψ is obtained from a Gaussian function via multiplication by a gauge factor e iθ .
The following theorem directly follows from results of Cazenave and Lions [17].
and M δ be a sublevel set Then for any ǫ > 0 there exists δ > 0 such that M δ ⊂ Ω ǫ .
Proof. Assume the contrary, namely that there exists ǫ > 0 and a sequence ψ j ∈ Ξ such that E ψ j → E ψ g and inf v∈Ω ψ j − v H 1 ≥ ǫ. According to Theorem II.1 and Remark II.3 in [17] a subsequence ψ j (· − r j ) is relatively compact in H 1 . Hence there is a subsequence that converges in H 1 to a global minimizer which belongs to Ω. This contradicts the assumption ǫ > 0.
Let us consider now one-particle equation (105) with logarithmic nonlinearity in the absence of external fields (a free particle) The existence and the uniqueness of solutions to the initial value problem for this equation is proven in [16]. Since the logarithm has a singularity at ψ = 0, to describe classes of functions for which the problem is well-posed we need to introduce special spaces. Following Cazenave [15] we introduce functions B (|ψ|) = |ψ| 2 ln |ψ| 2 + A (|ψ|) .
(136) So defined A is a convex C 1 function which is C 2 for |ψ| = 0. The function B (|ψ|) vanishes for small |ψ| and ψB (|ψ|) / |ψ| 2 is of class C 1 and satisfies uniform Lipschitz condition. Let A * be the convex conjugate function of A, it is also convex of class C 1 and positive except at zero. We introduce Banach spaces Note that M δ is bounded in W . By Lemma 9.3.2 from [15] the function ψ → ψ ln |ψ| 2 is continuous from X to X ′ and is bounded on bounded sets. According to Cazenave [14], [15], if ψ 0 ∈ W there exists unique solution ψ ∈ C (R, W ) ∩ C 1 (R, W * ) of (134) with initial data ψ (0) = ψ 0 . Properties of solutions of (134) are described in [14], [15], in particular Ξ ∩ W and M δ are invariant sets. According to [17] the solution ψ g is orbitally stable, namely invariance of M δ and Theorem 3.5 imply that if ψ (0, ·) ∈ M δ then ψ (0, t) ∈ Ω ǫ for t ≥ 0. Consequently, small initial perturbations of the Gaussian shape do not cause its large perturbations in the course of time evolution given by (134) but may cause considerable spatial shifts of the entire wave function. Such a time evolution of charges in electric field is described in Theorem 2.1. Note that if the potential ϕ ex (t, x) is linear as in Theorem 2.3, one can rewrite equation (72) in a moving frame with origin r (t) and corresponding phase correction for ψ as in (94) (see [5] for details) and obtain an equivalent equation of the form (72) with ϕ ex (t, x) = 0, namely (134). Wavecorpuscle solution (94) after the change of variables turns into ψ g . Therefore, using Theorem 3.5 we conclude that the wave-corpuscle solutions of (72) constructed in Theorem 2.3 are orbitally stable.

Two particle hydrogen-like system
In this section we consider a particular case of multiharmonic solutions of the system (113) with N = 2 which models a bound proton-electron system. For this system we discuss general properties and provide a heuristic analysis of the coupling between proton and electron. Using this analysis as a motivation, we then write equations and energy functional for one charge (electron) in the Coulomb field of proton. We study this problem in more detail along the lines of [9], [10].
4.1. Electron-proton system. To model states of a bound proton-electron system, the hydrogen atom, we consider the multiharmonic solutions of (30), (31) described by the system (113), (106) with N = 2, where indices ℓ take two values ℓ = 1 and ℓ = 2. The charges have opposite values q 1 = −q 2 = −q. For brevity we introduce notation The quantity a 1 turns into the Bohr radius if χ coincides with Planck constant , and m 1 , q are the electron mass and charge respectively. Using the above notation we rewrite the two-particle system (113), (106) in the form of the following nonlinear eigenvalue problem χ q 2 ω 1 ψ 1 + (143) The functions ψ 1 and ψ 2 are respectively the wave functions for the electron and the proton. As always, we look for solutions ψ ℓ ∈ Ξ with Ξ defined by (109). The nonlinearities G ′ 1 , G ′ 2 are assumed to be logarithmic as defined by (50) where the size parameter a = a ℓ is different for electron and proton. One can similarly consider other nonlinearities but in this paper we stay with the logarithmic ones.
In accordance with (107) we introduce the energy functional by the following formula where −∇ 2 −1 |ψ ℓ | 2 is defined by (32). Equations (141), (142), (143) can be derived from the Euler equations for a critical point of E (ψ 1 , ψ 2 ) with the constraint the frequencies ω 1 , ω 2 being proportional to the Lagrange multipliers. Using (32) one can see that hence the coupling term in (144) is symmetric with respect to ψ 1 , ψ 2 . Note also that the first term in formula (144) coincides with the energy E 01 of the first charge in the system field as given by (108) with ℓ = 1, and the second term with the expression for the energy E 2 (ψ 2 ) of a free second particle given by (123). Obviously, E (ψ 1 , ψ 2 ) also can be written as a sum of energy E 02 of the second charge in the system field plus the free energy of the first charge. According to Theorem 3.2 the energy E (ψ 1 , ψ 2 ) is bounded from below.
Remark 5. Note that system (141)-(143) is similar to the Hartree equations studied in [27], though it differs from it because of the presence of the nonlinearities G ℓ .

4.2.
Reduction to one charge in the Coulomb field. In this subsection we give a motivation for replacing the problem of determination of critical values of the energy functional E given by (144) for a system of two charges by a simpler problem for a single charge similarly to the Born-Oppenheimer approximation in quantum mechanics. To this end we use two changes of variables in two equations, x = a 1 y 1 in (141) and x = a 2 y 2 in (142), where a ℓ are defined by (140), namely and rescale the fields as follows: Hence (141), (142) takes the form of the following nonlinear hydrogen system (149) Here we use the same letter G to determine the function of rescaled variables. Note that the electron/proton mass ratio m1 m2 ≃ 1 1837 is small, therefore the parameter is small too. We rewrite (147), (148) as follows: We would like to show that if we are interested in electron frequencies ω 1 which correspond to lower energy levels of E (ψ 1 , ψ 2 ) then, in the case b → 0, it is reasonable to replace system (151), (152) by one equation with the Coulomb potential.
The discrete spectrum of the classical linear non-relativistic Schrödinger operator for hydrogen atom with the Coulomb potential can be recovered (if multiplicities are not taken into account) from radial eigenfunctions. Therefore here we restrict ourselves to the radial solutions of (151), (152) and restrict (144) to radial functions. Consider a radial solution (Ψ 1 , Ψ 2 ) of (151), (152). Equation for the electron frequency ω 1 given by (116) with ℓ = 1 takes the form (153) with the proton potential φ 2 . Then we would like to replace the potential φ 2 in (153) by the Coulomb potential 1 |y| as an approximation. To show that for small b one may expect the resulting perturbation of eigenvalues to be small, we use the following representation for the radial potential φ (r) = 4π −∇ 2 −1 |Ψ| 2 : The term which we expect to be small is Changing the order of integration and using substitution r/b = r 2 we obtain Hence The right-hand side of (156) involves the variance of Ψ 2 and Ψ 1 L ∞ . To show that |D prot | is small it is sufficient to show that the variance of Ψ 2 and Ψ 1 L ∞ are bounded uniformly for small b. Note that Ψ 1 , Ψ 2 in (156) are not arbitrary radial functions, but special solutions of (151), (152), (149) which correspond to lower critical energy levels of the energy functional E = E 01 + E 02 . Dependence of the coupling between (151) and (152) (151) and (152) respectively, according to (154) these potentials are bounded by 1 |y| uniformly for b > 0. Using the uniform boundedness of the potentials and estimate (182) we conclude that energy and H 1 -norms of Ψ 1 and Ψ 2 are bounded uniformly in b. As in the proof of Lemma 5.3 below, we conclude that radial Ψ 1 and Ψ 2 are bounded in L 5 , they also satisfy (210). From boundedness in L 5 we derive that 1 |y| Ψ 1 is bounded in L 7/4 and so does G ′ 1 |Ψ 1 | 2 Ψ 1 . Note also that since the energy is bounded from below, the negative values of ω = χa1 q 2 ω 1 in (151) which correspond to lower energy levels are bounded. Hence, the term (151) is bounded in L 7/4 , therefore the solution Ψ 1 of the elliptic equation (151) belongs to Sobolev space W 7/4,2 and by Sobolev imbedding is bounded in L ∞ . As for Ψ 2 , from (210) using the comparison argument as in the proof of Lemma 5.8 we derive exponential decay of solutions of the equation (152) for large r uniformly in b and conclude that the variance of Ψ 2 is uniformly bounded for small b. Hence, for such solutions |D prot | is small for small b. The smallness of b in (156) (recall that b ≃ 1 1837 for proton-electron system) gives us a motivation to set b = 0 in (152) and replace 1 b φ 2 1 b y by 1 |y| in (151). So, we replace the problem of finding frequencies ω 1 by formula (153) based on critical points of E (ψ 1 , ψ 2 ) which correspond to lower energy levels of E (ψ 1 , ψ 2 ) by the problem of finding critical points, lower critical levels and corresponding frequencies ω 1 (which for the logarithmic nonlinearity are expressed in terms of the critical levels by (122)) for the following energy functional with the Coulomb potential: This functional involves only Ψ 1 . In quantum mechanics a somewhat similar reduction (in a completely different setting) from proton-electron system to a single equation with the Coulomb potential is made via Born-Oppenheimer approximation. Additional motivation for the reduction is given in Remark 10.

Variational problem for a charge in the Coulomb field
In this section we consider in detail radial critical points of the functional E Cb (Ψ 1 ) defined by (157) and establish its basic properties. In this and following section we look for real solutions Ψ = Ψ * and all function spaces involve only real functions. In what follows we explicitly take into account the dependence of the nonlinearity on the size parameter a = a 1 (or on related parameter κ), namely where, for consistency with the notation (33) where a = a 1 /κ, we use notation /κ to denote the dependence of rescaled G on κ. Note that the dependence on a 1 and a = a 1 is factored out after the rescaling, and that G /κ in the above formula depends only on their ratio κ. We consider the case of small κ, namely when the electron size parameter a is much larger than the Bohr radius a 1 . To treat technical difficulties related to the singularity of logarithm at zero, together with the logarithmic nonlinearity given by (50) we use a regularized logarithmic nonlinearity defined as follows: where ξ ≤ 0, ln +,ξ (s) = max (ln (s) , ξ) for s > 0, ln +,ξ (0) = ξ.
The nonlinear eigenvalue problem with the Coulomb potential has the form similar to (151): where the dimensionless spectral parameter ω = ω κ,ξ is given by the formula Note that equation (172) can be obtained as the Euler equation by variation of the energy functional (157) which we write in the form where E 0 is the quadratic energy functional and G κ,ξ is the nonlinear functional As always we assume the charge normalization constraint and denote similarly to (109) the set of radial functions (which depend only on |x|) by Obviously, for Ψ ∈ Ξ ∩ X where X is defined by (137) (179) The spectral parameter ω is the Lagrange multiplier and it relates with corresponding critical energy levels E κ,ξ n of E κ,ξ on Ξ ∩ X by formula (122) which takes the form E κ,ξ n = ω κ,ξ,n + Let us introduce for κ > 0, ξ ≤ 0, Ψ ∈ C the function General properties of functional E κ,ξ (Ψ) with regularized logarithm are described in the following statements.
Proof. The estimate (183) together with the normalization condition imply the uniform boundedness of Ψ j in the Sobolev space H 1 R 3 . According to Theorem A.I' of [9], the space of radial functions H 1 rad R 3 is compactly imbedded into L p R 3 if 2 < p < 6. Hence, a subsequence Ψ j (x) → Ψ ∞ (x) strongly in L p R 3 and almost everywhere. Statement (ii) follows from Radial Lemma A.II in [9].
The functional E κ,ξ (Ψ) with regularized logarithm satisfies the Palais-Smale condition similar to condition (P-S + ) in [10]: Let a sequence Ψ j ∈ Ξ rad have the following properties: and strongly in H −1 R 3 . Then the sequence contains a subsequence which converges strongly in H 1 rad R 3 . Proof. The proof is similar to treatment of radial solutions in [9], [27]. We only sketch main steps. Since (182) holds, (183), (165) and (128) imply boundedness of Ψ j in H 1 R 3 . Using Lemma 5.3 we conclude that we can choose a subsequence such that Ψ j → Ψ ∞ in L 5 R 3 and L 8/3 R 3 , it converges weakly in H −1 R 3 and almost everywhere in R 3 . Multiplying (184) by Ψ j in L 2 R 3 = H with inner product ·, · we obtain We rewrite (185) in the form Using (164) we obtain for Ψ ∈ Ξ rad Since F Ψ j , Ψ j → 0, (183) and (186) imply that we can choose a subsequence Ψ j with ω j → ω ∞ where ω ∞ ≤ −β. Multiplying (184) by a smooth test function with compact support and passing to the limit we obtain in standard way that Ψ ∞ ∈ H 1 R 3 satisfies (172) with ω = ω ∞ .

5.1.
Nonlinear eigenvalues for a charge in the Coulomb field. In this section we prove that if κ is small enough, the nonlinear eigenvalue problem (172) with logarithmic nonlinearity has solutions with eigenvalues which are close to the eigenvalues of the corresponding linear problem. We look for real-valued solutions.
As a first step we prove existence of such solutions of the problem with regularized logarithm ln +,ξ . As the second step we pass to the limit as ξ → −∞ and obtain solutions of the eigenvalue problem with the original logarithmic nonlinearity. Consider the linear Schrödinger operator O with the Coulomb potential which corresponds to E 0 (Ψ), that is Note that O has the following well known negative eigenvalues which coincide with the negative eigenvalues of the operator O rad obtained by restriction of O to radial functions (without counting their multiplicity) and are equal to the negative critical energy levels of E 0 (Ψ). Every eigenvalue of the radial problem O rad Ψ 0 n = ω 0,n Ψ 0 n is simple and corresponding eigenfunctions Ψ 0 n (r) given by well-known explicit formulas decay exponentially as r → ∞, for example the linear ground state is given by the formula Let us denote by L 0 n an eigenspace of O rad which corresponds to the eigenvalue E 0 n and consists of radial functions, by L − n ⊂ H 1 rad a finite-dimensional invariant subspace of O rad which corresponds to eigenvalues λ ≤ E 0 n and by L + n ⊂ H 1 rad R 3 invariant infinite-dimensional subspace in H 1 rad R 3 of O rad which is orthogonal in L 2 to L − n . We have an orthogonal decomposition H 1 rad R 3 = L − n + L + n = L − n−1 + L 0 n + L + n . We look for eigenvalues ω = ω κ,ξ,n < 0 of the nonlinear eigenvalue problem (172) which can be considered as a perturbation of the eigenvalues E 0 n . The ground state of E κ,−∞ (Ψ) and the corresponding energy level ω κ,−∞,1 was found in [12], it is given by the formula where Hence the nonlinear correction for the ground level is of order κ 2 .
Since the Palais-Smale condition is satisfied according to Theorem 5.4, the existence of an infinite sequence of discrete negative energy values of E κ,ξ (Ψ) can be proven exactly along the lines of [9]. But we want to show, in addition, that the eigenvalues of the nonlinear problem with small κ are close to the eigenvalues of the linear problem (this is not surprising, since all the eigenvalues of the linear problem for radial functions are simple). Therefore we modify the construction of critical points of the functional E κ,ξ on Ξ.
We will use the following deformation lemma which is a minor modification and reformulation in our notation of Lemma 5 from [9] (with essentially the same proof).
Lemma 5.5. Suppose that the functional E κ,ξ ∈ C 1 (Ξ rad , R) satisfies Palais-Smale condition as in Theorem 5.4. Suppose also that E ′′ < 0 and that the segment E ′ ≤ E ≤ E ′′ does not contain critical values of E κ,ξ . Letǭ > 0 be arbitrary small number. Then there exists 0 < ǫ <ǭ and a mapping η ∈ C (Ξ rad , Ξ rad ) such that To prove existence of critical points of E κ,ξ (Ψ) on the infinite-dimensional sphere Ξ rad with critical values E κ,ξ n near E 0 n , n = 2, 3... when κ is small, we slightly modify the construction from [9]. As a first step we define finite-dimensional domain B = B n ⊂ L − n which lies in Ξ rad near eigenfunction Ψ 0 n of O rad as follows: where R 0 < 1/2 is a small fixed number, for simplicity we fix R 0 = 1/3. The domain B n is diffeomorphic to a n − 1 dimensional ball. The boundary of B n is the n − 2 dimensional sphere We consider continuous functions Γ defined on B n with values in Ξ. We consider the class {Γ} n of continuous mappings Γ : B n → Ξ rad such that Γ restricted to ∂B n is identical, The n-th min-max energy level E κ,ξ,n is defined as follows: Lemma 5.6. Let ξ ≤ 0, let n = 2, 3, ...., 0 < κ ≤ 1. Then E κ,ξ,n defined by (196) lies in the interval where C 1n and C 3n depend only on n.
Proof. Since we can take a particular function Γ (Ψ) = Γ 0 (Ψ) = Ψ 0 n + v for Ψ ∈ B n we have Note that B n lies in a ball in a finite-dimensional subspace L − n in H 1 , the basis of L − n consists of exponentially decaying functions with bounded first derivatives, hence B n is compact in H 1 R 3 and in L p R 3 for every p ≥ 1. Therefore Hence, using (165) and power estimates for Ψ 2 ln Ψ 2 we conclude that and we obtain max where C 1n depends only on n. Using Γ 0 (Ψ) and the estimate (198) of E κ,ξ (Ψ) we observe that where C 2n = C 1n κ 2 + E 0 n . Inequality (199) provides E ′′ in (197). Now we estimate E κ,ξ,n from below. Mappings Γ ∈ {Γ} n which satisfy (195) have the following topological property: the set Γ (B n ) must intersect the subspace L + n−1 . Indeed, assume that the contrary were true and for some Γ there were not intersection of Γ (Ψ) with L + n−1 , that is Π − n−1 Γ (Ψ) = 0 for all Ψ ∈ B n where Π − n−1 is the orthoprojection onto the finite-dimensional subspace L − n−1 . Then we could consider the following function on the ball v ≤ R 0 in L − n−1 : It is a well-known fact from the algebraic topology that such a retraction cannot exist. Therefore for every Γ we have Π − n−1 Γ (Ψ 0 ) = 0 for someΨ 0 ∈ B. Therefore there exists Ψ 0 ∈ B n such that To obtain estimate of E κ,ξ,n from below note that (201) implies Inequality E κ,ξ (Γ (B n )) ≤ C 2n by Lemma 5.2 implies the boundedness Γ (B n ) in H 1 . Using (iv) in Lemma 5.2 we conclude that Combining (202), (198) and (199) we obtain (197). Now we prove the existence of solutions of eigenvalue problem with regularized logarithm.
Now we show that solutions of (172) decay exponentially as r = |x| → ∞.
Proof. Since Ψ n H 1 are bounded uniformly in ξ ≤ 0 and κ ≤ 1 we can use the following result of Strauss (see Lemma A.II in [9]) where C ′ 0 andR are absolute constants,R ≥ 1. We apply this inequality to radial solutions Ψ = Ψ n of (172) with ω = ω κ,ξ,n < 0 and obtain where, according to Lemma 5.2, the constant C 0 is bounded uniformly in κ ≤ 1, ξ ≤ 0 and n. We write (172) in the form where According to (211) using monotonicity of the function ln +,ξ (s) with respect to s we obtain for |x| ≥ R ≥R ≥ 1, Therefore, for any 0 < κ ≤ κ 0 ,ξ ≤ 0 we can find R 1 ≥R such that ln κ −3 R −2 1 C 2 0 /C 2 g ≤ξ, and if ξ is so large that (208) holds, we obtain for ξ =ξ According to definition (160) V κ,ξ (x) is a monotonely increasing function of ξ, therefore this inequality holds for all ξ ≤ξ. Note that Ψ (x) cannot change sign for |x| > R 1 . Indeed, if Ψ (x) = 0 at |x| = R ′ 1 > R 1 we could multiply (212) by Ψ and integrate over {|x| > R ′ 1 } which would yield and Ψ (x) = 0 for all |x| > R ′ 1 which is impossible for a radial solution of (172) with Ψ = 1. So we assume Ψ (x) > 0 for |x| > R 1 (the case of negative Ψ is similar). We compare the solution Ψ of this equation in the domain |x| > R with the function which is a solution of the equation Subtracting half this equation from (212) we obtain We choose p 2 = κ 2 9 ξ and (213) implies We choose C Y > 0 so that C Y e −pR1 > C 0 , where C 0 is the same as in (211). We assert that This is true at |x| = R 1 according to the choice of C Y , the limit of Ψ − Y as |x| → ∞ is zero. Assume now that (216) does not hold for some |x| > R 1 . Then Ψ (x) − Y (x) must have a point of local positive maximum at some |x| = R ′ 1 > R 1 , at this point ∇ 2 (Ψ − Y ) ≤ 0, by (215) V κ,ξ + p 2 Y < 0 and V κ,ξ (Ψ − Y ) ≤ 0 by (213). This inequalities contradict (214), hence (216) holds. Inequality (209) follows from (216). Now we prove the main result of this section for the nonlinear eigenvalue problem with the original logarithmic inequality which corresponds to ξ = −∞.
Since the eigenfunctionsΨ j of O are bounded, continuous and decay exponentially, we see that operators (ω j − ω) Ψ,Ψ j Ψ j are bounded from L p R 3 to L p ′ R 3 for any given p, p ′ , ∞ ≥ p ≥ 1, ∞ ≥ p ′ ≥ 1. Solutions Ψ ∈ Ξ ∩ X ′ are bounded in H 1 , and by Sobolev imbedding in L 6 R 3 .
Remark 6. Theorem 5.9 states that for every energy level of the linear Schrödinger hydrogen operator there exists an energy level of nonlinear hydrogen equation with the Coulomb potential with an estimate of difference of order ∼ κ 2 ln κ. Based on (173) we see that if χ = where is the Planck constant and κ = 0 then ω 1 = ω 1 n coincide with standard expression for hydrogen frequencies, and a natural choice for parameter χ is χ = .
Remark 7. A theorem similar to Theorem 5.9 (without superexponential decay of solutions) can be proven not only for the logarithmic nonlinearity, but for more general subcritical nonlinearities which involve a small parameter κ. Moreover, if corresponding functions g κ are of class C 1 and their derivatives satisfy natural growth conditions, then the standard bifurcation analysis shows that every E 0 n generates a branch E κ,n of solutions Ψ n ∈ Ξ rad of (172) for small κ.
Remark 8. In Theorem 2.1 we assumed in the macroscopic case electron size a to be very small and in this section we assume κ = a 1 /a to be very small, which seem to contradict one another. But a closer look shows that in fact the both cases are consistent one with another, if we assume the possibility a 1 → 0. To show that this assumption is quite reasonable, we look at the values of physical constants. Note that a typical macroscopic scale is a macr ∼ 10 −3 m and the Bohr radius a 1 ∼ 5.3 × 10 −11 m. Note also that the error in (82) is of order a 2 /a 2 macr ≪ 1 and here κ 2 = a 2 1 /a 2 ≪ 1 which is possible if a 2 1 ≪ a 2 ≪ a 2 macr . Since a 2 1 /a 2 macr ∼ 2.8 × 10 −15 we consider a 2 1 /a 2 macr → 0 as a reasonable approximation. Remark 9. In the case when the energy functional with the Coulomb potential E Cb (Ψ) is considered not only on radial functions, the eigenvalues of the linear problem are not simple and have very high degeneracy. Note that if we start from the relativistic equation (13) and look for ω-harmonic solutions we arrive to a nonlinear perturbation of so-called relativistic Schrödinger equation, see [34]. The spectrum of this linear problem has a fine structure and the degeneracy with respect to the total angular momentum disappears. If we fix the value of z-component L z of the angular momentum L, for example by looking for solutions of the form Ψ 1 =e imzθ Ψ 10 (r, z) with a fixed integer m z , one can eliminate the remaining degeneracy in the linear problem. This provides a ground to expect that the lower energy levels of the nonlinear relativistic Schrödinger equation can be approximated by the lower energy levels of the linear relativistic Schrödinger equation if κ is small.
Remark 10. Note that the decoupled equation (152) with b = 0 for proton has a unique ground state solution (gausson) which corresponds to the minimum E p 1 of energy (see Theorem 3.4 ) and has a sequence of solutions with higher critical energy levels E p j , j ≥ 2 (see [14]). Since the ground state is unique, E p j > E p 1 if j > 1.