text:
# Generic singularities of implicit systems of first order differential equations on the plane

## 1 Introduction.

An implicit system of first order odinary differential equations on a smooth $`n`$-dimensional manifold is defined by a zero level, which is called system surface, of a smooth map $`F`$ from the tangent bundle of this manifold to the $`n`$-dimensional Cartesian space. In local coordinates $`x=(x_1,\mathrm{},x_n)`$ near a point of the manifold a system can be written in the standard form $`F(x,\dot{x})=0.`$

An implicit system is with locally bounded derivatives if the restriction of the bundle projection to the system surface is a proper map. This restriction is called the system folding. Only systems with locally bounded derivatives are considered below. We identify the space of systems with locally bounded derivatives with the space of respective maps $`F.`$ Then a generic system is a system from some open everywhere dense subset in this space, for the fine topology of Whitney.

A solution curve of an implicit system is defined as a differentiable map $`x:tx(t)`$ from an interval of the real line to the base manifold such that the image of its natural lifting $`(x(t),\dot{x}(t))`$ to the tangent bundle belongs to the system surface. A phase curve is the image of such a map $`x(t)`$ and a trajectory is the image of its lifting.

In this paper we study the point singularities of the family of phase curves provided by a germ of a system surface and present the complete list of generic singularities on the plane up to smooth orbital equivalence.

For a generic system its surface is a closed smooth $`n`$-dimensional submanifold in the tangent bundle space, owing to Thom transversality theorem. Then the system folding is a smooth map between $`n`$-dimensional manifolds. Therefore the folding of a generic system can have all generic singularities like a map between $`n`$-dimensional manifolds. In fact the kernel of the bundle projection is also $`n`$-dimensional and, due to Goryunov’s theorem, such dimension of the kernel permits all generic singularities for maps between $`n`$-dimensional manifolds.

A generic system near a regular point of its folding can be resolved with respect to derivatives. In this case, near such a point, the singularity theory of family of solutions of an implicit system is reduced to the well known theory of family of the phase curves for the generic smooth vector fields on $`n`$-dimensional manifolds .

For a generic system the velocity does not vanish at any singular point of the system folding. Consequently, near such a point, there is well-defined the system 1-folding which is the restriction of the projectivization of the tangent bundle to the system surface. Again Goryunov’s theorem implies that the 1-folding of a generic system can have all singularities like a generic map from an $`n`$-dimensional manifold to an $`(2n1)`$-dimensional one.

In particular, for the $`2`$-dimensional case the 1-folding of a generic system can have regular points and singular points providing Whitney umbrella singularities. That implies, for a generic implicit system, the classification of point singularities of families of phase curves (see Theorem 2.1, 2.4). Besides well known singularities of a generic vector fields on the plane and ones described by generic first order implicit differential equations there is only one singularity provided by the implicit differential equation on the Whitney umbrella generically embedded to the space of directions on the plane (Figures 1, 2, and 3). Up to smooth orbital equivalence the respective family of phase curve is the family of solutions of the implicit system

$$\dot{x}=\pm 1,(\dot{y})^2=x(xy)^2$$

near the origin.

The family of solutions of the implicit equation $`(dy/dx)^2=x(xy)^2`$ has been well-known in another situation: V.I.Arnol’d found it by the investigation of slow motion of generic relaxation type equations with one fast and two slow variables . However the last case and the case studied in this paper are different. At the first place, the equation surface in the space of directions on the plane is smooth in the theory of the relaxation type equation, while it has the Whitney umbrella singularity for the implicit system case. At the second place, the plane distribution on the space of directions on the plane has singularities and is not the contact structure in the theory of the relaxation type equation, while it is the contact structure in the implicit system case. However if we put into correspondence to a relaxation type equation $`\dot{x}=\epsilon f(x,y,z),\dot{y}=\epsilon g(x,y,z),\dot{z}=h(x,y,z)+ϵr(x,y,z)`$ (where $`f,g,h,r`$ are smooth functions and $`\epsilon `$ is a small parameter) the surface $`\dot{x}f(x,y,z)=0,\dot{y}g(x,y,z)=0,h(x,y,z)=0`$, then the restriction of the projection $`(x,y,z,\dot{x},\dot{y})(x,y,\dot{x}:\dot{y})`$ to this surface is an analogue of 1-folding. In a generic case this restriction is well-defined near any critical point of the folding being here the restriction of the projection $`(x,y,z)(x,y)`$ to the surface $`h=0`$. That reduces Arnold’s case to the considered one.

We also present the complete classification of generic singularities for Clairaut implicit systems on the plane (Figure 4). An implicit system is called of Clairaut type if the system surface is smooth and, for any critical point of its folding, the velocity corresponding to this point is non-zero and lies in the image of the tangent plane to the system surface under the derivative of the folding, like the classical Clairaut equations. Under a mild condition, any implicit system of Clairaut type can be approximated by a system of Clairaut type which is foliated by smooth trajectories projecting to non-sigular solution curves. Actually we give the generic classification of such systems. See section 2.2 for details.

Remark that several singularities of the systems are studied in the parametric forms, not in the implicit forms. The parametric method provides, as well as the simple process of getting normal forms, also the clear understanding of moduli spaces of the singularities and the easy drawing of pictures of the phase portaits.

In section 2, we present the normal forms both in generic case and generic Clairaut case. The proofs is given in section 3 for the generic systems. We give the proof in section 4 for the generic Clairaut systems, using the theory of integral diagrams (cf. ).

The work of this paper started during the visit of the first author to the Department of Mathematics of Hokkaido University. He is very thankful to the department staffs for a good scientific atmosphere and nice working conditions.

## 2 Classification of singularities.

Here we present the complete list of generic point singularities of first order implicit systems on two-dimensional manifolds up to smooth orbital equivalence for both general case and Clairaut one. Since we concern just local classification, we may treat systems on the plane $`𝐑^2`$.

### 2.1 Classification in general case.

###### Theorem 2.1

() For a generic implicit system with locally bounded derivatives on the plane and for any regular point of its folding, the respective point singularity takes one of the forms listed in the second column of Table 1 near the origin up to smooth orbital equivalence. Besides the parameters of the normal form from the second column have to satisfy the respective conditions from the third one.

###### Remark 2.2

In the first column of Table 1 there is pointed the standard name for the respective singularities of generic vector fields on the plane. The smooth orbital equivalence permits smooth change of coordinates and multiplication of vector fields by smooth positive function.

###### Remark 2.3

In  the list of the respective normal forms includes also ones for resonance node with exponent $`\lambda =n𝐍`$ ( $`\dot{x}=x,\dot{y}=ny+\epsilon x^n,\epsilon \{1,0,1\}`$ ), for resonance saddle with zero coefficients by few first resonance monomials ($`\dot{x}=x[1\pm u^k+au^{2k}],\dot{y}=\lambda y,\mathrm{\hspace{0.17em}1}<k𝐍`$) and for degenerate focus without formal first intergral starting from the positive definite quadratic form ($`\dot{x}=y\pm x(r^{2k}+ar^{4k}),\dot{y}=x\pm y(r^{2k}+ar^{4k})`$ ). All these subcases are not generic. Ones can be removed by small perturbations of the implicit system.

###### Theorem 2.4

For a generic implicit system on the plane with locally bounded derivatives and any singular point of its folding the respective point singularity takes one of the forms listed in the second column of Table 2 near the origin up to smooth orbital equivalence. Besides the parameters of the normal form from the second column have to satisfy the conditions from the third one.

###### Remark 2.5

By the topological orbital equivalence normal forms of folded saddle, node and focus are $`\dot{x}=1,(\dot{y})^2=ykx^2`$ near the origin with $`k`$ equals $`1,1/20`$ and $`k=1`$, respectively .

### 2.2 Classification in the Clairaut case.

Now we denote by $`T𝐑^2`$ the tangent bundle of the plane $`𝐑^2`$, and by $`\{0\}`$ its zero section. We have the canonical projection from the tangent bundle, outside of the zero section, to the manifold of contact elements on the plane. We denote it by $`\mathrm{\Pi }:T𝐑^2\{0\}PT^{}𝐑^2`$. Here $`PT^{}𝐑^2`$ means the fiber-wise projectivization of the cotangent bundle $`T^{}𝐑^2`$ on the plane $`𝐑^2`$. This projection $`\mathrm{\Pi }`$ actually induces the $`1`$-folding mappings from surfaces in $`T𝐑^2\{0\}`$. Remark that there does not exist a canonical isomorphism between the tangent bundle $`T𝐑^2`$ and the cotangent bundle $`T^{}𝐑^2`$, but there exists the canonical isomorphism between $`PT𝐑^2`$ and $`PT^{}𝐑^2`$, as the manifold of contact elements on the plane, defined by mapping a tangential direction on $`𝐑^2`$ to the co-direction on $`𝐑^2`$ having the direction as the kernel. Then the classification under the smooth orbital equivalence of surfaces in $`T𝐑^2\{0\}`$ is reduced, via its $`1`$-folding and up to orientation of orbits, to the classification under the contact diffeomorphisms on $`PT^{}𝐑^2`$ preserving the canonical fibration $`\pi :PT^{}𝐑^2𝐑^2`$.

Let us consider a system of Clairaut type in $`T𝐑^2\{0\}`$ and its 1-unfolding in $`PT^{}𝐑^2`$. Denote by $`\mathrm{\Sigma }_c`$, the locus of contact singular points, namely, the locus on the system surface consisting of points where the contact form vanishes on the corresponding point in $`PT^{}𝐑^2`$. Also denote by $`\mathrm{\Sigma }_\pi `$ the locus of singular points on the system surface for the projection $`\pi :T𝐑^2𝐑^2`$. Then, by definition, the system is of Clairaut type if and only if the folding mapping to $`𝐑^2`$ is of corank at most one and $`\mathrm{\Sigma }_c=\mathrm{\Sigma }_\pi `$.

Dara  gives the definition of Clairaut type equations for smooth surfaces in $`𝐑^3(PT^{}𝐑^2)`$ with coordinates $`x,y,p`$, $`p=\dot{y}/\dot{x}`$, as follows: An implicit system $`G(x,y,p)=0`$ is called of Clairaut type in the sense of Dara if $`G_x+pG_y=AG+BG_p`$ holds for some function-germs $`A(x,y,p)`$ and $`B(x,y,p)`$.

By the definition, if a system is of Clairaut type in the sense of Dara, then we have $`\mathrm{\Sigma }_c=\mathrm{\Sigma }_\pi `$. In fact, in , it is proved that a non-singular system $`G(x,y,p)=0`$ is of Clairaut type in the sense of Dara if and only if the system possesses the system of complete solutions consisting of classical (smooth) solutions. In particular, each trajectory projects to a non-singular curve via the folding.

The converse is not true in general: For example, the system $`G(x,y,p)=y2p^3=0`$ is not of Clairaut type in the sense of Dara, but it satisfes the condition $`\mathrm{\Sigma }_c=\mathrm{\Sigma }_\pi `$. Moreover it has the system of complete solutions $`\mathrm{\Gamma }(t,c)=(x,y,p)=(3t^2+c,2t^3,t)`$ and each cuspidal solution curve is tangent to the discriminant $`\pi (\mathrm{\Sigma }_\pi )=\{y=0\}`$ on the $`(x,y)`$-plane. Thus it is in fact a Clauraut type in the sense of this paper. Note that in this example the singular locus of the folding mapping is defined by $`p^2=0`$ on the system surface $`\{y=2p^3\}=\{(x,2p^3,p)(x,p)(𝐑^2,0)\}`$, and the singular locus has multiple components.

A system of Clauraut type is called reduced if the Jacobian of the folding mapping has no multiple components, or more exactly, if any differentiable function vanishing on the singular locus $`\mathrm{\Sigma }_\pi `$ is devided by the Jacobian of $`\pi `$ restricted to the system surface. Then we see that any reduced system of Clauraut type can be approximated by a system of Clauraut type with the property that each trajectory projects to a non-singular curve via the folding. Then we have

###### Theorem 2.6

A generic reduced system of Clairaut type on the plane with locally bounded derivatives takes locally one of the forms in the second column of Table 3 near the origin up to smooth orbital equivalence.

The Clairaut type systems can be classified by considering parametric surfaces in $`PT^{}𝐑^2`$.

A first order differential equation germ (or, briefly, equation) is defined to be a map germ $`f:(𝐑^2,0)J^1(𝐑,𝐑)PT^{}𝐑^2.`$ We also say that $`f`$ is completely integrable if there exists a submersion germ $`\mu :(𝐑^2,0)𝐑`$ such that $`d\mu f^{}\theta =0.`$ Here $`\theta =dypdx`$ denotes the canonical contact $`1`$-form on $`J^1(𝐑,𝐑)`$. We call $`\mu `$ an independent first (or, complete) integral of $`f`$ and the pair $`(\mu ,f):(𝐑^2,0)𝐑\times J^1(𝐑,𝐑)𝐑\times PT^{}𝐑^2`$ is called a holonomic system with independent first integral. We observe that $`f|\mu ^1(t)`$ is a Legendrian immersion whose image is contained in $`\mathrm{Image}f.`$ If $`\pi f|\mu ^1(t)`$ are non-singular map for each $`t(𝐑,\mu (0)),`$ then $`\{f|\mu ^1(t)\}_{t𝐑}`$ is the family of graphs of non-singular solutions of $`\mathrm{Image}f`$ by the previous arguments. We call such a system a Clairaut type equation. These situation lead us to the following definition. Let $`(\mu ,g)`$ be a pair of a map germ $`g:(𝐑^2,0)(𝐑^2,0)`$ and a submersion germ $`\mu :(𝐑^2,0)(𝐑,0).`$ Then the diagram

$$(𝐑,0)\stackrel{𝜇}{}(𝐑^2,0)\stackrel{𝑔}{}(𝐑^2,0)$$

or briefly $`(\mu ,g)`$, is called an integral diagram if there exists an equation $`f:(𝐑^2,0)PT^{}𝐑^2`$ such that $`(\mu ,f)`$ is an equation germ with independent first integral and $`\pi f=g,`$ and we say that the integral diagram $`(\mu ,g)`$ is induced by $`f.`$ If $`f`$ is a Clairaut type equation, then $`(\mu ,\pi f)`$ is called of Clairaut type. Furthermore we introduce an equivalence relation among integral diagrams. Let $`(\mu ,g)`$ and $`(\mu ^{},g^{})`$ be integral diagrams. Then $`(\mu ,g)`$ and $`(\mu ^{},g^{})`$ are equivalent (respectively, strictly equivalent) if the diagram

$$\begin{array}{ccccc}(𝐑,0)& \stackrel{\mu }{}& (𝐑^2,0)& \stackrel{g}{}& (𝐑^2,0)\\ \kappa & & \psi & & \varphi & & \\ (𝐑,0)& \underset{\mu ^{}}{}& (𝐑^2,0)& \underset{g^{}}{}& (𝐑^2,0)\end{array}$$

commutes for some diffeomorphism germs $`\kappa ,`$ $`\psi `$ and $`\varphi `$ (respectively, $`\kappa =id_𝐑).`$

We give a generic classification of Clairaut type equations in terms of the notion of integral diagrams which implies Theorem 2.6:

###### Theorem 2.7

For a generic Clairaut type equation

$$(\mu ,f):(𝐑^2,0)𝐑\times J^1(𝐑,𝐑),$$

the integral diagram $`(\mu ,\pi f)`$ is strictly equivalent to one of germs in the following list :

(1) $`\mu =v,`$ $`g=(u,v)`$ ; Nonsingular point.

(2) $`\mu =v{\displaystyle \frac{1}{2}}u,`$ $`g=(u,v^2)`$ ; Clairaut fold.

(3) $`\mu _\alpha =v+\alpha g\mathrm{for}\alpha 𝔐_{(x,y)},`$ $`g=(u,v^3+uv)`$ ; Clairaut cusp.

(4) $`\mu =v{\displaystyle \frac{1}{2}}u^2,`$ $`g=(u,{\displaystyle \frac{1}{4}}v^2)`$ ; Clairaut Whitney umbrella.

Note that the forms (1) (2) and (4) of Theorem 2.6 are obtained from those in Theorem 2.7 as follows:

(1) : The equation is given by $`f=(u,v,0)`$, namely, by $`p=0`$. Setting $`\dot{x}=1`$, we have $`\dot{y}=0`$.

(2) : The equation is given by $`f=(u,v^2,v)`$, namely, by $`p^2=y`$. Setting $`\dot{x}=1`$, we have $`(\dot{y})^2=y`$.

(4) : The equation is given by $`f=(u,{\displaystyle \frac{1}{4}}v^2,{\displaystyle \frac{1}{2}}uv)`$, namely by $`p^2=x^2y`$. Setting $`\dot{x}=1`$, we have $`(\dot{y})^2=x^2y`$.

The form (3) of Theorem 2.6 is obtained from the parametric form (3) of Theorem 2.7 as follows: The the $`1`$-folding map is given by $`(u,v)(x,y,p)=(u,v^3+uv,h(u,v))`$, where $`p={\displaystyle \frac{\dot{y}}{\dot{x}}}`$ and

$$h(u,v)=\frac{v(3v^2+u)\frac{\alpha }{x}(u,v^3+uv)}{1+(3v^2+u)\frac{\alpha }{y}(u,v^3+uv)}.$$

Then $`u`$ and $`h`$ form another parametrization, and we have the implicit form $`y=v^3+uv=v(x,p)^3+xv(x,p)=:\psi (x,p)`$. By using a cusp-preserving diffeomorphism which maps the phase curve passing through the cusp vertex to the $`x`$-axis, we can suppose that $`\psi (x,0)=0`$. So we have $`y=\psi (x,p)=p\phi (x,p)`$. By setting $`\dot{x}=1`$, we have $`y=\dot{y}\phi (x,\dot{y})`$.

Theorem 2.7 gives a generic classification of integral diagrams of Clairaut type under the strict equivalence. We remark that each germs from (1) to (4) are not equivalent as integral diagrams. Thus the problem is reduced to classify germs which are contained in the family (3) under the equivalence. This family is parametrized by function germs $`\alpha `$ which are called functional moduli. In , it is given a characterization of functional moduli relative to the equivalence. For a map $`g:x=u,y=v^3+uv`$, define in the $`x,y`$-plane the set $`\mathrm{\Delta }`$ of points where this map has three different preimages, and the set $`𝒟`$ like the boundary of the set $`\mathrm{\Delta }`$. In  Dufour has shown the following characterization theorem.

###### Theorem 2.8

Let $`(\mu _\alpha ,g)`$ be an integral diagram of the type (3) in Theorem 2.7. Then for any $`\alpha ,`$ there exists a function germ $`\alpha ^{}:(𝐑^2,0)(𝐑,0)`$ such that

(1) $`(\mu _\alpha ,g)`$ is equivalent to $`(\mu _\alpha ^{},g).`$

(2) $`\alpha ^{}|𝒟=0.`$

Dufour has also shown that the uniqueness of functional moduli relative to the equivalence. We say that $`\alpha `$ and $`\alpha ^{}`$ are equivalent as moduli if there exists $`a𝐑\{0\}`$ such that $`a\alpha (x,y)=\alpha ^{}(a^2x,a^3y)`$ for any $`(x,y)\mathrm{\Delta }.`$ We remark that the above definition of the equivalence among functional moduli is slightly different from Dufour’s original definition in . If we adopt his definition, we cannot assert the necessity of the condition that functional moduli are equivalent. Actually, in , we have introduced the above definition and shown the following theorem:

###### Theorem 2.9

Let $`(\mu _\alpha ,g)`$ and $`(\mu _\alpha ^{},g)`$ be integral diagrams of (3) such that $`\alpha |𝒟=\alpha ^{}|𝒟=0.`$ Then $`(\mu _\alpha ,g)`$ and $`(\mu _\alpha ^{},g)`$ are equivalent as integral diagram if and only if $`\alpha `$ and $`\alpha ^{}`$ are equivalent as moduli.

This theorem asserts that the equivalence classes of functional moduli $`\alpha `$ with $`\alpha |𝒟=0`$ are the complete invariant for generic classifications of Clairaut type equations under the equivalence relation given by the group of point transformations.

We define $`𝔐(𝒟)=\{\alpha 𝔐_{(x,y)}\alpha |𝒟=0\}`$ and $`_{\mathrm{cusp}}=𝔐(𝒟)/,`$ where $``$ denotes the equivalence relation as moduli. The above theorem asserts that the moduli space for generic Clairaut type equation is $`_{\mathrm{cusp}}.`$

###### Remark 2.10

In the paper , it is given the classification of first order implicit differential equations in $`PT^{}𝐑^2`$ endowed with independent first integrals. In , the pleated singular points are called the regular cusps and their normal forms are give in the parametric forms: $`(u,v)(u^3+uv,v)`$, the parametrization of the folding map, with the first integral $`\mu ={\displaystyle \frac{3}{4}}u^4+{\displaystyle \frac{1}{2}}u^2v+\alpha (u^3+uv,v)`$, where $`\alpha `$ is any function with $`\alpha (0,0)=0,{\displaystyle \frac{\alpha }{y}}(0,0)=\pm 1`$ (Theorem B (5) of ). Moreover, by a theorem of Kurokawa (Theorem A of ), we can take, up to the equivalence, the functional moduli $`\alpha `$ satisfying

$$\alpha (0,y)=\pm y+\frac{1}{2}\frac{^2\alpha }{y^2}(0,0)y^2,(y0).$$

Furthermore the sign $`\pm 1`$ and $`\chi _\alpha ={\displaystyle \frac{^2\alpha }{y^2}}(0,0)`$ are invariants of the equation (cf. ).

Then the $`1`$-folding map is given by $`(u,v)(x,y,p)=(u^3+uv,v,k(u,v))`$, where $`p={\displaystyle \frac{\dot{y}}{\dot{x}}}`$ and

$$k(u,v)=\frac{u+\frac{\alpha }{x}(u^3+uv,v)}{\frac{1}{2}u^2\frac{\alpha }{y}(u^3+uv,v)}.$$

Remark that functions $`v`$ and $`k`$ provide another parametrization of the system surface. Then we have $`x=u^3+uv=u(v,k)^3+u(v,k)v=u(y,p)^3+u(y,p)y=:\psi (y,p)`$. Note that the locus $`p=0`$ defines a smooth curve that is tangent to the $`y`$-axis at the cusp vertex on $`(x,y)`$-plane. By using a cusp-preserving diffeomorphism which maps this curve to the $`y`$-axis, we may suppose $`\psi (y,0)=0`$. Thus we have the form $`x=p\phi (y,p)`$. Then we have the normal form in table 2 by setting $`\dot{x}=1`$.

The first integral $`\mu _\alpha ={\displaystyle \frac{3}{4}}u^4+{\displaystyle \frac{1}{2}}u^2v+\alpha (u^3+uv,v)`$ for the system with the folding map $`(u,v)(u^3+uv,v)`$ is transformed to the first integral $`\mu _\alpha ^{}={\displaystyle \frac{3}{4}}u^4+{\displaystyle \frac{1}{2}}u^2v+\alpha ^{}(u^3+uv,v)`$ for the system with the same folding map $`(u,v)(u^3+uv,v)`$ if and only if $`\alpha ^{}(x,y)=\alpha (x,y)`$ or $`\alpha ^{}(x,y)=\alpha (x,y)`$ on the cuspidal open domain $`\mathrm{\Delta }:=\{y^3+{\displaystyle \frac{27}{4}}x^2<0\}`$ near the origin (Proposition 6.1 of ).

The exact description of the moduli space of pleated singular points remains open. However the reason of the existence of the functional moduli clearly comes from that for the $`3`$-webs on the plane by solution curves, and we see, as mentioned above, the moduli space is dominated by the function space

$$\left\{\alpha \right|\overline{\mathrm{\Delta }}|\alpha :(𝐑^2,0)(𝐑,0),\alpha (0,y)=\pm y+\frac{1}{2}\frac{^2\alpha }{y^2}(0,0)y^2,(y0)\}.$$

## 3 Singularities in a general case.

Here we prove Theorem 2.1, 2.4. At the first place we study typical singularities of system (1-)folding, and then on that base the theorems are proved.

### 3.1 Generic singularities of 1-folding.

###### Proposition 3.1

For a generic implicit system its surface is either empty or a smooth two dimensional submanifold of the tangent bundle space.

Proposition 3.1 follows immediately from Thom transversality theorem.

###### Proposition 3.2

The folding of a generic implicit system with locally bounded derivatives is $`LR`$-stable map and it can have singular points of type Whitney either fold or pleat only. In other words, near any its critical point, the folding takes respectively the form

$$\text{either}\{\begin{array}{c}x=u\hfill \\ y=v^2\hfill \end{array}\text{or}\{\begin{array}{c}x=u\hfill \\ y=v^3+uv\hfill \end{array}$$

in appropriate smooth coordinates near this point and its image under the folding with the origin at them.

Recall that the system folding is the restriction of the tangent bundle projection to the system surface. Thus due to Goryunov theorem such restriction in a generic case can have all generic singularities as a generic map between $`n`$-dimensional manifolds with the dimension of the kernel being no greater then the dimension of the kernel of this projection . But the last dimension is also equal to $`n`$, and that permits all generic singularities between the system surface and the phase space.

In the two dimensional case these singularities are Whitney fold and pleat , . Besides our systems are with locally bounded derivatives. Therefore system foldings are proper maps. Consequently, for a generic system its folding is $`LR`$-stable map , .

Thus Proposition 3.2 is true.

###### Remark 3.3

A map is called $`LR`$-stable (= left right stable) if for any map being sufficiently close to it these two maps can be carried one to another by diffeomorphisms of the image space and the preimage space which are close to the identities. For example the map $`(x,y)x`$ from the circle $`x^2+y^2=1`$ to the $`x`$-axis is $`LR`$-stable (under small perturbations in $`C^k`$-topology with $`k2`$.

###### Proposition 3.4

For a generic implicit system with locally bounded derivatives any critical point of its folding does not belong to the zero section of the tangent bundle to the phase space. In particular, near such a point the system 1-folding is well defined.

This proposition follows immediately from Thom transversality theorem because the conditions

$$F=0,rankF_{\dot{x}}<n,\dot{x}=0$$

defines in the jet space of systems the closed Whitney stratified manifold of the codimension $`2n+1`$ which is greater then the dimension $`2n`$ of the tangent bundle. Therefore for a generic system these conditions can not be satisfied semiltaniously.

###### Theorem 3.5

For a generic implicit system with locally bounded derivatives any critical point of its folding is either regular point or critical point of type Whitney umbrella for the 1-folding of this system. Besides the image of the set of critical points of the system folding under the 1-folding is generically replaced with respect to the standard contact structure in the space of directions on the phase space. In particular, it can have tangency with this structure only of the first order and only at the points being regular for the system 1-folding and critical one of type Whitney fold for the system folding.

###### Remark 3.6

A generic replacement of the image with respect to the direction axis in the projectivization of the tungent bundle follows immediately from Proposition 3.2 and the following corollary of the first statement of Theorem 3.5.

###### Corollary 3.7

For a generic implicit system with locally bounded derivatives any critical point of type Whitney umbrella this system 1-folding is the critical point of type Whitney fold for this system folding.

Theorem 3.5 also implies immediately

###### Corollary 3.8

For a generic implicit system with locally bounded derivatives its point singularities provided by the critical point of this system folding are described either by generic singularities of first order implicit differential equation or by such an equation provided by a generic replacement in the space of directions on the plane of germ of the Whitney umbrella at it vertex.

Let us prove Theorem 3.5. Let the surface of a generic implicit system with locally bounded derivatives be not empty.

For such a system its folding is $`LR`$-stable due to Proposition 3.2 and in the strength of Proposition 3.4 the set of critical points of this folding does not intersect the zero section of the tangent bundle to the phase space. Therefore for any system being sufficiently close to the given one the 1-folding is well defined near this set.

Note that any regular point of this system folding is also regular point of its 1-folding because the folding provides two components of the 1-folding.

Now we again can apply Goryunov theorem . Due to this theorem the 1-folding of a generic system can have all generic singularities like a generic map from two-dimensional manifold to 3-dimensional one. But any critical point of such a map is of Whitney umbrella type. Taking into account Corollary 3.7 we find that near such a point of a generic system this system 1-folding takes form

$$x=u,y=v^2,z=uv$$

in appropriate smooth coordinates near this critical point and local smooth coordinates in the image space fibered over the phase space ($`x,y`$-space) with the origins at this point and its image, respectively.

A typical replacement (of the image of the set of critical points of the folding under the 1-folding with respect to the standard contact structure in the space of directions on the phase space) can be obtained by small rotations of the tangent planes. Really such a rotations provides small perturbation of the system but they do not change the set of critical value of the system folding and can supply any small rotation of the field of direction defined by our system on this set. Finally if for a generic system the ”replacement” is typical then it is also typical for any system being sufficiently close to the chosen one due to $`LR`$-stability of the folding of a generic system in the strength of Proposition 3.2.

Theorem 3.5 is proved.

### 3.2 Proofs of Theorems 2.1, 2.4.

At the first place we prove Theorem 2.1 and then Theorem 2.4.

Let $`P`$ be a regular point of the folding of a generic system. Near such a point this system surface is smooth section of the tangent bundle. This section provides smooth vector field $`v`$ near the image $`\stackrel{~}{P}`$ of this point under the system folding.

If the point $`P`$ does not belong to the zero section of the tangent bundle then this field does not vanish at the point $`\stackrel{~}{P}`$. In that case the germ of the field $`v`$ at this point is $`C^{\mathrm{}}`$-diffeomorphic to the germ of the constant vector field $`(1,0)`$ at the origin . That gives the first singularity from Table 1.

If the point $`P`$ belongs to zero section of the tangent bundle then field $`v`$ vanishes at the point $`\stackrel{~}{P}`$. Small perturbations of the studied system implies smooth small changing of the field near the point $`P`$ due to $`LR`$-stability of this system folding in accordance with Proposition 3.4. But for a generic system this map (”small perturbation” $``$ ”small perturbations $`v`$ near $`\stackrel{~}{P}`$”) is continuous and small perturbations of the system provides all small perturbations of the system surface near the point $`P`$ and, hence, all small perturbations of the field $`v`$ near the point $`\stackrel{~}{P}`$. Therefore for a generic system vector field $`v`$ has at a point $`\stackrel{~}{P}`$ a generic singular point. Now the rest part of Theorem 2.1 follows from the classical results about normal forms of generic vector fields near singular points up to smooth orbital equivalence . Theorem 2.1 is proved.

Let us prove Theorem 2.4. Due to Corollary 3.8 local singularities of a generic implicit system at singular points of its folding are described by generic singularities of first order implicit ODE except the singular points of this system 1-folding of type Whitney umbrella. But here we need to take into account that a critical point of a generic implicit system never belongs the zero section of the tangent bundle due to Proposition 3.4. Consequently, generic singularities of implicit equations of type folded regular point, folded singular point and pleated point gives here first five singularities and the last one from Table 2.

Now due to Corollary 3.8 to finish the proof one need to get the normal form of an implicit first order ODE provided by a generic replacement of Whitney umbrella to the space of directions on the plane. For a generic system such a replacement has to have the following properties. At the Whitney umbrella vertex both the the contact plane and the vertical direction do not tangent to the image under a generic system 1-folding of set of critical points of this system folding and the ”handle” of this umbrella.

Consequently, smooth local coordinate systems $`u,v`$ on the system surface and $`x,y`$ on the phase space with the origins at the studied point and its image under the system folding can be choosen such that this system 1-folding takes the form

$$x=v^2,y=u,\frac{dy}{dx}=h(u,v)$$
(1)
where $`h`$ is a smooth function, $`h(0,0)=0`$, and $`dy/dx`$ is local coordinate along the direction axis and also in these coordinates the ”handle” of the Whitney umbrella is over the line $`xy=0`$.

Near the origin on the $`y`$-axis and the line $`\{xy=0\}\{x>0\},`$ the studied implicit first order equation gives two smooth direction fields $`dy/dx=f_1(y)`$ and $`dy/dx=f_2(y)`$ respectively, where $`f_1,f_2`$ are smooth functions; $`f_1(0)=f_2(0)=0`$ because $`h(0,0)=0`$. Due to Hadamard lemma these functions can be presented in the form $`f_1(y)=y\stackrel{~}{f}_1(y),f_2(y)=y\stackrel{~}{f}_2(y)`$ where $`\stackrel{~}{f}_1,\stackrel{~}{f}_2`$ are smooth functions. Near the origin the direction field

$$dy/dx=(yx)\stackrel{~}{f}_1(y)+x\stackrel{~}{f}_2(y)$$

provides the semiltanious smooth extension of these two fields.

Near the origin the extended field has a first integral of the form $`y+xI_1(x,y)`$ where $`I_1`$ is smooth function. Taking this integral and the function $`u+v^2I_1(v^2,u)`$ like new coordinates $`y`$ and $`u`$, respectively, one preserves the first two forms of the equation (1) but in new coordinates the function $`h`$ is identically zero on the $`u`$-axis and on the set corresponding to the ”handle”. In new coordinates near the origin the last set can be defined by some equation $`uv^2X(v^2)=0`$ where $`X`$ is a smooth function, $`X(0)>0`$. Due to Hadamard lemma the function $`h`$ in the last equation from (1) can be written in the form $`h(u,v)=v(uv^2X(v^2))H(u,v)`$ where $`H`$ is smooth function; $`H(0,0)0`$ because at the studied point the system 1-folding has singularity of type Whitney umbrella. Consequently near the origin the rescalings $`\stackrel{~}{v}=v\sqrt{X(v^2)}`$ and $`\stackrel{~}{x}=xX(x)`$ reduce the system (1) of equations to the form (”tilde” in notations of new coordinates is omitted)

$$x=v^2,y=u,\frac{dy}{dx}=v(uv^2)H(u,v)$$
(2)
with some new smooth function $`H`$ no vanishing at the origin. It is easy to see that the direction field provided by the last equation can be lifted to the smooth direction field

$$du/dv=2v^2(uv^2)H(u,v).$$
(3)
The last direction field has first integral of the form

$$I(u,v)=u+v^3J_1(u,v)$$

as it is easy to see.

Now it is sufficient to get normal form of this integral by the changing of coordinates commuting with the involution $`(u,v)(u,v)`$ defined by the system folding. Taking new coordinate $`u`$ in the form $`(I(u,v)+I(u,v))/2`$ (being even with respect to $`v`$ and so it is permitted by our involution) we reduce the integral to the form $`I(u,v)=u+v^3J_1(u,v^2)`$ or

$$I(u,v)=u+v^3a(u)+v^5b(u)+v^7c(u,v^2),$$

where $`a,b,c`$ are some smooth functions; $`a(0)=0a^{}(0)b(0)`$ because $`v^2u`$ is the term of the lower degree in the right hand side of the equation (3), $`b(0)0`$ due to $`H(0,0)0.`$

The following lemma completes the proof

###### Lemma 3.9

(, ) Near the origin the $`(u,v)`$-plane a function $`u+v^3a(u)+v^5b(u)+v^7c(u,v^2)`$ with smooth functions $`a,b`$ and $`c`$, $`a(0)=0a^{}(0)b(0)`$, is reduced to the form $`u+v^3u+v^5`$ by a smooth diffeomorphism preserving the origin and commuting with the involution $`(u,v)(u,v)`$.

Hence the function $`H`$ in the equation (3) can be reduced to $`1`$. After that this equation takes the form $`du/dv=2v^2(uv^2)`$. Near the points under consideration that implies the equation $`dy/dx=v(uv^2)`$ on the system surface or the equation $`(dy/dx)^2=x(yx)^2`$ on the $`(x,y)`$-plane. That gives the rest (fifth) singularity from the list of Table 2. Theorem 2.4 is proved.

## 4 Clairaut type equations and Legendre singularity theory.

### 4.1 Legendrian unfoldings.

We briefly review the theory of one-parameter Legendrian unfoldings. We now consider the $`1`$-jet bundle $`J^1(𝐑\times 𝐑,𝐑)`$ and the canonical $`1`$-form $`\mathrm{\Theta }`$ on the space. Let $`(t,x)`$ be the canonical coordinate on $`𝐑\times 𝐑`$ and $`(t,x,y,q,p)`$ be the corresponding coordinate on $`J^1(𝐑\times 𝐑,𝐑).`$ Then the canonical $`1`$-form is given by $`\mathrm{\Theta }=dypdxqdt=\theta qdt.`$ We also have the natural projection

$$\mathrm{\Pi }:J^1(𝐑\times 𝐑,𝐑)𝐑\times 𝐑\times 𝐑$$

defined by $`\mathrm{\Pi }(t,x,y,q,p)=(t,x,y).`$ We call the above $`1`$-jet bundle an unfolded $`1`$-jet bundle. Let $`(\mu ,f)`$ be an equation with complete integral. Then there exists a unique function germ $`h:(𝐑^2,0)𝐑`$ such that $`f^{}\theta =hd\mu .`$ Define a map germ

$$\mathrm{}_{(\mu ,f)}:(𝐑^2,0)J^1(𝐑\times 𝐑,𝐑)$$

by

$$\mathrm{}_{(\mu ,f)}(u)=(\mu (u),xf(u),yf(u),h(u),pf(u)).$$

Then we can easily show that if $`(\mu ,f)`$ is a Clairaut type equation, then $`\mathrm{}_{(\mu ,f)}`$ is a Legendrian immersion germ. We call $`\mathrm{}_{(\mu ,f)}`$ a complete Legendrian unfolding associated with $`(\mu ,f).`$ By the aid of the notion of Legendrian unfoldings, Clairaut type equations are characterized as follows :

###### Proposition 4.1

Let $`(\mu ,f):(𝐑^2,0)𝐑\times J^1(𝐑,𝐑)PT^{}𝐑^2`$ be an equation with complete integral. Then $`(\mu ,f)`$ is a Clairaut type equation if and only if $`\mathrm{}_{(\mu ,f)}`$ is a Legendrian non-singular Legendrian immersion germ.

A complete Legendrian unfolding $`\mathrm{}_{(\mu ,f)}`$ associated to $`(\mu ,f)`$ is called a Legendrian unfolding of Clairaut type if $`\mathrm{}_{(\mu ,f)}`$ is a Clairaut type equation.

### 4.2 Genericity.

Returning to the study of equations with complete integral, we now establish the notion of the genericity.

Let $`U𝐑^2`$ be an open set. We denote by $`\mathrm{Int}(U,𝐑\times J^1(𝐑,𝐑))`$ the set of systems with complete integral $`(\mu ,f):U𝐑\times J^1(𝐑,𝐑).`$ We also define $`L(U,J^1(𝐑\times 𝐑,𝐑))`$ to be the set of complete Legendrian unfoldings $`\mathrm{}_{(\mu ,f)}:UJ^1(𝐑\times 𝐑,𝐑).`$

These sets are topological spaces equipped with the Whitney $`C^{\mathrm{}}`$-topology. A subset of either spaces is said to be generic if it is an open dense subset in the space.

The genericity of a property of germs are defined as follows. Let $`P`$ be a property of equation germs with complete integral $`(\mu ,f):(𝐑^2,0)𝐑\times J^1(𝐑,𝐑)`$ (respectively, Legendrian unfoldings $`\mathrm{}_{(\mu ,f)}:(𝐑^2,0)J^1(𝐑\times 𝐑,𝐑)).`$ For an open set $`U𝐑^2,`$ we define $`𝒫(U)`$ to be the set of $`(\mu ,f)\text{Int}(U,𝐑\times J^1(𝐑,𝐑))`$ (respectively, $`\mathrm{}_{(\mu ,f)}L(U,J^1(𝐑\times 𝐑,𝐑)))`$ such that the germ at $`x`$ whose representative is given by $`(\mu ,f)`$ (respectively, $`\mathrm{}_{(\mu ,f)}`$) has property $`P`$ for any $`xU.`$

The property $`P`$ is said to be generic if for some neighbourhood $`U`$ of $`0`$ in $`𝐑^2,`$ the set $`𝒫(U)`$ is a generic subset in $`\text{Int}(U,𝐑\times J^1(𝐑,𝐑))`$ (respectively, $`L(U,J^1(𝐑\times 𝐑,𝐑)).`$

By the construction, we have a well-defined continuous mapping

$$(\mathrm{\Pi }_1)_{}:L(U,J^1(𝐑\times 𝐑,𝐑))\text{Int}(U,𝐑\times J^1(𝐑,𝐑))$$

defined by $`(\mathrm{\Pi }_1)_{}(\mathrm{}_{(\mu ,f)})=\mathrm{\Pi }_1\mathrm{}_{(\mu ,f)}=(\mu ,f),`$ where $`\mathrm{\Pi }_1:J^1(𝐑\times 𝐑,𝐑)J^1(𝐑,𝐑)`$ is the canonical projection. Then it has been shown the following fundamental theorem:

###### Theorem 4.2

The continuous map

$$(\mathrm{\Pi }_1)_{}:L(U,J^1(𝐑\times 𝐑,𝐑))\text{Int}(U,𝐑\times J^1(𝐑,𝐑))$$

is a homeomorphism.

As in , we define the equivalence relation among parametric systems under the group of point transformations: Two equations $`f,f^{}:(𝐑^2,0)PT^{}𝐑^2`$ are equivalent under the group of point transformations if there exists a diffeomorphism $`\varphi :(𝐑^2,\pi (f(0)))(𝐑^2,\pi (f^{}(0)))`$ such that the canonical lifting $`\widehat{\varphi }:(PT^{}𝐑^2,f(0))(PT^{}𝐑^2,f^{}(0))`$ transforms $`f`$ to $`f^{}`$, namely, $`\widehat{\varphi }f=f^{}\psi `$ for some diffeomorphism $`\psi :(𝐑^2,0)(𝐑^2,0)`$.

We need the following basic result:

###### Proposition 4.3

Let $`f,f^{}:(𝐑^2,0)PT^{}𝐑^2`$ be completely integrable equations with independent first integrals $`\mu ,\mu ^{}:(𝐑^2,0)(𝐑,0)`$ respectively. Assume $`\pi f`$ and $`\pi f^{}:(𝐑^2,0)(𝐑^2,0)`$ have nowhere dense singular sets. Then $`f`$ and $`f^{}`$ are equivalent under the group of point transformations if and only if the induced integral diagrams $`(\mu ,\pi f)`$ and $`(\mu ^{},\pi f^{})`$ are equivalent.

Proof.  Assume $`(\mu ,\pi f)`$ and $`(\mu ^{},\pi f^{})`$ are equivalent by diffeomorphisms $`(\kappa ,\psi ,\varphi )`$. Then $`\varphi `$ maps integral curves $`\pi (f(\mu ^1(\mu (u,v))))`$ through $`\pi (f(u,v))`$ to $`\pi (f^{}(\mu _{}^{}{}_{}{}^{1}(\kappa (\mu (u,v)))))=\pi (f^{}(\mu _{}^{}{}_{}{}^{1}(\mu ^{}(\psi (u,v)))))`$ through $`\pi (f^{}(\psi (u,v)))`$, so the tangent lines to them. Since the set of contact singular points is contained in the set of critical points of the projection $`\pi `$, we see $`\widehat{\varphi }f=f^{}\psi `$. This implies $`f`$ and $`f^{}`$ are equivalent under the group of point transformations. The converse implication is clear. $`\mathrm{}`$

## 5 Proofs for Theorem 2.7 and Theorem 2.6.

In the case when $`\mathrm{}_{(\mu ,f)}`$ is a Legendrian immersion germ, there exists a generating family of $`\mathrm{}_{(\mu ,f)}`$ by the Arnol’d-Zakalyukin’s theory (). In this case the generating family is naturally constructed by a one-parameter family of generating families associated with $`(\mu ,\mathrm{}).`$ Let $`F:((𝐑\times 𝐑)\times 𝐑^k,0)(𝐑,0)`$ be a function germ such that $`d_2F|0\times 𝐑\times 𝐑^k`$ is non-singular, where $`d_2F(t,x,q)=(\frac{F}{q_1}(t,x,q),\mathrm{},\frac{F}{q_k}(t,x,q)).`$ We call $`F`$ a Morse family. Then $`C(F)=d_2F^1(0)`$ is a smooth surface germ and $`\pi _F:(C(F),0)𝐑`$ is a submersion germ, where $`\pi _F(t,x,q)=t.`$ We call the submanifold $`C(F)`$ a catastrophe set of  $`F.`$ Define

$$\stackrel{~}{\mathrm{\Phi }}_F:(C(F),0)J^1(𝐑,𝐑)$$

by

$$\stackrel{~}{\mathrm{\Phi }}_F(t,x,q)=(x,F(t,x,q),\frac{F}{x}(t,x,q))$$

and

$$\mathrm{\Phi }_F:(C(F),0)J^1(𝐑\times 𝐑,𝐑)$$

by

$$\mathrm{\Phi }_F(t,x,q)=(t,x,F(t,x,q),\frac{F}{t}(t,x,q),\frac{F}{x}(t,x,q)).$$

Since $`\frac{F}{q_i}=0`$ on $`C(F),`$ we can easily show that $`(\stackrel{~}{\mathrm{\Phi }}_F)^{}\theta =\frac{F}{t}|C(F)dt|C(F)=0.`$ By definition, $`\mathrm{\Phi }_F`$ is a Legendrian unfolding associated with the Legendrian family $`(\pi _F,\stackrel{~}{\mathrm{\Phi }}_F).`$ By the same method of the theory of Arnol’d-Zakalyukin (), we can show the following proposition.

###### Proposition 5.1

All Legendrian unfolding germs are constructed by the above method.

Let $`(\mu ,f)`$ be a Clairaut type equation. By Proposition 4.1, $`\mathrm{}_{(\mu ,f)}`$ is a Legendrian immersion. Then we can choose a family of function germs

$$F:(𝐑\times 𝐑,0)(𝐑,0)$$

such that $`\mathrm{Image}j^1F_t=f(\mu ^1(t))`$ for any $`t𝐑,`$ where $`F_t(x)=F(t,x).`$ We remark that the map germ

$$j_1^1F:(𝐑\times 𝐑,0)J^1(𝐑,𝐑)$$

defined by $`j_1^1F(t,x)=j^1F_t(x)`$ is not necessary an immersion germ. In this case we have $`(C(F),0)=(𝐑\times 𝐑,0)`$ and

$$\mathrm{\Phi }_F=j^1F:(𝐑\times 𝐑,0)J^1(𝐑\times 𝐑,𝐑),$$

so that it is a complete Legendrian unfolding associated with $`(\pi _1,j_1^1F).`$ Thus the generating family of a Legendrian unfolding of Clairaut type is given by the above germ.

In order to prove Theorem 2.7, we now introduce equivalence relations among Legendrian unfoldings. Let $`(\mu ,g)`$ and $`(\mu ^{},g^{})`$ be integral diagrams. Then $`(\mu ,g)`$ and $`(\mu ^{},g^{})`$ are $`^+`$-equivalent if there exist a diffeomorphism germ $`\mathrm{\Psi }:(𝐑\times (𝐑\times 𝐑),0)(𝐑\times (𝐑\times 𝐑),0)`$ of the form $`\mathrm{\Psi }(t,x,y)=(t+\alpha (x,y),\psi (x,y))`$ and a diffeomorphism germ $`\mathrm{\Phi }:(𝐑^2,0)(𝐑^2,0)`$ such that $`\mathrm{\Psi }(\mu ,g)=(\mu ^{},g^{})\mathrm{\Phi }.`$ We remark that if $`(\mu ,g)`$ and $`(\mu ^{},g^{})`$ are $`^+`$-equivalent by the above diffeomorphisms, then we have $`\mu (u)+\alpha g(u)=\mu ^{}\mathrm{\Phi }(u)`$ and $`\psi g(u)=g^{}\mathrm{\Phi }(u)`$ for any $`u(𝐑^2,0).`$ Thus the diagram $`(\mu +\alpha g,g)`$ is strictly equivalent to $`(\mu ^{},g^{}).`$

We now define the corresponding equivalence relation among Legendrian unfoldings. Let $`\mathrm{}_{(\mu ,f)},\mathrm{}_{(\mu ^{},f^{})}:(𝐑^2,0)(J^1(𝐑\times 𝐑,𝐑),z_0)`$ be Legendrian unfoldings. We say that $`\mathrm{}_{(\mu ,f)}`$ and $`\mathrm{}_{(\mu ^{},f^{})}`$ are $`S.P^+`$-Legendrian equivalent (respectively, $`S.P`$-Legendrian equivalent) if there exist a contact diffeomorphism germ $`K:(J^1(𝐑\times 𝐑,𝐑),z_0)(J^1(𝐑\times 𝐑^n,𝐑),z_0^{}),`$ a diffeomorphism germ $`\mathrm{\Phi }:(𝐑^2,0)(𝐑^2,0)`$ and a diffeomorphism germ $`\mathrm{\Psi }:(𝐑\times (𝐑\times 𝐑),\mathrm{\Pi }(z_0))(𝐑\times (𝐑\times 𝐑),\mathrm{\Pi }(z_0^{}))`$ of the form $`\mathrm{\Psi }(t,x,y)=(t+\alpha (x,y),\psi (x,y))`$ (respectively, $`\mathrm{\Psi }(t,x,y)=(t,\psi (x,y))`$) such that $`\mathrm{\Pi }K=\mathrm{\Psi }\mathrm{\Pi }`$ and $`K=^{}\mathrm{\Phi }.`$ It is clear that if $`\mathrm{}_{(\mu ,f)}`$ and $`\mathrm{}_{(\mu ^{},f^{})}`$ are $`S.P^+`$-Legendrian equivalent (respectively, $`S.P`$-Legendrian equivalent), then $`(\mu ,\pi f)`$ and $`(\mu ^{},\pi f^{})`$ are $`^+`$-equivalent (respectively, strictly equivalent). By Theorem 1.1 in , the converse is also true for generic $`(\mu ,f)`$ and $`(\mu ^{},f^{}).`$ The notion of the stability of Legendrian unfoldings with respect to $`S.P^+`$-Legendrian equivalence (respectively, $`S.P`$-Legendrian equivalence) is analogous to the usual notion of the stability of Legendrian immersion germs with respect to Legendrian equivalence (cf. Part $`III`$ in ).

On the other hand, we can interpret the above equivalence relation in terms of generating families. For the purpose, we use some notations and results in . Let $`\stackrel{~}{F},\stackrel{~}{G}:(𝐑\times (𝐑\times 𝐑),0)(𝐑,0)`$ be generating families of Legendrian unfoldings of Clairaut type. We say that $`\stackrel{~}{F}`$ and $`\stackrel{~}{G}`$ are P-$`𝒞^+`$-equivalent (respectively, P-$`𝒞`$-equivalent) if there exists a diffeomorphism germ $`\mathrm{\Phi }:(𝐑\times (𝐑\times 𝐑),0)(𝐑\times (𝐑\times 𝐑),0)`$ of the form $`\mathrm{\Phi }(t,x,y)=(t+\alpha (x,y),\varphi _1(x,y),\varphi _2(x,y))`$ (respectively, $`\mathrm{\Phi }(t,x,y)=(t,\varphi _1(x,y),\varphi _2(x,y))`$) such that $`F\mathrm{\Phi }_{_{(t,x,y)}}=G_{_{(t,x,y)}}`$ where $`G_{_{(t,x,y)}}`$ is the ideal generated by $`G`$ in the local ring of function germs $`_{(t,x,y)}`$ of $`(t,x,y)`$-variables. We also say that $`\stackrel{~}{F}(t,x,y)`$ is $`𝒞^+`$ (respectively, $`𝒞`$)-versal deformation of $`f=F|𝐑\times 0`$ if

$$_t=\frac{df}{dt}_𝐑+f__t+\frac{F}{x}|𝐑\times \{(0,0)\},\frac{F}{y}|𝐑\times \{(0,0)\},1_𝐑$$

(respectively,

$$_t=f__t+\frac{F}{x_1}|𝐑\times \{0\},\mathrm{},\frac{F}{x_n}|𝐑\times \{0\},1_𝐑).$$

By the similar arguments like as those of Theorems 20.8 and 21.4 in , we can show the following :

###### Theorem 5.2

Let $`\stackrel{~}{F},\stackrel{~}{G}:(𝐑\times (𝐑\times 𝐑),0)(𝐑,0)`$ be generating families of Legendrian unfoldings of Clairaut type $`\mathrm{\Phi }_F,\mathrm{\Phi }_G`$ respectively. Then

(1) $`\mathrm{\Phi }_F`$ and $`\mathrm{\Phi }_G`$ are S.P<sup>+</sup> (respectively, S.P)-Legendrian equivalent if and only if $`\stackrel{~}{F}`$ and $`\stackrel{~}{G}`$ are P-$`𝒞^+`$ (respectively, $`𝒞`$)-equivalent.

(2) $`\mathrm{\Phi }_F`$ is S.P<sup>+</sup> (respectively, S.P)-Legendrian stable if and only if $`\stackrel{~}{F}`$ is a P-$`𝒞^+`$ (respectively, $`𝒞`$)-versal deformation of $`f=F|𝐑\times \{0\}.`$

The following theorem is a corollary of Damon’s general versality theorem in .

###### Theorem 5.3

Let $`\stackrel{~}{F},\stackrel{~}{G}:(𝐑\times (𝐑\times 𝐑),0)(𝐑,0)`$ be generating families of Legendrian unfoldings of Clairaut type such that $`\mathrm{\Phi }_F,\mathrm{\Phi }_G`$ are S.P<sup>+</sup> (respectively, S.P)-Legendrian stable. Then $`\mathrm{\Phi }_F,\mathrm{\Phi }_G`$ are S.P<sup>+</sup> (respectively, S.P)-Legendrian equivalent if and only if $`f=F|𝐑\times \{0\},g=G|𝐑\times \{0\}`$ are $`𝒞`$-equivalent (i.e. $`f__t=g__t`$).

Then the classification theory of function germs by the $`𝒞`$-equivalence is quite useful for our purpose. For each function germ $`f:(𝐑,0)(𝐑,0),`$ we set

$`𝒞\mathrm{cod}(f)`$ $`=`$ $`\mathrm{dim}_𝐑_t/f__t,`$
$`𝒞^+\mathrm{cod}(f)`$ $`=`$ $`\mathrm{dim}_𝐑_t/f__t+{\displaystyle \frac{df}{dt}}_𝐑,`$
$`𝒦\mathrm{cod}(f)`$ $`=`$ $`\mathrm{dim}_𝐑_t/f__t+{\displaystyle \frac{df}{dt}}__t.`$
Then we have the following well-known classification (cf. ).

###### Lemma 5.4

Let $`f:(𝐑,0)(𝐑,0)`$ be a function germ with $`𝒦\mathrm{cod}(f)<\mathrm{}.`$ Then $`f`$ is $`𝒞`$-equivalent to the map germ $`t^{\mathrm{}+1}`$ for some $`\mathrm{}.`$

By the direct calculation, we have

$`𝒞\mathrm{cod}(t^{\mathrm{}+1})`$ $`=`$ $`\mathrm{}+1,`$
$`𝒞^+\mathrm{cod}(t^{\mathrm{}+1})`$ $`=`$ $`\mathrm{}.`$
Thus we can easily determine $`𝒞`$ (respectively, $`𝒞^+`$)-versal deformations of the above germs by the usual method as follows :

The $`𝒞`$-versal deformation :

$$t^{\mathrm{}+1}+\underset{i=0}{\overset{\mathrm{}}{}}u_{i+1}t^i.$$

The $`𝒞^+`$-versal deformation :

$$t^{\mathrm{}+1}+\underset{i=0}{\overset{\mathrm{}1}{}}u_{i+1}t^i.$$

We now ready to give a proof of Theorem 2.7.

Proof of Theorem 2.7  Let $`(\mu ,f)`$ be a Clairaut type equation and the corresponding Legendrian unfolding $`\mathrm{}_{(\mu ,f)}`$. By the previous arguments, we have a function germ $`F:(𝐑\times 𝐑,0)(𝐑,0)`$ such that $`\mathrm{Image}j_1^1F=\mathrm{Image}\mathrm{}_{(\mu ,f)}.`$ Therefore we consider generic property of $`F(t,x).`$ By definition $`j_1^1F`$ is an immersion germ if and only if

$$(\frac{F}{t}(0),\frac{^2F}{tx}(0))0.$$

Under this condition, we have the characterization of the fold point and the cusp point of $`\pi j_1^1F`$ as follows (cf., ):

(A) $`\pi j_1^1F`$ is the fold germ if and only if $`{\displaystyle \frac{F}{t}}(0)=0`$ and $`{\displaystyle \frac{^2F}{t^2}}(0)0.`$

(B) $`\pi j_1^1F`$ is the cusp germ if and only if

$$\frac{F}{t}(0)=\frac{^2F}{t^2}(0)=0\text{ and }\frac{^2F}{tx}(0)\frac{^3F}{t^3}(0)0.$$

When $`j_1^1F`$ is not a immersion germ, we have the following characterization of the cross cap:

(C) $`j_1^1F`$ is a cross cap germ if and only if

$$\frac{F}{t}(0)=\frac{^2F}{tx}(0)=0\text{ and }\frac{^3F}{tx^2}(0)\frac{^2F}{t^2}(0)0.$$

In the first place, we give normal forms under the assumption that the conditions (a),(b),(c). Let assume that the condition (C) holds. In this case the function germ has the following form:

$$F(t,x)=at^2+bx^2+ctx^2+h(t,x),$$

where $`a0,c0`$ and $`h(0,0)=0.`$ Since $`F(t,0)=at^2+h(t,0)`$ is $`𝒞`$-equivalent to $`t^2,`$ $`F(t,x)`$ is $`P`$-$`𝒞`$-equivalent to a deformation of $`t^2.`$ By the previous arguments, the $`𝒞`$-versal deformation of $`t^2`$ is $`t^2+v_1t+v_2.`$ Therefore, $`F(t,x)`$ is $`P`$-$`𝒞`$-equivalent to the function germ of the form:

$$G(t,x)=t^2+t\varphi _1(x)+\varphi _2(x).$$

Since $`j_1^1G`$ is also a cross cap germ, we have

$$\varphi _1(x)=\alpha x^2+\mathrm{higher}\mathrm{order}\mathrm{term},$$

with $`\alpha 0.`$ By a local diffeomorphim of the variable $`x,`$ we have $`\varphi _1(x)=x^2.`$ This means that $`F(t,x)`$ is $`P`$-$`𝒞`$-equivalent to the germ of the form $`t^2+tx^2+\varphi (x).`$ Hence, we might put that $`F(t,x)=t^2+tx^2+\varphi (x).`$ In this case.

$$j^1F(t,x)=(t,x,t^2+tx^2+\varphi (x),2t+x^2,2tx+\varphi ^{}(x)).$$

The corresponding integral diagram is

$$\mu (u_1,u_2)=u_1,g(u_1,u_2)=(u_2,u_1^2+u_1u_2^2+\varphi (u_2)).$$

On the $`(x,y)`$-plane, we have a diffeomorphism germ $`\mathrm{\Psi }:(𝐑^2,0)(𝐑^2,0)`$ defined by $`\mathrm{\Psi }(x,y)=(x,{\displaystyle \frac{1}{4}}(y+{\displaystyle \frac{1}{4}}x^4\varphi (x)))`$. Then we have $`\mathrm{\Psi }g(u_1,u_2)=(u_2,{\displaystyle \frac{1}{4}}(u_1+{\displaystyle \frac{1}{2}}u_2^2)^2).`$ This is the normal form (4) in Theorem 2.7, after setting $`(u,v)=(u_2,u_1+{\displaystyle \frac{1}{2}}u_2^2)`$.

For the case (A), we can apply almost the same arguments as the above and get the normal form of (2) in Theorem 2.7.

For the case (B), the situation is a rather different. In this case the function $`F(t,0)`$ is $`𝒞`$-equivalent to $`t^3.`$ The $`𝒞`$-versal deformation of $`t^3`$ is $`t^3+v_1t^2+v_2t+v_3,`$ then the above arguments cannot work in this case. However, the $`𝒞^+`$-versal deformation of $`t^3`$ is $`t^3+v_1t+v_2.`$ Thus we can apply almost the same arguments as the above and the corrersponding integral diagram is $`^+`$-equvialent to

$$\mu (u_1,u_2)=u_2,g(u_1,u_2)=(u_1,u_2^3+u_1u_2).$$

This means that the diagram is strictly equivalent to the normal form (3) in Theorem 2.7.

We now show that the set of function $`F(t,x)`$ satisfying the conditions (A), (B), (C) or (R) at any point are generic in the space of all functions (equiped with the Whitney $`C^{\mathrm{}}`$-topology). Here the condition (R) is that $`{\displaystyle \frac{F}{t}}(0)0.`$ Let $`J^3(2,1)`$ be the set of $`3`$-jets of function germs $`h:(𝐑^2,0)(𝐑,0).`$ We consider the following two algebraic subset of $`J^3(2,1):`$

$$\mathrm{\Sigma }_1=\left\{j^3h(0)\right|\frac{h}{t}(0)=\frac{^2h}{tx}(0)=\frac{^2h}{t^2}(0)\frac{^3h}{tx^2}(0)=0\},$$

$$\mathrm{\Sigma }_2=\left\{j^3h(0)\right|\frac{h}{t}(0)=\frac{^2h}{t^2}(0)=\frac{^2h}{tx}(0)\frac{^3h}{t^3}(0)=0\}.$$

We consider the union $`W=\mathrm{\Sigma }_1\mathrm{\Sigma }_2,`$ then it is also an algebraic subset of $`J^3(2,1).`$ We can stratify the algebraic set $`W`$ by submanifolds whose codimensions are at least $`3.`$ By Thom’s jet transversality theorem, $`j^3F(𝐑^2)(𝐑^2\times 𝐑\times W)=\mathrm{}`$ for a generic function $`F(t,x).`$ We can easily show that the conditions (A),(B),(C) or (R) are satisfied for such a function $`F(t,x).`$ This completes the proof of Theorem 2.7 and in particular we have Therem 2.6. $`\mathrm{}`$

A.A. Davydov; Department of Mathematics, Vladimir State University, Gorkii street 87, 600026 Vladimir, Russia.

e-mail : davydov-m2@vpti.vladimir.su

G. Ishikawa; Department of Mathematics, Hokkaido University, Sapporo 060-0810, Japan.

e-mail : ishikawa@math.sci.hokudai.ac.jp

S. Izumiya; Department of Mathematics, Hokkaido University, Sapporo 060-0810, Japan.

e-mail : izumiya@math.sci.hokudai.ac.jp

W.-Z. Sun; Department of Mathematics, North East Normal University, Chang Chun 130024, P.R. China.

e-mail : wzsun@nenu.edu.cn


format: text

text:
# Quantum Mechanics with Trajectories: Quantum Trajectories and Adaptive Grids

## I Introduction

The concept of a path or trajectory plays a central role in our understanding of the motion of objects and fluids. Much like a route traced on a road map, a trajectory tells us where an object started, where it goes, and how it gets there. There may be alternate routes, some more likely than others. Hence, analysis of a set of trajectories provides us with an intuitive tool for understanding possible complex dynamics. Macroscopic objects obey Newton s equation of motion, $`m\ddot{q}=f(q(t))`$, where $`q`$ is the position of the particle at time $`t`$ and $`f(q(t))`$ is the force acting on it. Given values for both position and velocity at time $`t`$, we can compute the trajectory that the object will follow and as a result, we can predict with certainty where the object will wind up.

However, at the atomic and molecular level where objects obey the rules of quantum mechanics, Newton s equations of motion are no longer strictly valid and the concept of a unique trajectory given a set of initial conditions becomes murky at best. This is because, fundamentally, quantum mechanics is non-local, an issue to which we will return later. In addition, the Heisenberg uncertainty principle states that when measurements are made, we cannot simultaneously determine with infinite precision the exact position and velocity of a quantum particle, although in principle this is possible in Newtonian mechanics. Consequently, it seems as though we cannot speak in terms of the unique trajectory followed by a quantum mechanical object.

However, even in quantum mechanics, Feynman ref1  showed that we can talk in terms of paths, in fact, ensembles of them. For example, if we specify two fixed end-points, $`q(0)`$ and $`q(t)`$, then we can compute the probability of a particle starting at $`q(0)`$ and winding up at $`q(t)`$ by summing over all possible paths connecting the two points (these paths include the classical trajectory linking these points, if there is one) and weighting each path by the complex-valued factor $`\mathrm{exp}(iS/\mathrm{})`$ , where $`S(t)`$ is the classical action integral

$`S(t)={\displaystyle _0^t}{\displaystyle \frac{1}{m}}\dot{q}^2V(q(t))dt`$ (1)
in which the integrand is the classical Lagrangian. What we lose here is any indication of exactly which path the particle actually follows, although we can say that some paths are more likely than others, especially those close to the paths predicted by classical mechanics. In order to make a prediction, it is as if the particle must explore all possible routes between the two end points.

Another way that we can describe the motion of quantum mechanical objects is via the hydrodynamic formulation of quantum mechanics, which was first introduced in the late 1920 s and later explored and extensively developed by Bohm ref2  starting in the early 1950 s. Here we will define an ensemble of quantum trajectories, each with a precisely defined coordinate and velocity that uniquely characterizes the dynamical evolution of a quantum system. However, we still are not free to talk about independent trajectories as in classical mechanics; quantum trajectories are coupled together and evolve as a correlated ensemble. This correlation is a unique feature of quantum dynamics and is expressed as a non-local potential in the hydrodynamic formulation.

In this Article, we will explore the use of the hydrodynamic viewpoint of quantum mechanics to design new computational tools for predicting the evolution of quantum systems. Within the past few years, starting in 1999, it has become possible to directly solve the quantum hydrodynamic equations to predict the space-time dynamics of elements of the probability fluid.ref3 ; ref4 ; ref5  Elements of this quantum fluid are linked through the Bohm quantum potential, denoted $`Q`$, which is computed on-the-fly as the equations of motion are integrated to generate the hydrodynamic fields. The quantum potential introduces all quantum features into the dynamics, including interference effects, barrier tunneling, zero point energies, etc. The only approximation made in solving the hydrodynamic equations involves the use of a relatively small number of fluid elements. The equations of motion for these fluid elements are expressed in the moving with the fluid, Lagrangian, picture of fluid flow.

One implementation of these ideas, referred to as the quantum trajectory method ref3  (QTM), has been used to predict and analyze the dynamics of wavepackets in a number of scattering problems. An approach similar to the QTM is the quantum fluid dynamic method (QFD)ref5 . Within the past two years, there has been a surge of interest in the development and application of trajectory methods for solving both the time-dependent Schrödinger equation and density matrix equations of motion. Novel quantum trajectory methods for evolution of the reduced density matrix in both non-dissipative and dissipative systems have also been developedref6 . In order to circumvent computational problems associated with the propagation of Bohmian trajectories (especially in regions where nodes develop) adaptive grid strategies have been recently explored.ref4 ; ref7 ; ref8

In Sec. II of this article, the Bohmian formulation of the hydrodynamic equations of motion will be reviewed, and computational implementations will be described. Recent studies that employ cluster modeling to calculate the quantum potential will also be presented in this section ref9 . Section III continues with a discussion of adaptive dynamic grid techniques, and the transforms of the hydrodynamical equations appropriate for moving grids are presented. Finally, some suggestions for further study are presented in Sec.IV.

## II Implementation of the quantum hydrodynamic equations

### II.1 The equations of motion

In this Section, the equations needed to implement the quantum hydrodynamic formulation will be reviewed. This formulation is initiated by substituting the amplitude-phase decomposition of the time-dependent wavefunction, $`\psi (y,t)=R(y,t)\mathrm{exp}[iS(y,t)/\mathrm{})`$, into the time-dependent Schrödinger equation,

$`i\mathrm{}{\displaystyle \frac{d\psi }{dt}}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}^2+V\right)\psi `$ (2)
In terms of $`R`$ and $`S`$, the probability density and the local flow velocity are given by $`\rho =|\psi |^2=R^2`$ and $`v=j/\rho `$ , where $`j`$ is the probability current. The Lagrangian form of the hydrodynamic equations of motion resulting from this analysis are given by:

$`{\displaystyle \frac{d\rho }{dt}}=\rho v`$ (3)
$`{\displaystyle \frac{dv}{dt}}={\displaystyle \frac{1}{m}}(V+Q)`$ (4)
in which the derivative on the left side is appropriate for calculating the rate of change in a function along a fluid trajectory. Equation (3) is recognized as the continuity equation and Equation (4) is a Newtonian-type equation in which the flow acceleration is produced by the sum of the classical force, $`f_c=V`$, and the quantum force is $`f_q=Q`$, where $`V`$ is the potential energy function and $`Q`$ is quantum potential ref2 . The quantum potential measures the curvature induced internal stress and is given by

$`Q(y,t)={\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{1}{R(r,y)}}^2R(y,t)={\displaystyle \frac{\mathrm{}^2}{2m}}\rho ^{1/2}^2\rho ^{1/2}`$ (5)
Computation of the quantum potential is frequently rendered more accurate if derivatives are evaluated using the amplitude $`C=\mathrm{log}(R)`$ ($`C`$ is referred to as the $`C`$-amplitude). In terms of derivatives of this amplitude, the quantum potential is $`Q=\mathrm{}^2(^2C+(C)^2)/2m`$.

Closure of the set of dynamical equations is obtained by introducing the equation for the action function,

$`{\displaystyle \frac{dS}{dt}}={\displaystyle \frac{1}{2m}}(S)^2(V+Q)=L_q`$ (6)
This equation, referred to as the quantum Hamilton-Jacobi equation, relates the change in the action to the quantum Lagrangian.

The quantum trajectories obey two important noncrossing rules: (1) they cannot cross nodal surfaces (along which $`\psi =0`$ ); (2) they cannot cross each other. (In practice, because of numerical inaccuracies, these conditions may be violated.) The quantum trajectories are thus very different from the paths used to assemble the Feynman propagator.

No approximations were made in deriving these equations from the time-dependent Schrödinger equation.  However, in order to generate solutions to these equations, an approximation will be made. The initial wavepacket will be subdivided into $`N`$ fluid elements and the equations of motion will be used to find the hydrodynamic fields along the trajectories followed by these elements. As time proceeds, the fluid elements evolve into an unstructured mesh, and this presents problems for accurate derivative evaluation. For this purpose, least squares methods may be used.ref3 ; ref4

From the values of $`C`$ and $`S`$ carried by each evolving fluid elements, the complex-valued wavefunction may be synthesized. An example is shown in Fig. 1: this wavepacket has just propagated downhill on an 11 degree-of-freedom potential surface This complicated, oscillating wavefunction (only the real part of the wavefunction is shown) was synthesized from the information carried along 110 quantum trajectories. It is remarkable that this complicated structure can be built from the information that is propagated along so few quantum trajectories.

### II.2 Multidimensional dynamics with cluster modeling

For some systems, including those of low dimensionality, the methods introduced in the previous section allow us to solve the quantum hydrodynamic equations of motion to high accuracy. For higher dimensional systems, such as atomic clusters or larger molecules, an extremely accurate description of the quantum motion is very difficult to obtain and perhaps more difficult to understand. In this section, we will discuss how the hydrodynamic/Bohmian method can be used to develop an approximate statistical approach for high dimensional systems.ref9

We can specify the configuration of a collection of particles by a set of coordinates, $`R(t)=\{r_1(t),r_2(t)\mathrm{},r_N(t)\}`$. The probability of finding one of these configurations is given by the product of the volume element $`dV`$ weighted by the probability density $`\rho (R(t))`$. In essence, each configuration is an element in an ensemble of possible configurations and the act of measurement pulls one of these configurations out of the hat. In quantum mechanics, the probability density is related to square amplitude of the quantum wave function, $`\psi (R,t)`$, Finally, we apply the Bohmian postulate that the time evolution of a configuration, $`r_n`$, is given by ref2

$`\dot{r}_n={\displaystyle \frac{1}{m}}_nS(R,t)|_{R=r_n}`$ (7)
where $`S`$ is the action function for the quantum wave function as described in the previous section. In this section, we discuss how one can use a Bayesian statistical approach to approximate $`Q`$ by assuming that the density can be cast as a superposition of gaussian product states.

Any statistical distribution can be written as a superposition of a set of $`M`$ clusters ,

$`\rho (R)={\displaystyle \underset{m=1}{\overset{M}{}}}p(R,c_m)`$ (8)
where $`p(R,c_m)`$ is the probability of finding the system in configuration $`R`$ and being in the $`m`$-th cluster. We then use Bayesref10  theorem to break this joint probability into a conditional probability $`p(R|c_m)`$ which tells us the probability that a configuration is a member of $`R`$ knowing that it is also a member of the $`m`$-th cluster and a marginal probability, $`p(c_m)`$ which gives the likelihood of being in the $`mth`$ cluster. We then pick a functional form for the conditional probabilities by writing

$`p(R|c_m)={\displaystyle \frac{𝐜_m^1^{1/2}}{(2\pi )^{3N/2}}}\mathrm{exp}\left[(R\mu _m)^T𝐜_m^1(R\mu _m)\right],`$ (9)
where $`c_m`$ is the $`N\times N`$ dimensional covariance matrix and $`\mu _m`$ is the center of the gaussian in $`N`$ spatial dimensions. To determine the coefficients of the gaussians, we maximize the log-likelihood, $``$, of a given trial set of gaussians actually corresponding to the data by taking the variation $`\delta =0`$. Furthermore, since the density is now represented as a superposition of gaussians, it is a straightforward task to compute the quantum potential, $`Q`$, which is required to integrate the hydrodynamic equations (see Eq. 4). Note that the methodology itself is extremely general and can be used to estimate the probability distribution function (PDF) of any distribution of sample points.

In general, the $`N`$-dimensional covariance matrices specifying each gaussian cluster reflect the correlation between various degrees of freedom. A fully specified covariance will require that we solve $`N^2`$ simultaneous coupled equations per gaussian cluster. In practice, it is easier to use more gaussians with less covariance than few gaussian clusters with a high degree of covariance.

Our initial application of this method has been to determine ground states of high-dimensional systems, such as (He)n clusters for up to $`n=4`$. In Fig. 2, we demonstrate this by showing the convergence of an initial distribution of sampling points to the vibrational ground state for the anharmonic CH3-I bending/stretching modes of the methyl iodide molecule in its lowest electronic state. This system serves as a useful benchmark for our approach since the vibrational ground state can also be computed by diagonalization of the Hamiltonian matrix. Using this approach gives the ground state energy as 591.045 cm-1 above the bottom of the well. The discrete variable representation (DVR) grid used to converge this and the lowest few vibrational excited states is shown superimposed over potential energy surface in Fig.2(a). In the clustering approach for bound-state problems, the sample points remain localized in a small volume of configuration space as the algorithm relaxes to the minimal energy configuration. In Fig.2(d), we show the energy as the system relaxes. The horizontal solid grey line is the “exact” result. In Fig.2(b), our initial sampling is well away from the lowest energy point on the potential energy surface. The superimposed ovals give the location and widths of the $`M`$ gaussian clusters used in this calculation. For this calculation we use 4 fully factorized gaussians with no covariance. As the calculation progresses, the distribution follows the potential energy curve to its minimum and the distribution takes the form of the lowest energy quantum state. The final relaxed ground-state energy is $`665.205\pm 33.6`$ cm-1. At this point, the sample points are still in motion and we are considerably above the DVR ground state energy. Finally, we consider what happens if one takes fewer Gaussians each with full covariance in $`(x,y)`$. Starting from the same distribution as before, we propagate the particles applying a viscous force to bleed away the kinetic energy. Under this approach, the system relaxes to $`580.0\pm 10.1`$ cm-1 using only two Gaussians. Surprisingly, most of the amplitude in the final state is concentrated in the Gaussian located above the minimum in the potential energy surface. Looking at Fig.2(d), we notice that the energy fluctuates about the final average with a relatively small deviation. However, there are some spikes and these correspond to cases in which one of the clusters, typically the smaller one, suddenly jumped to a different position or the correlation matrix rotated slightly. Fortunately, the algorithm rapidly corrects for these excursions.

In many ways, the hydrodynamic formalism implemented with the clustering approach is akin to the guide function Monte Carlo (GFMC) widely used to calculate the ground state of multi-particle systems. For example, in GFMC one chooses an appropriate guide-function from a variational Monte Carlo calculation and uses the guide function to improve the Monte Carlo sampling. Much like the Bohmian quantum potential, the local kinetic energy determined by the curvature in the distribution function and iterative refinements to the sampling are made using the Monte Carlo algorithm. In our case, both the guide function and the sampling points are determined at each step in the calculation and the sampling points evolve according to hydrodynamic equations of motion. We expect that this approach will be most useful in high-dimensional problems where conventional variational basis set approaches are numerically impossible.

## III  Beyond Bohmian mechanics: Adaptive quantum paths

### III.1 What kinds of grids and paths are there?

Near nodes or nodal surfaces, places where , the propagation of Bohmian trajectories becomes problematic. On one hand, as mentioned earlier, these trajectories tend to avoid nodal regions and this leads to an under-sampling problem. For trajectories near the node, derivative evaluation of the hydrodynamic fields becomes increasingly inaccurate, thus leading to computational breakdown. A closely related problem is that the quantum potential near nodes acquires large values, and small errors in evaluating $`Q`$ in turn leads to large errors in the positions and momenta for the nearby quantum trajectories. In order to circumvent these difficulties, we will take a more general look at adaptive dynamic (moving) grids.ref11  The advantage of these grids is the flexibility afforded in directing the grid paths. Furthermore, these dynamically adaptive grids may be robust enough that accurate results can be obtained over long propagation times.

When the spatial coordinates associated with hydrodynamic equations of motion are discretized, each grid point follows a time-dependent path, $`x_j(t)`$, with an associated grid velocity, $`\dot{x}_j`$. The primary advantage of time-dependent grids is that many fewer points may be required because the grid can be chosen to adapt to the evolving hydrodynamic fields. In favorable cases, interesting features may be captured with the moving grid that would require significantly more fixed grid points in order to achieve comparable resolution. There are two sub-categories of moving grids.

* Lagrangian grids. For this frequently used type of dynamic grid, the grid point velocity is the same as the local flow speed of the fluid, . As a consequence, the grid points march in step with the fluid. Quantum trajectories obtained by integrating the hydrodynamic equations of motion fall within this category (Sec. II. A.) These trajectories follow along the fluid flow, but an observer has no control over the paths taken by the grid points.
* Arbitrary Lagrangian Eulerian (ALE) grids. For these intermediate grids, the grid point velocity is not the same as that of the fluid,ref11 . As a result, the grid either advances on the fluid, or the grid may lag the fluid flow. In this case, it is useful to introduce the slip velocity, , whose sign may be either positive or negative. The grid developer has considerable flexibility is creating these designer grids to satisfy certain objectives. This design issue will be taken up in more detail in the next section.

### III.2 Moving path transforms of the hydrodynamic equations

For an observer moving along a path $`x(t)`$, the rate of change of a function $`f`$ is denoted $`df/dt`$ (this is the total time derivative), whereas the time derivative at a fixed point (Eulerian frame) is $`f/t`$. These two time derivatives are related through the equation: $`df/dt=f/t+\dot{x}f`$, where the last term is referred to as the convective contribution. From this equation, the derivative at a space-fixed point is given by $`f/t=df/dt\dot{x}f`$. This equation allows us to transform equations of motion expressed in the Eulerian frame into the more general ALE frames.

When the three hydrodynamic equations given earlier for $`\rho `$, $`v`$, and $`S`$ are transformedref8 , we obtain the new equations (in the first equation, the $`C`$-amplitude is $`C=\mathrm{ln}\rho /2`$).

$`{\displaystyle \frac{dC}{dt}}`$ $`=`$ $`w{\displaystyle \frac{C}{x}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{v}{x}}`$ (10)
$`{\displaystyle \frac{dS}{dt}}`$ $`=`$ $`w(mv)+L_q(t)`$ (11)
$`m{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`mw{\displaystyle \frac{v}{x}}{\displaystyle \frac{}{x}}(V+Q)`$ (12)
where, for simplicity, a one dimensional problem has been assumed. It is important to note that the slip velocity appears in the first term in each of these equations. When $`w=0`$, the condition appropriate for a Lagrangian frame, these equations revert back to the usual ones that were presented earlier. However, in their more general form, these equations can be integrated along paths designed to capture features that develop during the course of the dynamics. How this can be done forms the subject of the next section.

### III.3 Grid adaptation using the monitor function

The local hydrodynamic fields surrounding each moving grid point can be monitored with a user designed function $`M(x)`$.ref11  There is considerable flexibility in crafting the monitor function and in principle different features can be captured at different times. Usually, the monitor captures the gradient and/or the curvature of the hydrodynamic fields, for example, $`M(x)=1+\alpha |u|^2+\beta ^2u`$. A constant (unity in this case) is added on the right side to prevent the monitor from becoming very small in regions where the fields are relatively flat. There are many possibilities and the grid designer must decide what features of the fields require monitoring. Once the monitor is introduced, the next issue is how to adjust the grid points so that they can adapt to these fields. In regions where the monitor takes on large values, we would like for the grid points to come closer together, but in such a way that crossing of paths is avoided. A commonly used way to do this is through use of the equidistribution principle.ref11

The equidistribution principle (EP) is easily stated for a one-dimensional problem. First, let $`\{M_j\}`$ denote the values of the monitor at the grid points $`\{x_j\}`$ , for $`j=1,..,N`$. The end points $`x_1`$ and $`x_N`$ are regarded as fixed during the adaption of the internal points. Finally, let $`M_{j+1/2}`$ denote the average value of the monitor between points $`j`$ and $`j+1`$. The EP then states, $`M_{j1/2}(x_jx_{j1})=M_{j+1/2}(x_{j+1}x_j)=\mathrm{constant}`$. A large value of the monitor then forces the spacing between adjacent points to be small. This equation can also be viewed in terms of the equilibration of a spring system, in which the local monitor functions act as the spring constants. These smart springs sense features in the hydrodynamic fields. The spring analogy has been used in the solution of classical fluid problems, and it has also been used recently in the solution of the quantum hydrodynamic equations. ref8  In order to prevent relatively sudden changes ( jerkiness ) in the grid, it is useful to consider a refined version of the EP that was first introduced by Dorfi and Drury ref12

### III.4 Application of dynamic grid algorithms

A dynamically adaptive grid was used to study the scattering of an initial gaussian wavepacket from a repulsive Eckart potential (whose shape is similar to that of a gaussian). ref7  When Bohmian trajectories are propagated, the calculation breaks down shortly after ripples start to develop in the reflected wavepacket. However, much longer propagation time are obtained through use of adaption techniques. In this application, the grid points were adapted using a monitor function designed to capture the local curvature of the wavefunction. The initial wavepacket has a translational energy $`E_{trans}=4000cm^1`$, the barrier height is $`V_o=6000cm^1`$, and the barrier is centered at $`x_b=6a.u.`$ This wavepacket was propagated using $`N=249`$ grid points between $`x_o=5`$ and $`x_N=25`$. The two edge points were held fixed during the entire propagation sequence (Eulerian frame), but more recent calculations have employed Lagrangian edge points. The internal points were adapted according to the Dorfi and Drury scheme ref12  that was mentioned earlier.

The paths followed by the grid points (only 1/4-th of which are shown) are plotted in Fig.3. At $`t=0`$, it is seen that the grid points are dominantly clustered in the region of large wavefunction density and curvature between $`x=1`$ and 3. At later times, the wavepacket spreads as it moves toward the barrier and this is reflected in the grid paths by the width of the clustered region spreading as the center heads toward the barrier. At around $`t=40`$fs, wavepacket bifurcation is manifest by the grid paths diverging near the barrier region to follow the reflected and transmitted parts of the wavepacket. After about 40fs, these two wavepackets are moving toward the lower right and upper right of the figure, respectively. By design, the grid points tend to congregate in these two regions of higher density.

## IV Future studies

The use of quantum Lagrangian trajectories and more general dynamic adaptive grids have the potential to solve quantum dynamical problems in multiple dimensions. The grid points follow the evolving probability fluid so that numerous points or basis functions are not wasted in regions of little activity. These methods have already yielded insights into a number of problems in chemical physics, including barrier tunneling, electronic nonadiabiatic dynamics, decoherence, and relaxation in dissipative environments.

In order to make these methods robust enough to tackle multi-dimensional dynamics on anharmonic potential surfaces, additional effort needs to be directed toward the following interconnected issues. (1) We have already mentioned that the quantum potential becomes large around wavefunction nodes and quasinodes and that fluid elements following Bohmian trajectories inflate away from these regions. Accurate computation of the quantum potential is very difficult and this in turn can lead to instabilities in the equations of motion. (2) Maintaining long time stability of the solutions is difficult, in part because of the large quantum potentials that arise locally. (3) Derivative evaluation on the unstructured mesh formed by the evolving fluid elements is difficult, although moving least squares or transformation to a structured grid for subsequent derivative evaluation are useful strategies. (4) Near caustics that form when classical trajectories are integrated, the quantum forces acting on the fluid elements are both large and rapidly changing. The hydrodynamic equations are stiff, and special integration techniques (so far untested) are required for stable propagation. All of the unresolved issues mentioned in this paragraph need further investigation.

###### Acknowledgements.

RW and EB both thank the National Science Foundation and the Robert Welch Foundation for financial support. In addition, many discussions with Jeremy Maddox, Kyungsun Na, Keith Hughes, and Corey Trahan are gratefully acknowledged.


text:
# Reentrant transitions in colloidal or dusty plasma bilayers

## Abstract

The phase diagram of crystalline bilayers of particles interacting via a Yukawa potential is calculated for arbitrary screening lengths and particle densities. Staggered rectangular, square, rhombic and triangular structures are found to be stable including a first-order transition between two different rhombic structures. For varied screening length at fixed density, one of these rhombic phases exhibits both a single and even a double reentrant transition. Our predictions can be verified experimentally in strongly confined charged colloidal suspensions or dusty plasma bilayers.

preprint: PREPRINT

Confined systems exhibit structural and dynamical behavior very much different from the corresponding bulk state Binderreview ; Bechingerreview . In particular freezing is strongly affected by the presence of a planar wall. In equilibrium, solidification near walls can occur at thermodynamic conditions where the bulk is still fluid (so called “prefreezing”) vanSwol ; Heni2 . In nonequilibrium, the wall may act as a center of heterogeneous nucleation Palberg1  in order to initiate crystal growth vanBlaaderen . A system confined between two parallel planar walls exhibits various layered crystalline states at low temperature if the plate distance gets comparable to the mean interparticle distance. For hard spheres between hard plates, geometric packing considerations lead to the stability of different crystalline lattices including multiple square and hexagonal layers Murray  as well as buckled Pansu , rhombic Pansu ; Schmidt  and prism superlattices Bechinger . On the other hand, for pure coulombic systems like (classical) electrons in quantum wells Kalman  or trapped ions Raizen , several crystalline bilayer-structures were reported Peeters .

Most of our experimental knowledge of freezing in confining slit-like geometry is based on real-space measurements of mesoscopic model systems such as charged colloidal suspensions between glass plates Murray ; Bechinger  or of multilayers of highly charged dusty plasmas plasma2 . The actual interaction between these mesoscopic “macroions” is neither hard-sphere like nor pure coulombic but is described by an intermediate screened Coulomb or Yukawa pair potential Robbins ; test  due to the screening via additional microions in the system. The screening length can be tailored by changing the microion concentration: for charged colloids, salt ions are conveniently added to the aqueous suspensions; the complex plasma, on the other hand, consists of electrons and impurity ions.

In the present letter, we study the stability of different crystal lattices in bilayers of Yukawa particles as motivated by the experimental model systems. The zero-temperature phase diagram is calculated for arbitrary screening lengths and particle densities note\_peeters . We find a variety of different staggered solid lattices to be stable which are separated by either first or second-order phase transitions. The two known extreme limits of zero or infinite screening length corresponding to hard-spheres Schmidt  and the plasma Peeters ; Goldoni  are recovered. For intermediate screening lengths, the phase behavior is strikingly different from a simple interpolation between these two limits. First, there is a first-order coexistence between two different staggered rhombic lattices differing in their relative shift of the two unit cells. Second, one of these staggered rhombic phases exhibits a novel reentrant effect for fixed density and varied screening length. Depending on the density, the reentrant transition can proceed via a staggered square or a staggered triangular solid including even a double reentrant transition of the rhombic phase. All our theoretical predictions can in principle be verified in real-space studies of confined charged suspensions or dusty plasmas.

In detail, our system consists of two layers containing in total $`N`$ particles in the $`(x,y)`$ plane. The total area density of the two layers is $`\rho =N/A`$ with $`A`$ denoting the system area in the $`(x,y)`$ plane. The distance $`D`$ between the layers in the $`z`$ direction is prescribed by the external potential confining the system. The particles are interacting via the Yukawa pair potential

$$V(r)=V_0\frac{\mathrm{exp}(\kappa r)}{\kappa r},$$
(1)
where $`r`$ is the central separation. The inverse screening length $`\kappa `$ which governs the range of the interaction is given in terms of the micro-ion concentration by Debye-Hückel screening theory. The energy amplitude $`V_0=Z^2\mathrm{exp}(2\kappa R)\kappa /ϵ(1+\kappa R)^2`$ scales with the square of the charges $`Z`$ of the particles of physical hard core radius $`R`$ point\_particle  reduced by the dielectric $`ϵ`$ permittivity of the solvent ($`ϵ=1`$ for the dusty plasma). Typically, $`Z`$ is of the order of $`100100000`$ elementary charges such that $`V(r)`$ at typical interparticle distances can be much larger than the thermal energy $`k_BT`$ at room temperature justifying formally zero-temperature calculations. Then the energy scale is set by $`V_0`$ alone and phase transitions in large bilayer systems are completely determined by two dimensionless parameters, namely the reduced layer density $`\eta =\rho D^2/2`$ and the relative screening strength $`\lambda =\kappa D`$. For zero temperature, the stable state is solid but different crystalline structures of the bilayers are conceivable. As possible candidate structures we assume that the two two-dimensional periodic lattices in the bilayers are the same, have a simple unit cell and are shifted relative to each other in the lateral direction by a displacement vector $`𝐜`$. If the two layers are labeled with $`A`$ and $`B`$, the particle positions in the $`(x,y)`$ plane of the two layers are given by

$`𝐑_A(m,n)`$ $`=`$ $`m𝐚_\mathrm{𝟏}+n𝐚_\mathrm{𝟐},`$
$`𝐑_B(m,n)`$ $`=`$ $`m𝐚_\mathrm{𝟏}+n𝐚_\mathrm{𝟐}+𝐜,`$ (2)
where $`𝐚_\mathrm{𝟏}`$ and $`𝐚_\mathrm{𝟐}`$ are the primitive vectors of the two-dimensional lattice and $`m,n`$ are integers. The total internal energy $`U`$ is obtained by the double lattice sum

$`U`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐑_A𝐑_A^{}}{}}V(|𝐑_A𝐑_A^{}|)+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐑_B𝐑_B^{}}{}}V(|𝐑_B𝐑_B^{}|)`$ (3)
$`+{\displaystyle \underset{𝐑_A,𝐑_B}{}}V([|𝐑_A𝐑_B|^2+D^2]^{1/2}).`$
In the limit $`N\mathrm{}`$, the stable crystalline structure minimizes the total internal energy per particle $`U/N`$.

We have minimized $`U/N`$ with respect to $`𝐚_\mathrm{𝟏}`$, $`𝐚_\mathrm{𝟐}`$ and $`𝐜`$ under the constraint of prescribed density $`\eta `$ for given $`\lambda `$ mapping out the phase diagram in the $`(\eta \lambda )`$ plane. As a result, five typical staggered lattice structures turn out to minimize $`U/N`$ for different $`\eta `$. Adopting the notation developed for plasma bilayers Goldoni , we label them by I, II, III, IV, and V. As summarized in Table 1, phase I is the staggered rectangular crystal with a fixed aspect ratio $`a_2/a_1`$ of $`\sqrt{3}`$; phase II is also staggered rectangular but with a different aspect ratio $`\gamma `$ interpolating continuously between phase I and the staggered-square phase III where $`a_2/a_1=1`$. The staggered rhombic phase IV has two non-orthogonal lattice unit vectors $`(𝐚_1`$ and $`𝐚_2)`$ forming an angle $`\theta `$ and contains a general lateral shift $`𝐜=\alpha (𝐚_\mathrm{𝟏}+𝐚_\mathrm{𝟐})`$ between the two rhombic lattices. In fact, we find two possibilities for $`\alpha `$ defining two variants of stable rhombic phases which we call IVA and IVB. For IVA, $`\alpha =1/2`$ while $`\alpha <1/2`$ for IVB. Finally, phase V is a staggered triangular crystal. Both phases III and V can be considered as special cases of the rhombic phase IV; the former has $`\theta =\pi /2`$ and $`\alpha =1/2`$ while the latter is characterized by $`\theta =\pi /3`$ and $`\alpha =1/3`$.

The result for the phase diagram for a wide range of screening strengths $`(0\lambda 100;\lambda \mathrm{})`$ and densities $`(0\eta 0.8)`$ is shown in Fig. 1.

At very low screening $`\lambda `$, we recover the known plasma limit Goldoni , with our labeling of the phases being in line with their sequence for increasing density $`\eta `$. Phase I has a finite but extremely small density region of stability up to $`\eta =3.6\times 10^5`$ at $`\lambda =0`$ Goldoni . For finite $`\lambda `$, the $`\text{I}\text{II}`$ transition stays second-order and occurs at even smaller densities which decrease monotonically to zero as a function of $`\lambda `$ until the hard-sphere limit $`\lambda \mathrm{}`$ is reached. In this latter case $`a_1`$ is playing the role of an effective particle diameter. This is sketched by the vertical line in Fig. 1. The $`\text{II}\text{III}`$ transition is second-order and the transition densities increase drastically with growing $`\lambda `$ and interpolating monotonically between the plasma and hard sphere limit. More details of the $`\text{I}\text{II}\text{III}`$ transition scenario are depicted in Fig. 2 where the aspect ratio $`\gamma `$ of phase II is shown versus $`\eta `$ for different $`\lambda `$. Phases I and III correspond to $`\gamma =\sqrt{3}`$ and $`\gamma =1`$, respectively. As can be clearly deduced from Fig. 2, the aspect ratio $`\gamma `$ interpolates continuously as a function of $`\eta `$ between $`\sqrt{3}`$ and $`1`$ for any $`\lambda `$ such that both the $`\text{I}\text{II}`$ and the $`\text{II}\text{III}`$ transitions are second-order. In the hard sphere limit, $`\gamma `$ approaches $`\gamma _{hs}=2\eta +\sqrt{4\eta ^2+3}`$ continuously which is also shown in Fig. 2.

Novel effects are observed for the $`\text{III}\text{IV}\text{V}`$ transitions. First, for small $`\lambda `$ the $`\text{III}\text{V}`$ transition proceeds via a IVA phase, the former being second and the latter first-order. For $`\lambda 8`$, however, there is a strong first-order transition directly from III to V with a large density jump as determined by Maxwell’s construction note\_maxwell . For even higher screening, $`\lambda 30`$, the $`\text{III}\text{V}`$ transition happens via the cascade $`\text{III}\text{VIA}\text{IVB}\text{V}`$. The stability range of the IVA phase becomes smaller for increasing $`\lambda `$ shrinking to zero in the hard sphere limit. Details of the $`\text{III}\text{V}`$ transition scenario can be detected via the order parameters $`\mathrm{sin}\theta `$ and $`\alpha `$ of the lattice minimizing the total potential energy at prescribed density $`\eta `$. Plotting $`\mathrm{sin}\theta `$ and $`\alpha `$ versus $`\eta `$ reveals the order of the transitions (see Fig. 3): A cusp, which is found for the $`\text{III}\text{IVA}`$ transformation, implies a second-order transition. All other transitions are first-order as signaled by discontinuous jumps in at least one of these order parameters. The corresponding coexistence density gap is not shown in Fig. 3 but included in Fig. 1(a). Across the $`\text{IVA}\text{IVB}`$ transition the order parameter jump is small yielding a tiny density gap which can not be resolved in Fig. 1(a).

Our most striking result is reentrant behavior of the IVA phase at fixed density upon varying $`\lambda `$ as indicated in Fig. 1 by the vertical arrow. For $`0.5<\eta <0.525`$, there is reentrance of the VIA phase via the III phase. The full sequence over the whole range of $`\lambda `$ is $`\text{IVA}\text{III}\text{IVA}\text{IVB}`$. For $`0.530<\eta <0.536`$ there is even a double reentrant behavior of the VIA phase via the sequence $`\text{IVA}\text{V}\text{IVA}\text{III}\text{IVA}\text{IVB}`$. This rich scenario is due to a subtle interplay of the range of the interaction in conjunction with the different bilayer lattice structures. Finally, at finite temperatures $`T`$, we performed extensive Monte-Carlo computer simulations with 800 particles in a rectangular-shaped box periodically-repeated in $`x`$ and $`y`$ direction and with hard walls of distance $`D`$ in $`z`$-direction allowing fluctuating $`z`$-positions of the particles (“buckling”). For fixed $`\eta =0.533`$, we investigated four states at $`\lambda =0.5`$, 3.0, 8.2, 10 \[see the four diamonds along the arrow in Fig. 1(b)\] for different $`T`$ up to melting. The melting point is detected via a modified Lindemann-type criterion involving differences of mean-square displacements of nearest-neighbors Peeters . We confirm that the reentrant behavior is stable with respect to increasing $`T`$ up to melting.

In summary, we have calculated the full phase diagram for a Yukawa bilayer at zero temperature by lattice sum minimizations. A competition between three length scales, namely the bilayer distance $`D`$, the averaged particle distance $`\rho ^{1/2}`$, and the range $`1/\kappa `$ of the interaction, induces a rich phase behavior which is different from a simple interpolation of the extreme limits of the confined plasma and the hard sphere system. We predict a coexistence of two different rhombic phases at finite screening and a single and double reentrant scenario for one of the rhombic phases for varied “softness” of the interaction. These effects are in principle detectable in real-space experiments of charged colloidal suspensions confined between plates and in layers of dusty plasmas by tuning the screening strength via the microion concentration. The reentrant effect as obtained here in equilibrium should also manifest itself as an interesting fingerprint in nonequilibrium situations. For example, bilayer crystal nucleation and growth could be greatly stimulated via structures which are energetically close to the stable ones Palberg\_review . Soft particle interactions different from the Yukawa type of Eq. (1), as e.g. inverse power potentials where $`V(r)r^n`$, will lead to similar reentrant effects as long as the softness of the potential (e.g. the exponent $`n`$) is varied. Different realizations of soft interactions occur in sterically-stabilized colloids, in spherical block-copolymer micelles and in star polymers and dendrimers, where the softness of the effective interaction can be tuned by the length and grafting density of the polymer chains or the solvent quality Likos\_PRL . Hence the reentrant scenario should also occur in foam films containing polymer bilayers kolaric . Finally, for a general external potential confining the particles to layers, the bilayer distance $`D`$ is not prescribed but the system will minimize its total energy realizing an optimal $`D`$. In this case, second-order phase transitions will still be described in terms of scaled parameters. This implies a universal behavior of our bilayer phase diagram. In a general external potential, however, the system has the additional possibility to split into tri- and higher order multilayers. This can happen either discontinuously or continuous via merging prism phases. Details have to be explored in future studies.

###### Acknowledgements.

We thank C. N. Likos and M. Schmidt for helpful discussions, and the DFG (SFB TR6) for financial support.