Analytical Study of Semi-Discrete Hamiltonian Dynamics in Diverse Potential Landscapes with Discrete Symmetries

1. Introduction

Semi-discretization offers a poIrful approach to modeling physical systems by treating some variables discretely while others remain continuous. Preserving the inherent discrete symmetries of these systems is crucial for accurate and physically relevant numerical methods. This paper presents an analytical study of semi-discrete time-stepping methods applied to a particle in two dimensions, exploring their behavior across five distinct potential functions exhibiting various discrete symmetries. This methodology combines Euler’s method for continuous evolution with probabilistic hopping for discrete spatial coordinates. For each potential, I analyze symmetry preservation, conserved quantities, fixed points, and the characteristics of the probabilistic transitions. This work aims to provide insights into the dynamics of hybrid discretization techniques for Hamiltonian systems with diverse discrete symmetries.

2. Analysis of Semi-Discrete Dynamics in a Potential with Mixed Polynomial and Sinusoidal Contributions

Solving Methodology:

The solving methodology for this potential involves first semi-discretizing the Hamiltonian by treating the x-coordinate as discrete and applying Euler’s method to the continuous variables. The discrete x-coordinate is updated probabilistically based on the x-momentum. Analytical investigations then focus on examining the symmetries of the potential and the semi-discrete map, analyzing the conservation of energy using the derived update rules, identifying fixed points by solving the resulting algebraic equations, and characterizing the expected behavior of the probabilistic stepping process [1,2]

2.1. Semi-Discretization Methodology

The Potential Funciton :

$$V(x,y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

2.1.1. Discretization of the x-Coordinate

For the discretization of the x-coordinate, a common and often convenient approach is to use a uniform grid. You can define the discrete x-coordinate as:

$$x_i=i\Delta x$$

where:

• i is an integer ($...,-2,-1,0,1,2,...$) representing the index of the discrete point.
• $\Delta x$ is the grid spacing, which is a positive constant representing the distance betIen adjacent discrete points.

Implications of this Discretization:

• The particle’s position in the x-direction will only be alloId to take on these discrete values.
• $\Delta x$ becomes a parameter in yThis system that you will need to choose. The choice of $\Delta x$ can influence the accuracy of yThis semi-discrete approximation compared to the original continuous system. A smaller $\Delta x$ typically leads to a finer resolution and potentially better accuracy but might also increase the computational effort if you Ire to simulate this numerically. For analytical work, $\Delta x$ will appear as a parameter in yThis equations.

2.1.2. Continuous Evolution Using Euler’s Method

Let the Hamiltonian of the system be given by:

$$H(x,y,p_x,p_y)={\displaystyle \frac{p_x^2}{2m}}+{\displaystyle \frac{p_y^2}{2m}}+A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

During the continuous evolution step, the x-coordinate is held fixed at a discrete value $x_i=i\Delta x$. Thus, the effective Hamiltonian for this step is:

$$H(x_i,y,p_x,p_y)={\displaystyle \frac{p_x^2}{2m}}+{\displaystyle \frac{p_y^2}{2m}}+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

I use Hamilton’s equations of motion: $\dot{q}={\displaystyle \frac{\partial H}{\partial p}}$ and $\dot{p}=-{\displaystyle \frac{\partial H}{\partial q}}$.

For the y-coordinate:

$$\dot{y}={\displaystyle \frac{\partial H}{\partial p_y}}={\displaystyle \frac{\partial }{\partial p_y}}\left({\displaystyle \frac{p_x^2}{2m}}+{\displaystyle \frac{p_y^2}{2m}}+\dots \right)={\displaystyle \frac{p_y}{m}}$$

$$\dot{p_y}=-{\displaystyle \frac{\partial H}{\partial y}}=-{\displaystyle \frac{\partial }{\partial y}}\left(\dots +D{sin}^2\left(by\right)+E(y^4+F{y}^2)\right)$$

Using the chain rule:

$${\displaystyle \frac{\partial }{\partial y}}{sin}^2\left(by\right)=2sin\left(by\right)cos\left(by\right)·b=bsin\left(2by\right)$$

$${\displaystyle \frac{\partial }{\partial y}}(y^4+F{y}^2)=4y^3+2Fy=2y(2y^2+F)$$

So,

$$\dot{p_y}=-(bDsin\left(2by\right)+2Ey(2y^2+F))$$

For the x-coordinate and its momentum:

$$\dot{p_x}=-{\displaystyle \frac{\partial H}{\partial x}}{|}_{x=x_i}=-{\displaystyle \frac{\partial }{\partial x}}\left(A{sin}^2\left(ax\right)+B(x^4+C{x}^2)\right){|}_{x=x_i}$$

$${\displaystyle \frac{\partial }{\partial x}}{sin}^2\left(ax\right)=2sin\left(ax\right)cos\left(ax\right)·a=asin\left(2ax\right)$$

$${\displaystyle \frac{\partial }{\partial x}}(x^4+C{x}^2)=4x^3+2Cx$$

So,

$$\dot{p_x}=-(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

Now, applying Euler’s method for a time step h:

$$y_{n+1}\approx y_n+h\dot{y}\left(t_n\right)=y_n+{\displaystyle \frac{p_{yn}}{m}}h$$

$$p_{yn+1}\approx p_{yn}+h\dot{p_y}\left(t_n\right)=p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))$$

$$p_{xn+1}\approx p_{xn}+h\dot{p_x}\left(t_n\right)=p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

Here, $x_i$ remains constant during this continuous update from time n to $n+1$.

2.1.3. Probabilistic Update of the Discrete x-Coordinate

The probabilistic update of the discrete x-coordinate i from $i_n$ to $i_{n+1}$ is governed by the value of ${\eta }_n$, defined as:

$${\eta }_n={\displaystyle \frac{p_{xn+1}·h}{m·\Delta x}}$$

where $p_{xn+1}$ is the x-momentum after the continuous Euler update, h is the time step, m is the mass of the particle, and $\Delta x$ is the grid spacing.

The probability of hopping to neighboring sites is then determined as follows:

**Case 1: ${\eta }_n>0$** (Indicates a tendency for the particle to move in the positive x-direction)

The probability of hopping to the next grid point $i_n+1$ is:

$$P(i_{n+1}=i_n+1)=min(1,|{\eta }_n\left|\right)=min\left(1,{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right)$$

The probability of staying at the current grid point $i_n$ is:

$$P(i_{n+1}=i_n)=1-P(i_{n+1}=i_n+1)=1-min\left(1,{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right)$$

The probability of hopping to the previous grid point $i_n-1$ is:

$$P(i_{n+1}=i_n-1)=0$$

**Case 2: ${\eta }_n<0$** (Indicates a tendency for the particle to move in the negative x-direction)

The probability of hopping to the next grid point $i_n+1$ is:

$$P(i_{n+1}=i_n+1)=0$$

The probability of staying at the current grid point $i_n$ is:

$$P(i_{n+1}=i_n)=1-P(i_{n+1}=i_n-1)=1-min\left(1,\left|{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right|\right)=1-min\left(1,-{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right)$$

The probability of hopping to the previous grid point $i_n-1$ is:

$$P(i_{n+1}=i_n-1)=min(1,|{\eta }_n\left|\right)=min\left(1,\left|{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right|\right)=min\left(1,-{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}\right)$$

**Case 3: ${\eta }_n=0$** (Indicates no tendency for the particle to move in either x-direction based on momentum)

The probability of hopping to the next grid point $i_n+1$ is:

$$P(i_{n+1}=i_n+1)=0$$

The probability of staying at the current grid point $i_n$ is:

$$P(i_{n+1}=i_n)=1$$

The probability of hopping to the previous grid point $i_n-1$ is:

$$P(i_{n+1}=i_n-1)=0$$

This probabilistic update completes one time step of the semi-discretization methodology for the first potential function.

2.2. Analysis

2.2.1. Symmetry Analysis

I start with the reflection symmetry about the y-axis, where I consider the transformation $x\to -x$. In the context of This semi-discretization, this corresponds to the discrete index $i\to -i$, so $x_i=i\Delta x$ becomes $x_{-i}=-i\Delta x=-(i\Delta x)=-x_i$.

The first step is to see how the potential function behaves under this transformation. Can you substitute $-x$ for x in the expression for $V(x,y)$ and see if $V(-x,y)$ is equal to $V(x,y)$?

The potential function is:

$$V(x,y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

Substitute $-x$ for x:

$$V(-x,y)=A{sin}^2\left(a(-x)\right)+B({(-x)}^4+C{(-x)}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

To show that the potential $V(x,y)$ is symmetric under the reflection $x\to -x$, I substitute $-x$ into the potential function:

$$V(-x,y)=A{sin}^2\left(a(-x)\right)+B({(-x)}^4+C{(-x)}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

Using the properties that $sin(-z)=-sin\left(z\right)$ and thus ${sin}^2(-z)={(-sin\left(z\right))}^2={sin}^2\left(z\right)$, and that ${(-x)}^4=x^4$ and ${(-x)}^2=x^2$, I get:

$$V(-x,y)=A{(sin\left(ax\right))}^2+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

$$V(-x,y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

Thus, I can see that $V(-x,y)=V(x,y)$, which confirms that the potential is symmetric under reflection about the y-axis.

Next I analyze how the probabilistic update of the discrete x-coordinate behaves under the transformations $i_n\to -i_n$ and $p_{xn+1}\to -p_{xn+1}$. Recall that ${\eta }_n={\displaystyle \frac{p_{xn+1}·h}{m·\Delta x}}$. Under these transformations, ${\eta }_n$ becomes ${\eta }_n^{\prime }={\displaystyle \frac{(-p_{xn+1})·h}{m·\Delta x}}=-{\eta }_n$.

Now, I consider the probabilities:

If, before transformation, ${\eta }_n>0$:

$$P(i_{n+1}=i_n+1)=min(1,{\eta }_n)$$

$$P(i_{n+1}=i_n-1)=0$$

After transformation, ${\eta }_n^{\prime }<0$, and the index $i_n$ becomes $-i_n$. A hop to $i_n+1$ corresponds to a hop to $-i_n-1$, and a hop to $i_n-1$ corresponds to a hop to $-i_n+1$. For ${\eta }_n^{\prime }<0$:

$$P(i_{n+1}=-i_n+1)=0$$

$$P(i_{n+1}=-i_n-1)=min(1,|-{\eta }_n\left|\right)=min(1,{\eta }_n)$$

I see that the probability of hopping from $i_n$ to $i_n+1$ is the same as the probability of hopping from $-i_n$ to $-i_n-1=-(i_n+1)$, which respects the reflection.

Now consider if, before transformation, ${\eta }_n<0$:

$$P(i_{n+1}=i_n+1)=0$$

$$P(i_{n+1}=i_n-1)=min(1,|{\eta }_n\left|\right)=min(1,-{\eta }_n)$$

After transformation, ${\eta }_n^{\prime }>0$:

$$P(i_{n+1}=-i_n+1)=min(1,|-{\eta }_n\left|\right)=min(1,-{\eta }_n)$$

$$P(i_{n+1}=-i_n-1)=0$$

Here, the probability of hopping from $i_n$ to $i_n-1$ is the same as the probability of hopping from $-i_n$ to $-i_n+1=-(i_n-1)$, again respecting the reflection.

Finally, if ${\eta }_n=0$, then ${\eta }_n^{\prime }=0$, and the particle stays at the current index in both cases ($i_{n+1}=i_n$ and $i_{n+1}=-i_n$).

I have shown that the potential function $V(x,y)$ is symmetric under the reflection $x\to -x$. In the semi-discretized system, this corresponds to the transformation of the discrete index $i\to -i$ and the momentum $p_x\to -p_x$. This analysis of the probabilistic update rule for the x-coordinate revealed that the probabilities of hopping to reflected states are equal, thus the semi-discretization respects the reflection symmetry about the y-axis in a statistical sense.

To show that the potential $V(x,y)$ is symmetric under the reflection $y\to -y$, I substitute $-y$ into the potential function:

$$V(x,-y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(b(-y)\right)+E({(-y)}^4+F{(-y)}^2)$$

Using the properties that $sin(-z)=-sin\left(z\right)$ and thus ${sin}^2(-z)={(-sin\left(z\right))}^2={sin}^2\left(z\right)$, and that ${(-y)}^4=y^4$ and ${(-y)}^2=y^2$, I get:

$$V(x,-y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{(sin\left(by\right))}^2+E(y^4+F{y}^2)$$

$$V(x,-y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

Thus, I can see that $V(x,-y)=V(x,y)$, which confirms that the potential is symmetric under reflection about the x-axis.

Next I look at the continuous evolution with the transformed variables $y\to -y$ and $p_y\to -p_y$, starting from a state $(i_n,-y_n,p_{xn},-p_{yn})$ at time n.

Using the Euler update rules:

For the y-coordinate:

$$y_{n+1}^{\prime }=(-y_n)+{\displaystyle \frac{(-p_{yn})}{m}}h=-(y_n+{\displaystyle \frac{p_{yn}}{m}}h)=-y_{n+1}$$

For the y-momentum:

$$p_{yn+1}^{\prime }=(-p_{yn})-h(bDsin\left(2b(-y_n)\right)+2E(-y_n)(2{(-y_n)}^2+F))$$

Since $sin(-z)=-sin\left(z\right)$ and ${(-y)}^2=y^2$, this becomes:

$$p_{yn+1}^{\prime }=-p_{yn}-h(bD(-sin\left(2b{y}_n\right))-2E{y}_n(2y_n^2+F))$$

$$p_{yn+1}^{\prime }=-p_{yn}+hbDsin\left(2b{y}_n\right)+2hE{y}_n(2y_n^2+F)$$

Now I look at $-p_{yn+1}$:

$$-p_{yn+1}=-(p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F)))$$

$$-p_{yn+1}=-p_{yn}+hbDsin\left(2b{y}_n\right)+2hE{y}_n(2y_n^2+F)$$

I see that $p_{yn+1}^{\prime }=-p_{yn+1}$.

For the x-momentum, the transformation $y\to -y$ does not directly affect the update rule, so $p_{xn+1}^{\prime }=p_{xn+1}$.

To verify the symmetry of the continuous evolution under the reflection $y\to -y$ and $p_y\to -p_y$, I recall the Euler update rules:

$$y_{n+1}=y_n+{\displaystyle \frac{p_{yn}}{m}}h$$

$$p_{yn+1}=p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))$$

$$p_{xn+1}=p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

Now, consider a state at time n transformed as $(i_n,-y_n,p_{xn},-p_{yn})$. Applying the update rules to this transformed state yields:

$$y_{n+1}^{\prime }=(-y_n)+{\displaystyle \frac{(-p_{yn})}{m}}h=-(y_n+{\displaystyle \frac{p_{yn}}{m}}h)=-y_{n+1}$$

$$p_{yn+1}^{\prime }=(-p_{yn})-h(bDsin\left(2b(-y_n)\right)+2E(-y_n)(2{(-y_n)}^2+F))$$

Using $sin(-z)=-sin\left(z\right)$, I have:

$$p_{yn+1}^{\prime }=-p_{yn}-h(bD(-sin\left(2b{y}_n\right))-2E{y}_n(2y_n^2+F))$$

$$p_{yn+1}^{\prime }=-p_{yn}+hbDsin\left(2b{y}_n\right)+2hE{y}_n(2y_n^2+F)$$

I also found that $-p_{yn+1}$ is:

$$-p_{yn+1}=-(p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F)))=-p_{yn}+hbDsin\left(2b{y}_n\right)+2hE{y}_n(2y_n^2+F)$$

Thus, $p_{yn+1}^{\prime }=-p_{yn+1}$. The update rule for $p_{xn+1}$ remains unchanged as it does not depend on y or $p_y$. These results show that the continuous evolution using Euler’s method respects the reflection symmetry about the x-axis.

I have established that the potential and the continuous evolution are symmetric under reflection about the x-axis. Now, I consider the probabilistic update of the x-coordinate.

The probability of hopping in the x-direction depends on the parameter ${\eta }_n$, defined as:

$${\eta }_n={\displaystyle \frac{p_{xn+1}·h}{m·\Delta x}}$$

I need to examine how $p_{xn+1}$ is affected by the transformation $y\to -y$ and $p_y\to -p_y$. From the Euler update rule for $p_x$, I have:

$$p_{xn+1}=p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

I can observe that this equation for $p_{xn+1}$ depends on:

• $p_{xn}$: the x-momentum at the current time step.
• h: the time step size.
• $a,A,B,C$: parameters of the potential function related to the x-component.
• $x_i$: the discrete x-coordinate, which is not transformed under $y\to -y$.

The transformation $y\to -y$ and $p_y\to -p_y$ does not appear in the update rule for $p_{xn+1}$. Therefore, the value of $p_{xn+1}$ remains unchanged under this transformation.

Since ${\eta }_n$ depends directly on $p_{xn+1}$ and other constant parameters ($h,m,\Delta x$), ${\eta }_n$ will also remain unchanged. As the probabilities of hopping in the x-direction are determined solely by ${\eta }_n$, these probabilities will also be unaffected by the transformation $y\to -y$ and $p_y\to -p_y$.

Therefore, the probabilistic update for the x-coordinate respects the reflection symmetry about the x-axis.

To analyze the translational symmetry in x, I substitute $x+a_x$ for x in the potential function:

$$V(x+a_x,y)=A{sin}^2\left(a(x+a_x)\right)+B({(x+a_x)}^4+C{(x+a_x)}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

Expanding the terms involving x:

$${sin}^2\left(a(x+a_x)\right)={sin}^2(ax+a{a}_x)$$

$${(x+a_x)}^4=x^4+4x^3{a}_x+6x^2{a}_x^2+4x{a}_x^3+a_x^4$$

$${(x+a_x)}^2=x^2+2x{a}_x+a_x^2$$

Substituting these back into the potential:

$$V(x+a_x,y)=A{sin}^2(ax+a{a}_x)+B(x^4+4x^3{a}_x+6x^2{a}_x^2+4x{a}_x^3+a_x^4+C{x}^2+2Cx{a}_x+C{a}_x^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$$

For $V(x+a_x,y)$ to be equal to $V(x,y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$ for any arbitrary $a_x$, the following conditions would need to hold:

$$A{sin}^2(ax+a{a}_x)=A{sin}^2\left(ax\right)$$

$$B(x^4+4x^3{a}_x+6x^2{a}_x^2+4x{a}_x^3+a_x^4+C{x}^2+2Cx{a}_x+C{a}_x^2)=B(x^4+C{x}^2)$$

The first condition requires $a{a}_x=n\pi$ for some integer n, which means $a_x$ cannot be arbitrary unless $A=0$. The second condition requires all the terms involving $a_x$ to be zero for any x, which is only true if $B=0$ and $a_x=0$.

Therefore, in general, $V(x+a_x,y)\ne V(x,y)$ for an arbitrary constant $a_x$. This means the potential does not exhibit translational symmetry in the x-direction.

To analyze the translational symmetry in y, I substitute $y+a_y$ for y in the potential function:

$$V(x,y+a_y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(b(y+a_y)\right)+E({(y+a_y)}^4+F{(y+a_y)}^2)$$

Expanding the terms involving y:

$${sin}^2\left(b(y+a_y)\right)={sin}^2(by+b{a}_y)$$

$${(y+a_y)}^4=y^4+4y^3{a}_y+6y^2{a}_y^2+4y{a}_y^3+a_y^4$$

$${(y+a_y)}^2=y^2+2y{a}_y+a_y^2$$

Substituting these back into the potential:

$$V(x,y+a_y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2(by+b{a}_y)+E(y^4+4y^3{a}_y+6y^2{a}_y^2+4y{a}_y^3+a_y^4+F{y}^2+2Fy{a}_y+F{a}_y^2)$$

For $V(x,y+a_y)$ to be equal to $V(x,y)=A{sin}^2\left(ax\right)+B(x^4+C{x}^2)+D{sin}^2\left(by\right)+E(y^4+F{y}^2)$ for any arbitrary $a_y$, the following conditions would need to hold:

$$D{sin}^2(by+b{a}_y)=D{sin}^2\left(by\right)$$

$$E(y^4+4y^3{a}_y+6y^2{a}_y^2+4y{a}_y^3+a_y^4+F{y}^2+2Fy{a}_y+F{a}_y^2)=E(y^4+F{y}^2)$$

Similar to the x-translation case, these conditions are generally not met for an arbitrary $a_y$ unless $D=0$ and $E=0$.

Therefore, in general, $V(x,y+a_y)\ne V(x,y)$ for an arbitrary constant $a_y$. This means the potential does not exhibit translational symmetry in the y-direction.

2.2.2. Conservation of Energy

I analyze the change in the continuous part of the Hamiltonian $H_c$ over one time step. At time n, I have:

$$H_c\left(n\right)={\displaystyle \frac{p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn}^2}{2m}}+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(b{y}_n\right)+E(y_n^4+F{y}_n^2)$$

At time $n+1$, using the Euler update rules, I have:

$$y_{n+1}=y_n+{\displaystyle \frac{p_{yn}}{m}}h$$

$$p_{yn+1}=p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))$$

$$p_{xn+1}=p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

The Hamiltonian at time $n+1$ is:

$$H_c(n+1)={\displaystyle \frac{p_{xn+1}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2}{2m}}+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)$$

Substituting the update rules:

$$H_c(n+1)={\displaystyle \frac{1}{2m}}{\left(p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))\right)}^2$$

$$+{\displaystyle \frac{1}{2m}}{\left(p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))\right)}^2$$

$$+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(b\left(y_n+{\displaystyle \frac{p_{yn}}{m}}h\right)\right)+E\left({\left(y_n+{\displaystyle \frac{p_{yn}}{m}}h\right)}^4+F{\left(y_n+{\displaystyle \frac{p_{yn}}{m}}h\right)}^2\right)$$

Expanding the squared terms:

$${\displaystyle \frac{p_{xn+1}^2}{2m}}={\displaystyle \frac{p_{xn}^2-2p_{xn}h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))+h^2{(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))}^2}{2m}}$$

$${\displaystyle \frac{p_{yn+1}^2}{2m}}={\displaystyle \frac{p_{yn}^2-2p_{yn}h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))+h^2{(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))}^2}{2m}}$$

By comparing $H_c(n+1)$ with $H_c\left(n\right)$, it is clear that due to the terms proportional to h and $h^2$ in the kinetic energy part, and also due to the non-linearity of the potential terms in y, $H_c(n+1)\ne H_c\left(n\right)$ in general. The change in energy over one time step is of the order of h. Therefore, Euler’s method is not energy-conserving for this system.

The expected potential energy at time $n+1$, denoted as $\langle V_{n+1}\rangle$, given the particle was at discrete position i at time n, can be expressed as a Iighted sum of the potential energies at the possible next positions ($i+1,i-1,i$), with the Iights being the corresponding probabilities of hopping. Let $y_{n+1}$ be the y-coordinate after the continuous Euler update. Then:

$$\langle V_{n+1}\rangle =P(i_{n+1}=i+1|i_n=i)V(x_{i+1},y_{n+1})+P(i_{n+1}=i-1|i_n=i)V(x_{i-1},y_{n+1})+P(i_{n+1}=i|i_n=i)V(x_i,y_{n+1})$$

I have the following cases for the hopping probabilities based on ${\eta }_n={\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$:

**Case 1: ${\eta }_n>0$**

$$\langle V_{n+1}\rangle =min(1,{\eta }_n)V(x_{i+1},y_{n+1})+0·V(x_{i-1},y_{n+1})+(1-min(1,{\eta }_n))V(x_i,y_{n+1})$$

**Case 2: ${\eta }_n<0$**

$$\langle V_{n+1}\rangle =0·V(x_{i+1},y_{n+1})+min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})$$

**Case 3: ${\eta }_n=0$**

$$\langle V_{n+1}\rangle =0·V(x_{i+1},y_{n+1})+0·V(x_{i-1},y_{n+1})+1·V(x_i,y_{n+1})=V(x_i,y_{n+1})$$

where $V(x_j,y_{n+1})$ for $j\in \{i-1,i,i+1\}$ is given by:

$$V(x_j,y_{n+1})=A{sin}^2\left(a{x}_j\right)+B(x_j^4+C{x}_j^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)$$

To analyze the conservation of average total energy, I would then need to look at $\langle H_{n+1}\rangle ={\displaystyle \frac{p_{xn+1}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2}{2m}}+\langle V_{n+1}\rangle$ and compare it with the energy at time n, potentially also in an average sense if the initial x-position was considered probabilistic. Since the continuous part of the energy is generally not conserved by Euler’s method, and the potential energy changes due to the hop, the average energy is also likely not conserved in general. A more detailed analysis would involve substituting the specific form of the potential and the Euler update rules.

I delve into the detailed algebra for the conservation of average total energy, starting with the case where ${\eta }_n>0$.

The Hamiltonian at time n is:

$$H_n={\displaystyle \frac{p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn}^2}{2m}}+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(b{y}_n\right)+E(y_n^4+F{y}_n^2)$$

The average Hamiltonian at time $n+1$ for ${\eta }_n>0$ is:

$$\langle H_{n+1}\rangle ={\displaystyle \frac{p_{xn+1}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2}{2m}}+min(1,{\eta }_n)V(x_{i+1},y_{n+1})+(1-min(1,{\eta }_n))V(x_i,y_{n+1})$$

Substitute the Euler update for $p_{xn+1}$:

$$p_{xn+1}^2={\left(p_{xn}-h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))\right)}^2=p_{xn}^2-2p_{xn}h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))+O\left(h^2\right)$$

Substitute the Euler update for $p_{yn+1}$:

$$p_{yn+1}^2={\left(p_{yn}-h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))\right)}^2=p_{yn}^2-2p_{yn}h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))+O\left(h^2\right)$$

Substitute the Euler update for $y_{n+1}$ in the potential terms. This will lead to terms involving $y_n$ and poIrs of h. For instance:

$${sin}^2\left(b{y}_{n+1}\right)={sin}^2\left(b(y_n+{\displaystyle \frac{p_{yn}}{m}}h)\right)={sin}^2(b{y}_n+{\displaystyle \frac{b{p}_{yn}}{m}}h)={sin}^2\left(b{y}_n\right)+O\left(h\right)$$

$$y_{n+1}^4={\left(y_n+{\displaystyle \frac{p_{yn}}{m}}h\right)}^4=y_n^4+O\left(h\right)$$

and similarly for other terms in the potential.

Now, I look at $\langle H_{n+1}\rangle -H_n$:

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& =\left[{\displaystyle \frac{p_{xn}^2-2p_{xn}h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))}{2m}}-{\displaystyle \frac{p_{xn}^2}{2m}}\right]\hfill \\ & +\left[{\displaystyle \frac{p_{yn}^2-2p_{yn}h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))}{2m}}-{\displaystyle \frac{p_{yn}^2}{2m}}\right]\hfill \\ & +min(1,{\eta }_n)V(x_{i+1},y_{n+1})+(1-min(1,{\eta }_n))V(x_i,y_{n+1})-V(x_i,y_n)+O\left(h^2\right)\hfill \end{array}$$

I can see terms like:

$$-{\displaystyle \frac{p_{xn}h}{m}}(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))$$

$$-{\displaystyle \frac{p_{yn}h}{m}}(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))$$

and terms arising from the difference in potential energy due to the hop in x and the evolution of y to $y_{n+1}$.

Consider the difference in potential energy terms:

$$\begin{array}{cc}& min(1,{\eta }_n)[A{sin}^2(a(i+1)\Delta x)+B({((i+1)\Delta x)}^4+C{((i+1)\Delta x)}^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)]\hfill \\ & +(1-min(1,{\eta }_n))[A{sin}^2(ai\Delta x)+B({(i\Delta x)}^4+C{(i\Delta x)}^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)]\hfill \\ & -[A{sin}^2(ai\Delta x)+B({(i\Delta x)}^4+C{(i\Delta x)}^2)+D{sin}^2\left(b{y}_n\right)+E(y_n^4+F{y}_n^2)]\hfill \end{array}$$

Even if $y_{n+1}\approx y_n$ for small h, the terms involving $x_{i+1}=x_i+\Delta x$ will generally not cancel out with the terms involving $x_i$, unless the potential is translationally invariant in x, which I have shown is not the case.

Therefore, $\langle H_{n+1}\rangle -H_n\ne 0$ in general. The average energy is also not conserved by this semi-discretization method. A similar analysis holds for ${\eta }_n<0$.

I now consider the case where ${\eta }_n<0$. Recall that $|{\eta }_n|=-{\eta }_n=-{\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$.

The Hamiltonian at time n remains the same:

$$H_n={\displaystyle \frac{p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn}^2}{2m}}+A{sin}^2\left(a{x}_i\right)+B(x_i^4+C{x}_i^2)+D{sin}^2\left(b{y}_n\right)+E(y_n^4+F{y}_n^2)$$

The average Hamiltonian at time $n+1$ for ${\eta }_n<0$ is:

$$\langle H_{n+1}\rangle ={\displaystyle \frac{p_{xn+1}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2}{2m}}+min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})$$

Substituting the Euler update terms for $p_{xn+1}^2$ and $p_{yn+1}^2$ as I did before:

$$\langle H_{n+1}\rangle ={\displaystyle \frac{p_{xn}^2-2p_{xn}h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))+O\left(h^2\right)}{2m}}$$

$$+{\displaystyle \frac{p_{yn}^2-2p_{yn}h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))+O\left(h^2\right)}{2m}}$$

$$+min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})$$

Now consider $\langle H_{n+1}\rangle -H_n$:

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& =\left[{\displaystyle \frac{p_{xn}^2-2p_{xn}h(aAsin\left(2a{x}_i\right)+B(4x_i^3+2C{x}_i))}{2m}}-{\displaystyle \frac{p_{xn}^2}{2m}}\right]\hfill \\ & +\left[{\displaystyle \frac{p_{yn}^2-2p_{yn}h(bDsin\left(2b{y}_n\right)+2E{y}_n(2y_n^2+F))}{2m}}-{\displaystyle \frac{p_{yn}^2}{2m}}\right]\hfill \\ & +min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})-V(x_i,y_n)+O\left(h^2\right)\hfill \end{array}$$

Again, I have terms proportional to h arising from the kinetic energy updates. The difference in potential energy is now:

$$\begin{array}{cc}& min(1,|{\eta }_n\left|\right)[A{sin}^2(a(i-1)\Delta x)+B({((i-1)\Delta x)}^4+C{((i-1)\Delta x)}^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)]\hfill \\ & +(1-min(1,|{\eta }_n\left|\right))[A{sin}^2(ai\Delta x)+B({(i\Delta x)}^4+C{(i\Delta x)}^2)+D{sin}^2\left(b{y}_{n+1}\right)+E(y_{n+1}^4+F{y}_{n+1}^2)]\hfill \\ & -[A{sin}^2(ai\Delta x)+B({(i\Delta x)}^4+C{(i\Delta x)}^2)+D{sin}^2\left(b{y}_n\right)+E(y_n^4+F{y}_n^2)]\hfill \end{array}$$

Similar to the ${\eta }_n>0$ case, the terms involving $x_{i-1}=x_i-\Delta x$ will generally not cancel out with the terms involving $x_i$, and the evolution of y also contributes to energy change. Thus, even for ${\eta }_n<0$, $\langle H_{n+1}\rangle -H_n\ne 0$ in general, and the average energy is not conserved.

2.2.3. Properties of the Probabilistic Step

The probability of transitioning from a discrete x-position i at time n to a position j at time $n+1$ is given by:

$$P(i_{n+1}=j|i_n=i)=\left\{\begin{array}{cc}min(1,{\eta }_n)\hfill & \mathrm{if}j=i+1\mathrm{and}{\eta }_n>0\hfill \\ 0\hfill & \mathrm{if}j=i+1\mathrm{and}{\eta }_n\le 0\hfill \\ min(1,|{\eta }_n\left|\right)\hfill & \mathrm{if}j=i-1\mathrm{and}{\eta }_n<0\hfill \\ 0\hfill & \mathrm{if}j=i-1\mathrm{and}{\eta }_n\ge 0\hfill \\ 1-min(1,|{\eta }_n\left|\right)\hfill & \mathrm{if}j=i\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

where ${\eta }_n={\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$.

The expected change in the discrete x-coordinate over one time step, $\langle \Delta i\rangle =\langle i_{n+1}-i_n\rangle$, is $\langle \Delta i\rangle =\mathrm{sgn}\left({\eta }_n\right)min(1,|{\eta }_n\left|\right)$.

The variance in the x-position after one time step is $min(1,|{\eta }_n\left|\right)(1-min(1,|{\eta }_n\left|\right))$ if ${\eta }_n\ne 0$, and 0 if ${\eta }_n=0$.

The probabilistic step can be vieId as a biased random walk on a 1D lattice. The bias is determined by ${\eta }_n$.

Conclusion:

This investigation of the first potential function using the semi-discretization method revealed that while the potential exhibits reflection symmetry about both the x and y axes, it lacks translational symmetry. The standard Euler method used for the continuous evolution did not conserve energy, and due to the probabilistic hopping in the x-coordinate, the average energy of the system was also found to be generally non-conserved. The probabilistic step itself introduces a biased random walk on the discrete x-lattice, driven by the x-momentum, which allows the system to explore different spatial regions.

3. Analysis of a Semi-Discrete System in a Potential with Exponentially Damped Oscillations

Solving Methodology:

For the potential with coupled trigonometric dependencies, the semi-discretization methodology remains the same. The analytical investigation, hoIver, requires careful handling of the coupled terms when examining symmetries and deriving the change in energy. Finding fixed points involves solving a system of equations that includes these coupled terms, and the probabilistic step analysis remains consistent with the previous cases [3]

The Potential Function :

$$V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

3.1. Semi-Discretization Methodology

3.1.1. Discretization of the x-Coordinate

The semi-discretization involves considering the x-coordinate at discrete points $x_i=i\Delta x$. The Hamiltonian becomes:

$$H(x_i,y,p_x,p_y)={\displaystyle \frac{p_x^2}{2m}}+{\displaystyle \frac{p_y^2}{2m}}+e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y^2}cos\left(by\right)$$

The partial derivatives calculated are:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial H}{\partial p_y}}& ={\displaystyle \frac{p_y}{m}}\hfill \\ \hfill -{\displaystyle \frac{\partial H}{\partial y}}& =2\beta y{e}^{-\beta y^2}cos\left(by\right)+b{e}^{-\beta y^2}sin\left(by\right)\hfill \\ \hfill {\displaystyle \frac{\partial H}{\partial p_x}}& ={\displaystyle \frac{p_x}{m}}\hfill \\ \hfill -{\displaystyle \frac{\partial H}{\partial x}}& =2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\hfill \end{array}$$

3.1.2. Continuous Evolution Using Euler’s Method

$$\begin{array}{cc}\hfill y_{n+1}& =y_n+{\displaystyle \frac{p_{yn}}{m}}h\hfill \\ \hfill p_{yn+1}& =p_{yn}+\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h\hfill \\ \hfill p_{xn+1}& =p_{xn}+\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h\hfill \end{array}$$

where ${\eta }_n={\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$.

The probabilistic update for i is as follows: If ${\eta }_n>0$, then $P(i_{n+1}=i_n+1)=min(1,{\eta }_n)$, $P(i_{n+1}=i_n-1)=0$, and $P(i_{n+1}=i_n)=1-min(1,{\eta }_n)$. If ${\eta }_n<0$, then $P(i_{n+1}=i_n+1)=0$, $P(i_{n+1}=i_n-1)=min(1,|{\eta }_n\left|\right)$, and $P(i_{n+1}=i_n)=1-min(1,|{\eta }_n\left|\right)$. If ${\eta }_n=0$, then $P(i_{n+1}=i_n+1)=0$, $P(i_{n+1}=i_n-1)=0$, and $P(i_{n+1}=i_n)=1$.

3.1.3. Probabilistic Update of the Discrete x-coordinate

The probabilistic update of the discrete x-coordinate i (and hence $x_i$) is based on ${\eta }_n={\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$:

If ${\eta }_n>0$:

$$P(i_{n+1}=i_n+1)=min(1,{\eta }_n)$$

$$P(i_{n+1}=i_n-1)=0$$

$$P(i_{n+1}=i_n)=1-min(1,{\eta }_n)$$

If ${\eta }_n<0$:

$$P(i_{n+1}=i_n+1)=0$$

$$P(i_{n+1}=i_n-1)=min(1,|{\eta }_n\left|\right)$$

$$P(i_{n+1}=i_n)=1-min(1,|{\eta }_n\left|\right)$$

If ${\eta }_n=0$:

$$P(i_{n+1}=i_n+1)=0$$

$$P(i_{n+1}=i_n-1)=0$$

$$P(i_{n+1}=i_n)=1$$

3.2. Analysis

3.2.1. Symmetry Analysis

Let the potential function be

$$V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Now, I substitute x with $-x$:

$$V(-x,y)=e^{-\alpha {(-x)}^2}sin\left(a(-x)\right)+e^{-\beta y^2}cos\left(by\right)$$

I know that ${(-x)}^2=x^2$ and $sin(-ax)=-sin\left(ax\right)$. Therefore,

$$V(-x,y)=e^{-\alpha x^2}(-sin\left(ax\right))+e^{-\beta y^2}cos\left(by\right)$$

$$V(-x,y)=-e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Comparing $V(-x,y)$ with $V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$, I can see that $V(-x,y)\ne V(x,y)$ unless $e^{-\alpha x^2}sin\left(ax\right)=0$ for all x, which is not generally true.

Therefore, the potential is not symmetric under reflection about the y-axis ($x\to -x$), unless specific conditions on the parameters a and $\alpha$ are met such that $e^{-\alpha x^2}sin\left(ax\right)=0$.

Let the potential function be

$$V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Now, I substitute y with $-y$:

$$V(x,-y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta {(-y)}^2}cos\left(b(-y)\right)$$

I know that ${(-y)}^2=y^2$ and $cos(-by)=cos\left(by\right)$. Therefore,

$$V(x,-y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Comparing $V(x,-y)$ with $V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$, I can see that

$$V(x,-y)=V(x,y)$$

Therefore, the potential is symmetric under reflection about the x-axis ($y\to -y$).

Let the potential function be

$$V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Now, I substitute x with $x+a_x$:

$$V(x+a_x,y)=e^{-\alpha {(x+a_x)}^2}sin\left(a(x+a_x)\right)+e^{-\beta y^2}cos\left(by\right)$$

Expanding the terms, I get:

$$V(x+a_x,y)=e^{-\alpha (x^2+2a_{x}x+a_x^2)}sin(ax+a{a}_x)+e^{-\beta y^2}cos\left(by\right)$$

For translational symmetry in the x-direction, I require $V(x+a_x,y)=V(x,y)$ for any arbitrary $a_x$. This would mean:

$$e^{-\alpha (x^2+2a_{x}x+a_x^2)}sin(ax+a{a}_x)+e^{-\beta y^2}cos\left(by\right)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Subtracting the second term from both sides, I need:

$$e^{-\alpha (x^2+2a_{x}x+a_x^2)}sin(ax+a{a}_x)=e^{-\alpha x^2}sin\left(ax\right)$$

This equality does not hold for arbitrary $a_x$ and x unless $\alpha =0$ and $a=0$, in which case the first term vanishes. If $\alpha =0$ and $a\ne 0$, then I would need $sin\left(a(x+a_x)\right)=sin\left(ax\right)$, which implies $a{a}_x=2n\pi$ for some integer n. This must hold for any $a_x$, which is not possible. If $\alpha \ne 0$, the exponential term breaks the translational symmetry as $e^{-\alpha (x^2+2a_{x}x+a_x^2)}\ne e^{-\alpha x^2}$ for $a_x\ne 0$.

Therefore, the potential is not symmetric under translation in the x-direction.

Let the potential function be

$$V(x,y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Now, I substitute y with $y+a_y$:

$$V(x,y+a_y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta {(y+a_y)}^2}cos\left(b(y+a_y)\right)$$

Expanding the terms, I get:

$$V(x,y+a_y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta (y^2+2a_{y}y+a_y^2)}cos(by+b{a}_y)$$

For translational symmetry in the y-direction, I require $V(x,y+a_y)=V(x,y)$ for any arbitrary $a_y$. This would mean:

$$e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta (y^2+2a_{y}y+a_y^2)}cos(by+b{a}_y)=e^{-\alpha x^2}sin\left(ax\right)+e^{-\beta y^2}cos\left(by\right)$$

Subtracting the first term from both sides, I need:

$$e^{-\beta (y^2+2a_{y}y+a_y^2)}cos(by+b{a}_y)=e^{-\beta y^2}cos\left(by\right)$$

This equality does not hold for arbitrary $a_y$ and y unless $\beta =0$ and $b=0$, in which case the second term is constant. If $\beta =0$ and $b\ne 0$, then I would need $cos\left(b(y+a_y)\right)=cos\left(by\right)$, which implies $b{a}_y=2n\pi$ for some integer n. This must hold for any $a_y$, which is not possible. If $\beta \ne 0$, the exponential term breaks the translational symmetry as $e^{-\beta (y^2+2a_{y}y+a_y^2)}\ne e^{-\beta y^2}$ for $a_y\ne 0$.

Therefore, the potential is not symmetric under translation in the y-direction.

3.2.2. Conservation of Energy

The potential energy at time $n+1$ is $V_{n+1}=e^{-\alpha x_{n+1}^2}sin\left(a{x}_{n+1}\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)$. The expected potential energy $\langle V_{n+1}\rangle$ depends on the probability of the hop in the x-direction.

Case 1: ${\eta }_n>0$ The probability of hopping to $i+1$ is $min(1,{\eta }_n)$, and the probability of staying at i is $1-min(1,{\eta }_n)$.

$$\langle V_{n+1}\rangle =min(1,{\eta }_n)V(x_{i+1},y_{n+1})+(1-min(1,{\eta }_n))V(x_i,y_{n+1})$$

$$\langle V_{n+1}\rangle =min(1,{\eta }_n)\left(e^{-\alpha x_{i+1}^2}sin\left(a{x}_{i+1}\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)\right)+(1-min(1,{\eta }_n))\left(e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)\right)$$

Case 2: ${\eta }_n<0$ The probability of hopping to $i-1$ is $min(1,|{\eta }_n\left|\right)$, and the probability of staying at i is $1-min(1,|{\eta }_n\left|\right)$.

$$\langle V_{n+1}\rangle =min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})$$

$$\langle V_{n+1}\rangle =min(1,|{\eta }_n\left|\right)\left(e^{-\alpha x_{i-1}^2}sin\left(a{x}_{i-1}\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)\right)+(1-min(1,|{\eta }_n\left|\right))\left(e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)\right)$$

Case 3: ${\eta }_n=0$ The probability of staying at i is 1.

$$\langle V_{n+1}\rangle =V(x_i,y_{n+1})=e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)$$

In all cases, $y_{n+1}$ is given by the Euler update $y_{n+1}=y_n+{\displaystyle \frac{p_{yn}}{m}}h$.

The Hamiltonian at time n is:

$$H_n={\displaystyle \frac{p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn}^2}{2m}}+e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y_n^2}cos\left(b{y}_n\right)$$

The average Hamiltonian at time $n+1$ is given by:

$$\langle H_{n+1}\rangle =\langle {\displaystyle \frac{p_{xn+1}^2}{2m}}\rangle +{\displaystyle \frac{p_{yn+1}^2}{2m}}+\langle V_{n+1}\rangle$$

where

$$p_{yn+1}=p_{yn}+\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h$$

and

$$p_{xn+1}=p_{xn}+\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h$$

For the case ${\eta }_n=0$, I have $\langle V_{n+1}\rangle =V(x_i,y_{n+1})=e^{-\alpha x_i^2}sin\left(a{x}_i\right)+e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)$. The difference in Hamiltonians is:

$$\langle H_{n+1}\rangle -H_n={\displaystyle \frac{p_{xn+1}^2-p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2-p_{yn}^2}{2m}}+V(x_i,y_{n+1})-V(x_i,y_n)$$

I expand the squared momenta using the Euler update rules:

$$p_{yn+1}=p_{yn}+\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h$$

$$p_{xn+1}=p_{xn}+\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h$$

So,

$$\begin{array}{cc}\hfill p_{yn+1}^2-p_{yn}^2& ={\left(p_{yn}+\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h\right)}^2-p_{yn}^2\hfill \\ & =p_{yn}^2+2p_{yn}\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h+\mathcal{O}\left(h^2\right)-p_{yn}^2\hfill \\ & =2p_{yn}\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h+\mathcal{O}\left(h^2\right)\hfill \end{array}$$

And,

$$\begin{array}{cc}\hfill p_{xn+1}^2-p_{xn}^2& ={\left(p_{xn}+\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h\right)}^2-p_{xn}^2\hfill \\ & =p_{xn}^2+2p_{xn}\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h+\mathcal{O}\left(h^2\right)-p_{xn}^2\hfill \\ & =2p_{xn}\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h+\mathcal{O}\left(h^2\right)\hfill \end{array}$$

Now consider the potential energy difference:

$$V(x_i,y_{n+1})-V(x_i,y_n)=e^{-\beta y_{n+1}^2}cos\left(b{y}_{n+1}\right)-e^{-\beta y_n^2}cos\left(b{y}_n\right)$$

Expanding $y_{n+1}=y_n+{\displaystyle \frac{p_{yn}}{m}}h+\mathcal{O}\left(h^2\right)$ using Taylor series, I would generally find a term proportional to h.

Therefore, $\langle H_{n+1}\rangle -H_n$ will generally not be zero but will have terms of order h, indicating that the average energy is not conserved by the Euler method for this semi-discrete system as Ill.

Case 1: ${\eta }_n>0$

The average Hamiltonian difference is:

$$\langle H_{n+1}\rangle -H_n={\displaystyle \frac{p_{xn+1}^2-p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2-p_{yn}^2}{2m}}+\langle V_{n+1}\rangle -V(x_i,y_n)$$

$$\langle H_{n+1}\rangle -H_n={\displaystyle \frac{p_{xn+1}^2-p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2-p_{yn}^2}{2m}}+min(1,{\eta }_n)V(x_{i+1},y_{n+1})+(1-min(1,{\eta }_n))V(x_i,y_{n+1})-V(x_i,y_n)$$

Substituting the expressions for the kinetic energy differences (up to order h):

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& \approx {\displaystyle \frac{1}{2m}}\left(2p_{xn}\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h\right)\hfill \\ & +{\displaystyle \frac{1}{2m}}\left(2p_{yn}\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h\right)\hfill \\ & +min(1,{\eta }_n)(A{sin}^2\left(a{x}_{i+1}\right)+B{sin}^4\left(a{x}_{i+1}\right)+C{sin}^2\left(b{y}_{n+1}\right))\hfill \\ & +(1-min(1,{\eta }_n))(A{sin}^2\left(a{x}_i\right)+B{sin}^4\left(a{x}_i\right)+C{sin}^2\left(b{y}_{n+1}\right))\hfill \\ & -(A{sin}^2\left(a{x}_i\right)+B{sin}^4\left(a{x}_i\right)+C{sin}^2\left(b{y}_n\right))\hfill \end{array}$$

Case 2: ${\eta }_n<0$

The average Hamiltonian difference is:

$$\langle H_{n+1}\rangle -H_n={\displaystyle \frac{p_{xn+1}^2-p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2-p_{yn}^2}{2m}}+\langle V_{n+1}\rangle -V(x_i,y_n)$$

$$\langle H_{n+1}\rangle -H_n={\displaystyle \frac{p_{xn+1}^2-p_{xn}^2}{2m}}+{\displaystyle \frac{p_{yn+1}^2-p_{yn}^2}{2m}}+min(1,|{\eta }_n\left|\right)V(x_{i-1},y_{n+1})+(1-min(1,|{\eta }_n\left|\right))V(x_i,y_{n+1})-V(x_i,y_n)$$

Substituting the expressions for the kinetic energy differences (up to order h):

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& \approx {\displaystyle \frac{1}{2m}}\left(2p_{xn}\left(2\alpha x_i{e}^{-\alpha x_i^2}sin\left(a{x}_i\right)-a{e}^{-\alpha x_i^2}cos\left(a{x}_i\right)\right)h\right)\hfill \\ & +{\displaystyle \frac{1}{2m}}\left(2p_{yn}\left(2\beta y_n{e}^{-\beta y_n^2}cos\left(b{y}_n\right)+b{e}^{-\beta y_n^2}sin\left(b{y}_n\right)\right)h\right)\hfill \\ & +min(1,|{\eta }_n\left|\right)(A{sin}^2\left(a{x}_{i-1}\right)+B{sin}^4\left(a{x}_{i-1}\right)+C{sin}^2\left(b{y}_{n+1}\right))\hfill \\ & +(1-min(1,|{\eta }_n\left|\right))(A{sin}^2\left(a{x}_i\right)+B{sin}^4\left(a{x}_i\right)+C{sin}^2\left(b{y}_{n+1}\right))\hfill \\ & -(A{sin}^2\left(a{x}_i\right)+B{sin}^4\left(a{x}_i\right)+C{sin}^2\left(b{y}_n\right))\hfill \end{array}$$

In both cases, I see that $\langle H_{n+1}\rangle -H_n$ is generally not zero. The Euler method introduces errors of order h in the continuous part, and the probabilistic hop introduces differences in the potential energy at $x_i,x_{i+1},$ or $x_{i-1}$. Therefore, the average energy is not conserved for this semi-discrete system using Euler’s method.

3.2.3. Properties of the Probabilistic Step

The probability of transitioning from a discrete x-position i at time n to a position j at time $n+1$ is given by:

$$P(i_{n+1}=j|i_n=i)=\left\{\begin{array}{cc}min(1,{\eta }_n)\hfill & \mathrm{if}j=i+1\mathrm{and}{\eta }_n>0\hfill \\ 0\hfill & \mathrm{if}j=i+1\mathrm{and}{\eta }_n\le 0\hfill \\ min(1,|{\eta }_n\left|\right)\hfill & \mathrm{if}j=i-1\mathrm{and}{\eta }_n<0\hfill \\ 0\hfill & \mathrm{if}j=i-1\mathrm{and}{\eta }_n\ge 0\hfill \\ 1-min(1,|{\eta }_n\left|\right)\hfill & \mathrm{if}j=i\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

where ${\eta }_n={\displaystyle \frac{p_{xn+1}h}{m\Delta x}}$.

The expected change in the discrete x-coordinate over one time step, $\langle \Delta i\rangle =\langle i_{n+1}-i_n\rangle$, can be calculated as: If ${\eta }_n>0$:

$$\langle \Delta i\rangle =(i_n+1-i_n)min(1,{\eta }_n)+(i_n-1-i_n)·0+(i_n-i_n)(1-min(1,{\eta }_n))=min(1,{\eta }_n)$$

If ${\eta }_n<0$:

$$\langle \Delta i\rangle =(i_n+1-i_n)·0+(i_n-1-i_n)min(1,|{\eta }_n\left|\right)+(i_n-i_n)(1-min(1,|{\eta }_n\left|\right))=-min(1,|{\eta }_n\left|\right)$$

If ${\eta }_n=0$:

$$\langle \Delta i\rangle =0$$

In summary:

$$\langle \Delta i\rangle =\mathrm{sgn}\left({\eta }_n\right)min(1,|{\eta }_n\left|\right)$$

The variance in the x-position after one time step can be expressed as $\langle {(\Delta i)}^2\rangle -{\langle \Delta i\rangle }^2$. I calculate $\langle {(\Delta i)}^2\rangle =\langle {(i_{n+1}-i_n)}^2\rangle$:

If ${\eta }_n>0$:

$$\langle {(\Delta i)}^2\rangle ={\left(1\right)}^{2}min(1,{\eta }_n)+{(-1)}^2·0+{\left(0\right)}^2(1-min(1,{\eta }_n))=min(1,{\eta }_n)$$

If ${\eta }_n<0$:

$$\langle {(\Delta i)}^2\rangle ={\left(1\right)}^2·0+{(-1)}^{2}min(1,|{\eta }_n\left|\right)+{\left(0\right)}^2(1-min(1,|{\eta }_n\left|\right))=min(1,|{\eta }_n\left|\right)$$

If ${\eta }_n=0$:

$$\langle {(\Delta i)}^2\rangle =0$$

So, $\langle {(\Delta i)}^2\rangle =min(1,|{\eta }_n\left|\right)$.

The variance is then: If ${\eta }_n>0$: $min(1,{\eta }_n)-{(min(1,{\eta }_n))}^2=min(1,{\eta }_n)(1-min(1,{\eta }_n))$ If ${\eta }_n<0$: $min(1,|{\eta }_n\left|\right)-(-min(1,|{\eta }_n{\left|\right))}^2=min(1,|{\eta }_n\left|\right)(1-min(1,|{\eta }_n\left|\right))$ If ${\eta }_n=0$: $0-0^2=0$

The probabilistic step can be vieId as a biased random walk on a 1D lattice. The bias is determined by the parameter ${\eta }_n$. When ${\eta }_n>0$, there is a higher probability of moving to the right, and when ${\eta }_n<0$, there is a higher probability of moving to the left. The magnitude $|{\eta }_n|$ controls the strength of this bias (up to a maximum probability of 1). When ${\eta }_n=0$, the walk has no bias in either direction (although the probability of hopping is zero in this specific implementation).

Conclusion:

Investigation of the second potential function $V(x,y)=A{sin}^2\left(ax\right)+B{sin}^4\left(ax\right)+C{sin}^2\left(by\right)$ using the semi-discretization method shoId that I could successfully apply the methodology. HoIver, similar to the first potential, the average energy of the semi-discrete system is generally not conserved when evolved using the Euler method due to the inherent non-conservation of Euler’s method for continuous Hamiltonian systems and the probabilistic nature of the hopping in the x-direction. The probabilistic step itself functions as a biased random walk on the discrete x-lattice, driven by the x-momentum.

4. Semi-Discrete Dynamics on a 2D Lattice with a Discrete Cosine Potential

Solving Methodology:

The methodology for this case involves a full discretization of the spatial coordinates, while momenta remain continuous. I use finite differences to approximate the forces from the discrete potential and Euler’s method to update the momenta. The position update is probabilistic, with hopping on the 2D lattice. The analytical investigation includes analyzing the discrete symmetries of the lattice and potential, examining the conservation of average energy considering the probabilistic hops, determining conditions for a stationary average position on the lattice, and investigating potential diffusion-like behavior [4].

The Potential Function :

$$V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$$

4.1. Semi-Discretization Methodology

4.1.1. Discrete Positions on the 2D Lattice

I are given the discrete cosine potential function:

$$V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$$

I approximate the forces using finite differences:

$$F_i(i,j)\approx -{\displaystyle \frac{V(i+1,j)-V(i-1,j)}{2\Delta x}}$$

$$F_j(i,j)\approx -{\displaystyle \frac{V(i,j+1)-V(i,j-1)}{2\Delta y}}$$

First, I evaluate the potential at $(i+1,j)$ and $(i-1,j)$:

$$V(i+1,j)=Acos(k_x(i+1)\Delta x+k_{y}j\Delta y+\varphi )=Acos(k_{x}i\Delta x+k_x\Delta x+k_{y}j\Delta y+\varphi )$$

$$V(i-1,j)=Acos(k_x(i-1)\Delta x+k_{y}j\Delta y+\varphi )=Acos(k_{x}i\Delta x-k_x\Delta x+k_{y}j\Delta y+\varphi )$$

Now, substitute these into the expression for $F_i(i,j)$:

$$F_i(i,j)\approx -{\displaystyle \frac{Acos(k_{x}i\Delta x+k_x\Delta x+k_{y}j\Delta y+\varphi )-Acos(k_{x}i\Delta x-k_x\Delta x+k_{y}j\Delta y+\varphi )}{2\Delta x}}$$

I use the trigonometric identity: $cos(a+b)-cos(a-b)=-2sin\left(a\right)sin\left(b\right)$. Let $a=k_{x}i\Delta x+k_{y}j\Delta y+\varphi$ and $b=k_x\Delta x$.

$$F_i(i,j)\approx -{\displaystyle \frac{A(-2sin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_x\Delta x))}{2\Delta x}}$$

$$F_i(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_x\Delta x)}{\Delta x}}$$

Next, I evaluate the potential at $(i,j+1)$ and $(i,j-1)$:

$$V(i,j+1)=Acos(k_{x}i\Delta x+k_y(j+1)\Delta y+\varphi )=Acos(k_{x}i\Delta x+k_{y}j\Delta y+k_y\Delta y+\varphi )$$

$$V(i,j-1)=Acos(k_{x}i\Delta x+k_y(j-1)\Delta y+\varphi )=Acos(k_{x}i\Delta x+k_{y}j\Delta y-k_y\Delta y+\varphi )$$

Now, substitute these into the expression for $F_j(i,j)$:

$$F_j(i,j)\approx -{\displaystyle \frac{Acos(k_{x}i\Delta x+k_{y}j\Delta y+k_y\Delta y+\varphi )-Acos(k_{x}i\Delta x+k_{y}j\Delta y-k_y\Delta y+\varphi )}{2\Delta y}}$$

I use the trigonometric identity: $cos(a+b)-cos(a-b)=-2sin\left(a\right)sin\left(b\right)$. Let $a=k_{x}i\Delta x+k_{y}j\Delta y+\varphi$ and $b=k_y\Delta y$.

$$F_j(i,j)\approx -{\displaystyle \frac{A(-2sin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_y\Delta y))}{2\Delta y}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_y\Delta y)}{\Delta y}}$$

The forces from the discrete cosine potential, using finite difference approximations, are:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_x\Delta x)}{\Delta x}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_y\Delta y)}{\Delta y}}$$

4.1.2. Momentum Update Using Euler’s Method

The Euler update rules for the momenta are given by:

$$p_{i,n+1}=p_{i,n}+F_i(i_n,j_n)h$$

$$p_{j,n+1}=p_{j,n}+F_j(i_n,j_n)h$$

I found the forces to be:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_x\Delta x)}{\Delta x}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_y\Delta y)}{\Delta y}}$$

Substituting these into the Euler update rules, I get:

$$p_{i,n+1}=p_{i,n}+\left({\displaystyle \frac{Asin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_x\Delta x)}{\Delta x}}\right)h$$

$$p_{j,n+1}=p_{j,n}+\left({\displaystyle \frac{Asin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_y\Delta y)}{\Delta y}}\right)h$$

The momenta at the next time step $n+1$ are updated from the momenta at the current time step n according to the following rules:

$$p_{i,n+1}=p_{i,n}+{\displaystyle \frac{Ah}{\Delta x}}sin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_x\Delta x)$$

$$p_{j,n+1}=p_{j,n}+{\displaystyle \frac{Ah}{\Delta y}}sin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_y\Delta y)$$

4.1.3. Probabilistic Position Update (Hopping)

Let the current lattice site at time step n be $(i_n,j_n)$. The parameters influencing the hopping probabilities are:

$${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$$

$${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$$

The probabilities of transitioning to a new site $(i_{n+1},j_{n+1})$ at time step $n+1$ are:

• Probability of hopping in the $+i$ direction:

  $$P(i_{n+1}=i_n+1,j_{n+1}=j_n)=max(0,min(0.25,{\eta }_{i,n}))$$
• Probability of hopping in the $-i$ direction:

  $$P(i_{n+1}=i_n-1,j_{n+1}=j_n)=max(0,min(0.25,-{\eta }_{i,n}))$$
• Probability of hopping in the $+j$ direction:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n+1)=max(0,min(0.25,{\eta }_{j,n}))$$
• Probability of hopping in the $-j$ direction:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n-1)=max(0,min(0.25,-{\eta }_{j,n}))$$
• Probability of staying at the same site:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n)=1-P_{+i}-P_{-i}-P_{+j}-P_{-j}$$

The position on the 2D lattice is updated probabilistically at each time step based on the values of ${\eta }_{i,n}$ and ${\eta }_{j,n}$, which depend on the momenta. The particle has a chance to move one step in the positive or negative i or j direction, or to stay at its current location. The maximum probability of hopping in any single direction is limited to 0.25.

4.2. Analysis

4.2.1. Discrete Symmetry Analysis

For the potential $V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$ to be invariant under the lattice translations $(i,j)\to (i+n_x,j+n_y)$, where $n_x$ and $n_y$ are integers, I must have:

$$V(i+n_x,j+n_y)=V(i,j)$$

Substituting the potential function, this condition becomes:

$$Acos(k_x(i+n_x)\Delta x+k_y(j+n_y)\Delta y+\varphi )=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$$

This equality holds if and only if the arguments of the cosine functions differ by an integer multiple of $2\pi$:

$$(k_x(i+n_x)\Delta x+k_y(j+n_y)\Delta y+\varphi )-(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )=2\pi m$$

where m is an integer. Simplifying the left side, I get:

$$k_x{n}_x\Delta x+k_y{n}_y\Delta y=2\pi m$$

This equation must hold for all integers $n_x$ and $n_y$.

Consider the case where $n_y=0$. Then I need $k_x{n}_x\Delta x=2\pi m$. This must hold for any integer $n_x$. The smallest non-zero integer is $n_x=1$, so I must have $k_x\Delta x=2\pi m_x$ for some integer $m_x$. For this to hold for all $n_x$, it implies that $k_x\Delta x$ must be an integer multiple of $2\pi$.

Similarly, consider the case where $n_x=0$. Then I need $k_y{n}_y\Delta y=2\pi m$. For this to hold for any integer $n_y$, it implies that $k_y\Delta y$ must be an integer multiple of $2\pi$.

Therefore, the conditions for the potential to be invariant under the translational symmetries of the lattice are:

$$k_x\Delta x={\displaystyle \frac{2\pi m_x}{L_x}}{L}_x=2\pi m_x$$

$$k_y\Delta y={\displaystyle \frac{2\pi m_y}{L_y}}{L}_y=2\pi m_y$$

where $m_x$ and $m_y$ are integers, and $L_x=1,L_y=1$ as I are considering a shift by one lattice unit. More generally, if I consider periodicity over $L_x$ and $L_y$ lattice units, then $k_x\Delta x={\displaystyle \frac{2\pi m_x}{L_x}}$ and $k_y\Delta y={\displaystyle \frac{2\pi m_y}{L_y}}$ where $m_x,m_y,L_x,L_y$ are integers. For invariance under a shift by one lattice unit ($n_x=1,n_y=1$), I need $k_x\Delta x=2\pi m_x$ and $k_y\Delta y=2\pi m_y$.

The discrete cosine potential $V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$ exhibits the translational symmetries of the 2D lattice if and only if $k_x\Delta x$ is an integer multiple of $2\pi$ and $k_y\Delta y$ is an integer multiple of $2\pi$.

4.2.2. Conservation of Average Energy

The expected potential energy at time step $n+1$ is given by:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =P(i_{n+1}=i_n+1,j_{n+1}=j_n)V(i_n+1,j_n)\hfill \\ & +P(i_{n+1}=i_n-1,j_{n+1}=j_n)V(i_n-1,j_n)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n+1)V(i_n,j_n+1)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n-1)V(i_n,j_n-1)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n)V(i_n,j_n)\hfill \end{array}$$

Substituting the probabilities I defined:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =max(0,min(0.25,{\eta }_{i,n}))V(i_n+1,j_n)\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))V(i_n-1,j_n)\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))V(i_n,j_n+1)\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))V(i_n,j_n-1)\hfill \\ & +1-max(0,min(0.25,{\eta }_{i,n}))-max(0,min(0.25,-{\eta }_{i,n}))-max(0,min(0.25,{\eta }_{j,n}))\hfill \\ & -max(0,min(0.25,-{\eta }_{j,n}))\hfill \\ & V(i_n,j_n)\hfill \end{array}$$

Now, I substitute the potential function $V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =max(0,min(0.25,{\eta }_{i,n}))Acos(k_x(i_n+1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))Acos(k_x(i_n-1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n+1)\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n-1)\Delta y+\varphi )\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)Acos(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \end{array}$$

The expected potential energy at the next time step is a Iighted average of the potential evaluated at the current site and its fThis nearest neighbors, with the Iights given by the probabilities of hopping or staying.

Difference in Average Energy $\langle H_{n+1}\rangle -H_n$

The difference in average energy betIen time step $n+1$ and n is:

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& =\left({\displaystyle \frac{1}{2m}}{\left(p_{i,n}+{\displaystyle \frac{Ah}{\Delta x}}sin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_x\Delta x)\right)}^2\right.\hfill \\ & +{\displaystyle \frac{1}{2m}}{\left(p_{j,n}+{\displaystyle \frac{Ah}{\Delta y}}sin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )sin(k_y\Delta y)\right)}^2\hfill \\ & +max(0,min(0.25,{\eta }_{i,n}))Acos(k_x(i_n+1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))Acos(k_x(i_n-1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n+1)\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n-1)\Delta y+\varphi )\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)Acos(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & \left.-{\displaystyle \frac{p_{i,n}^2}{2m}}-{\displaystyle \frac{p_{j,n}^2}{2m}}-Acos(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )\right)\hfill \end{array}$$

where ${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$, ${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$, and $P_{+i}=max(0,min(0.25,{\eta }_{i,n}))$, $P_{-i}=max(0,min(0.25,-{\eta }_{i,n}))$, $P_{+j}=max(0,min(0.25,{\eta }_{j,n}))$, $P_{-j}=max(0,min(0.25,-{\eta }_{j,n}))$.

Expanding the squared momentum terms:

$$\begin{array}{cc}\hfill \langle H_{n+1}\rangle -H_n& =\left({\displaystyle \frac{p_{i,n}^2}{2m}}+{\displaystyle \frac{p_{i,n}Ah}{m\Delta x}}sin(...)sin(k_x\Delta x)+{\displaystyle \frac{1}{2m}}{\left({\displaystyle \frac{Ah}{\Delta x}}sin(...)sin(k_x\Delta x)\right)}^2\right.\hfill \\ & +{\displaystyle \frac{p_{j,n}^2}{2m}}+{\displaystyle \frac{p_{j,n}Ah}{m\Delta y}}sin(...)sin(k_y\Delta y)+{\displaystyle \frac{1}{2m}}{\left({\displaystyle \frac{Ah}{\Delta y}}sin(...)sin(k_y\Delta y)\right)}^2\hfill \\ & +max(0,min(0.25,{\eta }_{i,n}))Acos(k_x(i_n+1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))Acos(k_x(i_n-1)\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n+1)\Delta y+\varphi )\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))Acos(k_x{i}_n\Delta x+k_y(j_n-1)\Delta y+\varphi )\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)Acos(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )\hfill \\ & \left.-{\displaystyle \frac{p_{i,n}^2}{2m}}-{\displaystyle \frac{p_{j,n}^2}{2m}}-Acos(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )\right)\hfill \end{array}$$

where $sin(...)$ represents $sin(k_x{i}_n\Delta x+k_y{j}_n\Delta y+\varphi )$.

Summary of Energy Conservation

By observing this expression, I can see that $\langle H_{n+1}\rangle -H_n$ is generally not equal to zero. The terms proportional to h arise from the Euler method, indicating that the continuous part of the dynamics is not energy-conserving. Furthermore, the expected potential energy at time $n+1$ is generally different from the potential energy at time n due to the probabilistic hopping, introducing further terms that do not cancel out. Therefore, the average energy is not conserved for this semi-discrete system using this methodology.

Conclusion Regarding Conservation of Average Energy:

The average energy of the system is generally not conserved under the described semi-discretization methodology and the use of Euler’s method for updating the momenta. This lack of conservation arises from two primary sThisces:

1. Non-Conservation of Euler’s Method: The Euler method, being a first-order approximation, does not generally conserve energy for continuous-time Hamiltonian systems. This introduces terms proportional to the time step h in the energy difference. 2. Probabilistic Hopping and Potential Energy Changes: The probabilistic movement of the particle on the lattice leads to the expected potential energy at the next time step being different from the potential energy at the current time step. This difference, along with the kinetic energy changes, does not generally sum to zero.

Therefore, while this semi-discretization method allows for the evolution of the system, it does not inherently conserve the average energy.

4.2.3. Conditions for a Stationary Average Position

For a stationary average position on the lattice, I require that the expected change in both the i and j coordinates over one time step is zero:

$$\langle \Delta i\rangle =\langle i_{n+1}-i_n\rangle =0$$

$$\langle \Delta j\rangle =\langle j_{n+1}-j_n\rangle =0$$

The expected change in the i-coordinate is given by:

$$\langle \Delta i\rangle =(i_n+1-i_n)P_{+i}+(i_n-1-i_n)P_{-i}+(i_n-i_n)(1-P_{+i}-P_{-i}-P_{+j}-P_{-j})$$

$$\langle \Delta i\rangle =P_{+i}-P_{-i}$$

Substituting the expressions for $P_{+i}$ and $P_{-i}$:

$$\langle \Delta i\rangle =max(0,min(0.25,{\eta }_{i,n}))-max(0,min(0.25,-{\eta }_{i,n}))$$

For $\langle \Delta i\rangle =0$, I must have:

$$max(0,min(0.25,{\eta }_{i,n}))=max(0,min(0.25,-{\eta }_{i,n}))$$

This condition is satisfied if and only if ${\eta }_{i,n}=0$.

Similarly, the expected change in the j-coordinate is given by:

$$\langle \Delta j\rangle =(j_n+1-j_n)P_{+j}+(j_n-1-j_n)P_{-j}+(j_n-j_n)(1-P_{+i}-P_{-i}-P_{+j}-P_{-j})$$

$$\langle \Delta j\rangle =P_{+j}-P_{-j}$$

Substituting the expressions for $P_{+j}$ and $P_{-j}$:

$$\langle \Delta j\rangle =max(0,min(0.25,{\eta }_{j,n}))-max(0,min(0.25,-{\eta }_{j,n}))$$

For $\langle \Delta j\rangle =0$, I must have:

$$max(0,min(0.25,{\eta }_{j,n}))=max(0,min(0.25,-{\eta }_{j,n}))$$

This condition is satisfied if and only if ${\eta }_{j,n}=0$.

Recall that ${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$ and ${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$. Therefore, the conditions for a stationary average position are:

$${\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}=0\Rightarrow p_{i,n+1}=0$$

$${\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}=0\Rightarrow p_{j,n+1}=0$$

The average position on the 2D lattice remains stationary if and only if the momenta at the next time step in both the i and j directions are zero. This implies that there is no net drift of the particle on the lattice in either direction, on average.

4.2.4. Analysis of Diffusion-like Behavior

To explore diffusion-like behavior, I consider the mean squared displacement (MSD). The expected squared change in the i and j coordinates at each step are:

$$\langle {(\Delta i_n)}^2\rangle =P_{+i}+P_{-i}=min(0.25,max(0,{\eta }_{i,n}))+min(0.25,max(0,-{\eta }_{i,n}))$$

$$\langle {(\Delta j_n)}^2\rangle =P_{+j}+P_{-j}=min(0.25,max(0,{\eta }_{j,n}))+min(0.25,max(0,-{\eta }_{j,n}))$$

The MSD after N time steps is $\langle R_N^2\rangle =\langle i_N^2\rangle +\langle j_N^2\rangle ={\sum }_{n=0}^{N-1}(\langle {(\Delta i_n)}^2\rangle +\langle {(\Delta j_n)}^2\rangle )$, assuming the average position starts at zero and the changes in i and j are uncorrelated and have zero mean on average.

The momenta evolve as:

$$p_{i,n+1}=p_{i,n}+F_i(i_n,j_n)h$$

$$p_{j,n+1}=p_{j,n}+F_j(i_n,j_n)h$$

where the forces are:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_x\Delta x)}{\Delta x}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Asin(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )sin(k_y\Delta y)}{\Delta y}}$$

Assuming average initial momentum is zero and the average force over the lattice is zero, the squared momentum in the i-direction evolves as:

$$\langle p_{i,n+1}^2\rangle \approx \langle p_{i,n}^2\rangle +{\displaystyle \frac{h^2}{2}}{\left({\displaystyle \frac{Asin(k_x\Delta x)}{\Delta x}}\right)}^2$$

This implies that the average of the squared momentum grows linearly with time step n. Consequently, $\langle {\eta }_{i,n}^2\rangle$ will also grow.

As $\langle {\eta }_{i,n}^2\rangle$ increases, $\langle {(\Delta i_n)}^2\rangle$ tends towards saturation. If $\langle {(\Delta i_n)}^2\rangle$ and $\langle {(\Delta j_n)}^2\rangle$ reach approximately constant values, then the MSD, $\langle R_N^2\rangle$, grows linearly with N, indicating diffusion-like behavior.

HoIver, the periodic potential might introduce trapping or oscillations, complicating the simple linear growth of MSD. A rigorous analysis requires further investigation.

The semi-discrete dynamics on the 2D lattice with a discrete cosine potential can exhibit diffusion-like behavior if the mean squared displacement of the particle grows linearly with time. This is expected to occur when the average squared momentum increases over time due to the fluctuating forces from the potential, leading to a roughly constant rate of hopping. HoIver, the periodic nature of the potential can also lead to deviations from simple diffusion, such as trapping or oscillations. A detailed characterization of the long-time behavior would require more advanced analytical techniques or numerical simulations.

Summary of Observations and Results for Semi-Discrete Dynamics on a 2D Lattice with a Discrete Cosine Potential

• Semi-Discretization: I successfully applied a semi-discretization method where the spatial coordinates are discrete, momenta are continuous, forces are approximated using finite differences, momenta are updated using Euler’s method, and the position on the 2D lattice evolves probabilistically with defined hopping rules.
• Discrete Symmetries: The potential

  $$V(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$$

  exhibits the translational symmetries of the 2D lattice if $k_x\Delta x$ and $k_y\Delta y$ are integer multiples of $2\pi$.
• Conservation of Average Energy: This analysis indicated that the average energy of the system is generally not conserved due to the inherent non-conservation of the Euler method for the continuous momenta and the changes in potential energy arising from the probabilistic hopping on the lattice.
• Stationary Average Position: The condition for a stationary average position on the lattice (i.e., no net drift) is that the average momenta in both the i and j directions are zero, leading to ${\eta }_{i,n}=0$ and ${\eta }_{j,n}=0$.
• Potential Diffusion-Like Behavior: I explored the possibility of diffusion-like behavior by looking at the mean squared displacement (MSD). The analysis suggests that if the average squared momentum increases over time due to the fluctuating forces from the potential, the MSD could grow linearly with time, indicating diffusion. HoIver, the periodic nature of the cosine potential might also lead to more complex phenomena like trapping or oscillations, which could deviate from simple diffusion. A more detailed investigation would likely require further analytical techniques or numerical simulations.

5. Analysis of Semi-Discrete Motion in a 2D Lattice Potential Exhibiting Permutation Symmetry

Solving Methodology:

The methodology for this case involves a full discretization of the spatial coordinates, while momenta remain continuous. I use finite differences to approximate the forces from the discrete potential and Euler’s method to update the momenta. The position update is probabilistic, with hopping on the 2D lattice. The analytical investigation includes analyzing the discrete symmetries of the lattice and potential, examining the conservation of average energy considering the probabilistic hops, determining conditions for a stationary average position on the lattice, and investigating potential diffusion-like behavior [5].

The Potential Function :

$$V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$$

5.1. Semi-discretization Methodology

5.1.1. Discrete Positions on the 2D Lattice

Forces from the Discrete Potential

I are given the discrete cosine potential function:

$$V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$$

I approximate the forces using finite differences:

$$F_i(i,j)\approx -{\displaystyle \frac{V(i+1,j)-V(i-1,j)}{2\Delta x}}$$

$$F_j(i,j)\approx -{\displaystyle \frac{V(i,j+1)-V(i,j-1)}{2\Delta y}}$$

Calculation of $F_i(i,j)$

First, I evaluate the potential at $(i+1,j)$ and $(i-1,j)$:

$$V(i+1,j)=Acos\left(k\right(i+1)\Delta x)cos(kj\Delta y)=Acos(ki\Delta x+k\Delta x)cos(kj\Delta y)$$

$$V(i-1,j)=Acos\left(k\right(i-1)\Delta x)cos(kj\Delta y)=Acos(ki\Delta x-k\Delta x)cos(kj\Delta y)$$

Now, substitute these into the expression for $F_i(i,j)$:

$$F_i(i,j)\approx -{\displaystyle \frac{Acos(ki\Delta x+k\Delta x)cos(kj\Delta y)-Acos(ki\Delta x-k\Delta x)cos(kj\Delta y)}{2\Delta x}}$$

Factor out the common term $Acos(kj\Delta y)$:

$$F_i(i,j)\approx -{\displaystyle \frac{Acos(kj\Delta y)(cos(ki\Delta x+k\Delta x)-cos(ki\Delta x-k\Delta x\left)\right)}{2\Delta x}}$$

I use the trigonometric identity: $cos(a+b)-cos(a-b)=-2sin\left(a\right)sin\left(b\right)$. Let $a=ki\Delta x$ and $b=k\Delta x$.

$$F_i(i,j)\approx -{\displaystyle \frac{Acos(kj\Delta y)(-2sin(ki\Delta x)sin(k\Delta x\left)\right)}{2\Delta x}}$$

$$F_i(i,j)\approx {\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

Calculation of $F_j(i,j)$

Next, I evaluate the potential at $(i,j+1)$ and $(i,j-1)$:

$$V(i,j+1)=Acos(ki\Delta x)cos\left(k\right(j+1)\Delta y)=Acos(ki\Delta x)cos(kj\Delta y+k\Delta y)$$

$$V(i,j-1)=Acos(ki\Delta x)cos\left(k\right(j-1)\Delta y)=Acos(ki\Delta x)cos(kj\Delta y-k\Delta y)$$

Now, substitute these into the expression for $F_j(i,j)$:

$$F_j(i,j)\approx -{\displaystyle \frac{Acos(ki\Delta x)cos(kj\Delta y+k\Delta y)-Acos(ki\Delta x)cos(kj\Delta y-k\Delta y)}{2\Delta y}}$$

Factor out the common term $Acos(ki\Delta x)$:

$$F_j(i,j)\approx -{\displaystyle \frac{Acos(ki\Delta x)(cos(kj\Delta y+k\Delta y)-cos(kj\Delta y-k\Delta y\left)\right)}{2\Delta y}}$$

I use the trigonometric identity: $cos(a+b)-cos(a-b)=-2sin\left(a\right)sin\left(b\right)$. Let $a=kj\Delta y$ and $b=k\Delta y$.

$$F_j(i,j)\approx -{\displaystyle \frac{Acos(ki\Delta x)(-2sin(kj\Delta y)sin(k\Delta y\left)\right)}{2\Delta y}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

Summary of Results

The forces from the discrete cosine potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$, using finite difference approximations, are:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

5.1.2. Continuous Momenta Update using Euler’s Method

Momentum Update using Euler’s Method

The Euler update rules for the momenta are given by:

$$p_{i,n+1}=p_{i,n}+F_i(i_n,j_n)h$$

$$p_{j,n+1}=p_{j,n}+F_j(i_n,j_n)h$$

I found the forces for the potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ to be:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

Substituting these into the Euler update rules, I get:

$$p_{i,n+1}=p_{i,n}+\left({\displaystyle \frac{Asin(k{i}_n\Delta x)cos(k{j}_n\Delta y)sin(k\Delta x)}{\Delta x}}\right)h$$

$$p_{j,n+1}=p_{j,n}+\left({\displaystyle \frac{Acos(k{i}_n\Delta x)sin(k{j}_n\Delta y)sin(k\Delta y)}{\Delta y}}\right)h$$

Summary of Momentum Update Rules

The momenta at the next time step $n+1$ are updated from the momenta at the current time step n according to the following rules:

$$p_{i,n+1}=p_{i,n}+{\displaystyle \frac{Ahsin(k\Delta x)}{\Delta x}}sin(k{i}_n\Delta x)cos(k{j}_n\Delta y)$$

$$p_{j,n+1}=p_{j,n}+{\displaystyle \frac{Ahsin(k\Delta y)}{\Delta y}}cos(k{i}_n\Delta x)sin(k{j}_n\Delta y)$$

5.1.3. Probabilistic Position Update (Hopping)

Probabilistic Position Update on the 2D Lattice

Let the current lattice site at time step n be $(i_n,j_n)$. The parameters influencing the hopping probabilities are:

$${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$$

$${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$$

The probabilities of transitioning to a new site $(i_{n+1},j_{n+1})$ at time step $n+1$ are:

• Probability of hopping in the $+i$ direction:

  $$P(i_{n+1}=i_n+1,j_{n+1}=j_n)=max(0,min(0.25,{\eta }_{i,n}))$$
• Probability of hopping in the $-i$ direction:

  $$P(i_{n+1}=i_n-1,j_{n+1}=j_n)=max(0,min(0.25,-{\eta }_{i,n}))$$
• Probability of hopping in the $+j$ direction:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n+1)=max(0,min(0.25,{\eta }_{j,n}))$$
• Probability of hopping in the $-j$ direction:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n-1)=max(0,min(0.25,-{\eta }_{j,n}))$$
• Probability of staying at the same site:

  $$P(i_{n+1}=i_n,j_{n+1}=j_n)=1-P_{+i}-P_{-i}-P_{+j}-P_{-j}$$

  where

  $$P_{+i}=max(0,min(0.25,{\eta }_{i,n}))$$

  $$P_{-i}=max(0,min(0.25,-{\eta }_{i,n}))$$

  $$P_{+j}=max(0,min(0.25,{\eta }_{j,n}))$$

  $$P_{-j}=max(0,min(0.25,-{\eta }_{j,n}))$$

Summary of Probabilistic Update Rule

The position on the 2D lattice is updated probabilistically at each time step based on the values of ${\eta }_{i,n}$ and ${\eta }_{j,n}$, which depend on the momenta. The particle has a chance to move one step in the positive or negative i or j direction, or to stay at its current location. The maximum probability of hopping in any single direction is limited to 0.25.

5.2. Analysis

5.2.1. Permutation Symmetry Analysis

Translational Symmetry of the Potential

For the potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ to be invariant under the lattice translations $(i,j)\to (i+n_x,j+n_y)$, where $n_x$ and $n_y$ are integers, I must have:

$$V(i+n_x,j+n_y)=V(i,j)$$

Substituting the potential function, this condition becomes:

$$Acos(k(i+n_x)\Delta x)cos(k(j+n_y)\Delta y)=Acos(ki\Delta x)cos(kj\Delta y)$$

This equality holds if and only if $k{n}_x\Delta x$ is an integer multiple of $2\pi$ and $k{n}_y\Delta y$ is an integer multiple of $2\pi$. That is:

$$k{n}_x\Delta x=2\pi m_x$$

$$k{n}_y\Delta y=2\pi m_y$$

where $m_x$ and $m_y$ are integers. These conditions must hold for all integers $n_x$ and $n_y$. The smallest non-zero integer is $n_x=1$ and $n_y=1$, so for the potential to be invariant under a shift by one lattice unit in both directions, I require:

$$k\Delta x=2\pi m_x$$

$$k\Delta y=2\pi m_y$$

where $m_x$ and $m_y$ are integers.

More generally, if I consider the periodicity of the lattice, the potential will be periodic if $k\Delta x={\displaystyle \frac{2\pi m_x}{L_x}}$ and $k\Delta y={\displaystyle \frac{2\pi m_y}{L_y}}$ where $L_x$ and $L_y$ are the periods in terms of the number of lattice units, and $m_x,m_y$ are integers.

Summary of Translational Symmetry

The discrete cosine potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ exhibits the translational symmetries of the 2D lattice (by one lattice unit in each direction) if and only if $k\Delta x$ is an integer multiple of $2\pi$ and $k\Delta y$ is an integer multiple of $2\pi$.

Permutation Symmetry of the Potential

I will now examine the permutation symmetry of the potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$. Permutation symmetry requires that the potential is invariant under the exchange of indices i and j, i.e., $V(i,j)=V(j,i)$.

The potential function is:

$$V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$$

Swapping the indices i and j, I get:

$$V(j,i)=Acos(kj\Delta x)cos(ki\Delta y)$$

For permutation symmetry, I need $V(i,j)=V(j,i)$, so:

$$Acos(ki\Delta x)cos(kj\Delta y)=Acos(kj\Delta x)cos(ki\Delta y)$$

This equality holds for all integer values of i and j if the lattice spacing is the same in both directions, i.e., $\Delta x=\Delta y$.

I consider a square lattice where $\Delta x=\Delta y=a$. In this case, the potential becomes:

$$V(i,j)=Acos\left(kia\right)cos\left(kja\right)$$

Now, if I exchange i and j:

$$V(j,i)=Acos\left(kja\right)cos\left(kia\right)$$

Since multiplication is commutative, $cos\left(kja\right)cos\left(kia\right)=cos\left(kia\right)cos\left(kja\right)$, which means $V(j,i)=V(i,j)$.

Summary of Permutation Symmetry

The potential function $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ exhibits permutation symmetry ($V(i,j)=V(j,i)$) when the lattice spacing is the same in both the x and y directions ($\Delta x=\Delta y$), which corresponds to a square lattice.

5.2.2. Conservation of Average Energy

Expected Potential Energy at Time $n+1$

The expected potential energy at time step $n+1$ is given by:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =P(i_{n+1}=i_n+1,j_{n+1}=j_n)V(i_n+1,j_n)\hfill \\ & +P(i_{n+1}=i_n-1,j_{n+1}=j_n)V(i_n-1,j_n)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n+1)V(i_n,j_n+1)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n-1)V(i_n,j_n-1)\hfill \\ & +P(i_{n+1}=i_n,j_{n+1}=j_n)V(i_n,j_n)\hfill \end{array}$$

Substituting the probabilities:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =max(0,min(0.25,{\eta }_{i,n}))V(i_n+1,j_n)\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))V(i_n-1,j_n)\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))V(i_n,j_n+1)\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))V(i_n,j_n-1)\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)V(i_n,j_n)\hfill \end{array}$$

Now, I substitute the potential function $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =max(0,min(0.25,{\eta }_{i,n}))Acos(k(i_n+1)\Delta x)cos(k{j}_n\Delta y)\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))Acos(k(i_n-1)\Delta x)cos(k{j}_n\Delta y)\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))Acos(k{i}_n\Delta x)cos(k(j_n+1)\Delta y)\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))Acos(k{i}_n\Delta x)cos(k(j_n-1)\Delta y)\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)Acos(k{i}_n\Delta x)cos(k{j}_n\Delta y)\hfill \end{array}$$

where $P_{+i}=max(0,min(0.25,{\eta }_{i,n}))$, $P_{-i}=max(0,min(0.25,-{\eta }_{i,n}))$, $P_{+j}=max(0,min(0.25,{\eta }_{j,n}))$, $P_{-j}=max(0,min(0.25,-{\eta }_{j,n}))$.

Summary of Expected Potential Energy

The expected potential energy at the next time step is a Iighted average of the potential evaluated at the current site and its fThis nearest neighbors, with the Iights determined by the probabilities of hopping or staying, which depend on the momenta.

Difference in Average Energy $\langle H_{n+1}\rangle -H_n$

The difference in average energy betIen time step $n+1$ and n is:

$$\langle H_{n+1}\rangle -H_n=\left({\displaystyle \frac{p_{i,n+1}^2}{2m}}+{\displaystyle \frac{p_{j,n+1}^2}{2m}}+\langle V(i_{n+1},j_{n+1})\rangle \right)-\left({\displaystyle \frac{p_{i,n}^2}{2m}}+{\displaystyle \frac{p_{j,n}^2}{2m}}+V(i_n,j_n)\right)$$

where

$$p_{i,n+1}=p_{i,n}+{\displaystyle \frac{Ahsin(k\Delta x)}{\Delta x}}sin(k{i}_n\Delta x)cos(k{j}_n\Delta y)$$

$$p_{j,n+1}=p_{j,n}+{\displaystyle \frac{Ahsin(k\Delta y)}{\Delta y}}cos(k{i}_n\Delta x)sin(k{j}_n\Delta y)$$

and $\langle V(i_{n+1},j_{n+1})\rangle$ is:

$$\begin{array}{cc}\hfill \langle V(i_{n+1},j_{n+1})\rangle & =max(0,min(0.25,{\eta }_{i,n}))Acos(k(i_n+1)\Delta x)cos(k{j}_n\Delta y)\hfill \\ & +max(0,min(0.25,-{\eta }_{i,n}))Acos(k(i_n-1)\Delta x)cos(k{j}_n\Delta y)\hfill \\ & +max(0,min(0.25,{\eta }_{j,n}))Acos(k{i}_n\Delta x)cos(k(j_n+1)\Delta y)\hfill \\ & +max(0,min(0.25,-{\eta }_{j,n}))Acos(k{i}_n\Delta x)cos(k(j_n-1)\Delta y)\hfill \\ & +\left(1-P_{+i}-P_{-i}-P_{+j}-P_{-j}\right)Acos(k{i}_n\Delta x)cos(k{j}_n\Delta y)\hfill \end{array}$$

Substituting these into the energy difference equation and expanding the squared momentum terms will result in a lengthy expression. HoIver, similar to the previous case with the single coupled cosine potential, I can expect that $\langle H_{n+1}\rangle -H_n$ will generally not be zero. The Euler method introduces terms of order h, and the probabilistic hopping leads to changes in the expected potential energy that will not perfectly balance the changes in kinetic energy.

Summary of Energy Conservation

As with the previous potential, the average energy is generally not conserved for this semi-discrete system when using Euler’s method for momentum updates and a probabilistic hopping rule. The approximations inherent in Euler’s method and the statistical nature of the position updates lead to a non-zero difference betIen the average energy at consecutive time steps.

5.2.3. Stationary Average Position on the Lattice

Conditions for a Stationary Average Position

For a stationary average position on the lattice, I require that the expected change in both the i and j coordinates over one time step is zero:

$$\langle \Delta i\rangle =\langle i_{n+1}-i_n\rangle =0$$

$$\langle \Delta j\rangle =\langle j_{n+1}-j_n\rangle =0$$

The expected change in the i-coordinate is given by:

$$\langle \Delta i\rangle =(i_n+1-i_n)P_{+i}+(i_n-1-i_n)P_{-i}+(i_n-i_n)(1-P_{+i}-P_{-i}-P_{+j}-P_{-j})$$

$$\langle \Delta i\rangle =P_{+i}-P_{-i}$$

Substituting the expressions for $P_{+i}$ and $P_{-i}$:

$$\langle \Delta i\rangle =max(0,min(0.25,{\eta }_{i,n}))-max(0,min(0.25,-{\eta }_{i,n}))$$

For $\langle \Delta i\rangle =0$, I must have:

$$max(0,min(0.25,{\eta }_{i,n}))=max(0,min(0.25,-{\eta }_{i,n}))$$

This condition is satisfied if and only if ${\eta }_{i,n}=0$.

Similarly, the expected change in the j-coordinate is given by:

$$\langle \Delta j\rangle =(j_n+1-j_n)P_{+j}+(j_n-1-j_n)P_{-j}+(j_n-j_n)(1-P_{+i}-P_{-i}-P_{+j}-P_{-j})$$

$$\langle \Delta j\rangle =P_{+j}-P_{-j}$$

Substituting the expressions for $P_{+j}$ and $P_{-j}$:

$$\langle \Delta j\rangle =max(0,min(0.25,{\eta }_{j,n}))-max(0,min(0.25,-{\eta }_{j,n}))$$

For $\langle \Delta j\rangle =0$, I must have:

$$max(0,min(0.25,{\eta }_{j,n}))=max(0,min(0.25,-{\eta }_{j,n}))$$

This condition is satisfied if and only if ${\eta }_{j,n}=0$.

Recall that ${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$ and ${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$. Therefore, the conditions for a stationary average position are:

$${\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}=0\Rightarrow p_{i,n+1}=0$$

$${\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}=0\Rightarrow p_{j,n+1}=0$$

Summary of Conditions for Stationary Average Position

The average position on the 2D lattice remains stationary if and only if the momenta at the next time step in both the i and j directions are zero. This implies that there is no net drift of the particle on the lattice in either direction, on average.

5.3. Potential Diffusion-like BehaviThis

Analysis of Potential Diffusion-Like Behavior

To explore diffusion-like behavior, I again consider the mean squared displacement (MSD). The expected squared change in the i and j coordinates at each step are:

$$\langle {(\Delta i_n)}^2\rangle =P_{+i}+P_{-i}=min(0.25,max(0,{\eta }_{i,n}))+min(0.25,max(0,-{\eta }_{i,n}))$$

$$\langle {(\Delta j_n)}^2\rangle =P_{+j}+P_{-j}=min(0.25,max(0,{\eta }_{j,n}))+min(0.25,max(0,-{\eta }_{j,n}))$$

The momenta evolve as:

$$p_{i,n+1}=p_{i,n}+{\displaystyle \frac{Ahsin(k\Delta x)}{\Delta x}}sin(k{i}_n\Delta x)cos(k{j}_n\Delta y)$$

$$p_{j,n+1}=p_{j,n}+{\displaystyle \frac{Ahsin(k\Delta y)}{\Delta y}}cos(k{i}_n\Delta x)sin(k{j}_n\Delta y)$$

I consider the average force over the lattice. The average of $sin(ki\Delta x)cos(kj\Delta y)$ and $cos(ki\Delta x)sin(kj\Delta y)$ over a sufficiently large lattice is zero (assuming $k\Delta x$ and $k\Delta y$ are not integer multiples of $2\pi$). Thus, if the initial average momentum is zero, it will remain around zero.

Now, I look at the evolution of the average squared momentum in the i-direction:

$$\langle p_{i,n+1}^2\rangle =\langle {(p_{i,n}+F_i(i_n,j_n)h)}^2\rangle =\langle p_{i,n}^2\rangle +2h\langle p_{i,n}{F}_i(i_n,j_n)\rangle +h^2\langle F_i{(i_n,j_n)}^2\rangle$$

Assuming $\langle p_{i,n}\rangle =0$, the term $\langle p_{i,n}{F}_i(i_n,j_n)\rangle$ might also be zero on average. The last term is:

$$\langle F_i{(i_n,j_n)}^2\rangle ={\left({\displaystyle \frac{Asin(k\Delta x)}{\Delta x}}\right)}^2\langle {sin}^2(k{i}_n\Delta x){cos}^2(k{j}_n\Delta y)\rangle$$

Using the averages $\langle {sin}^2\left(\theta \right)\rangle =1/2$ and $\langle {cos}^2\left(\theta \right)\rangle =1/2$ over a cycle:

$$\langle F_i{(i_n,j_n)}^2\rangle \approx {\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{Asin(k\Delta x)}{\Delta x}}\right)}^2$$

Similarly,

$$\langle F_j{(i_n,j_n)}^2\rangle \approx {\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{Asin(k\Delta y)}{\Delta y}}\right)}^2$$

Thus, $\langle p_{i,n+1}^2\rangle \approx \langle p_{i,n}^2\rangle +{\displaystyle \frac{h^2}{4}}{\left({\displaystyle \frac{Asin(k\Delta x)}{\Delta x}}\right)}^2$, and similarly for $p_j$. This shows that the average squared momenta grow over time. Consequently, $\langle {\eta }_{i,n}^2\rangle$ and $\langle {\eta }_{j,n}^2\rangle$ will also grow.

As $\langle {\eta }^2\rangle$ increases, the values of $\langle {(\Delta i_n)}^2\rangle$ and $\langle {(\Delta j_n)}^2\rangle$ will tend towards saturation. If these values become roughly constant, the MSD will grow linearly with time, indicating diffusion-like behavior.

Summary of Potential Diffusion-Like Behavior

Similar to the previous potential, the semi-discrete dynamics with $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ can also exhibit diffusion-like behavior. The fluctuating forces from the potential cause the average squared momenta to grow, leading to a non-zero probability of hopping. Over time, this can result in a linear growth of the mean squared displacement. The periodic nature of the potential can still influence the details of this diffusion, potentially leading to anisotropic diffusion if $\Delta x\ne \Delta y$ or if the parameters in the i and j directions are different.

5.3.1. Motion Along the Line i = j

I consider the motion along the line $i=j$ for the potential function $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$. I assume a square lattice where $\Delta x=\Delta y=a$, so the potential is:

$$V(i,j)=Acos\left(kia\right)cos\left(kja\right)$$

Potential Along $i=j$

Substituting $j=i$, the potential along this line is:

$$V(i,i)=Acos\left(kia\right)cos\left(kia\right)=A{cos}^2\left(kia\right)$$

This can also be expressed using the double angle formula for cosine:

$$V(i,i)={\displaystyle \frac{A}{2}}(1+cos\left(2kia\right))={\displaystyle \frac{A}{2}}+{\displaystyle \frac{A}{2}}cos\left(2kia\right)$$

This represents a 1D periodic potential along the line $i=j$.

Forces Along $i=j$

For a square lattice, the force components are:

$$F_i(i,j)\approx {\displaystyle \frac{Asin\left(kia\right)cos\left(kja\right)sin\left(ka\right)}{a}}$$

$$F_j(i,j)\approx {\displaystyle \frac{Acos\left(kia\right)sin\left(kja\right)sin\left(ka\right)}{a}}$$

Along the line $i=j$, these become:

$$F_i(i,i)={\displaystyle \frac{Asin\left(kia\right)cos\left(kia\right)sin\left(ka\right)}{a}}={\displaystyle \frac{A}{2a}}sin\left(2kia\right)sin\left(ka\right)$$

$$F_j(i,i)={\displaystyle \frac{Acos\left(kia\right)sin\left(kia\right)sin\left(ka\right)}{a}}={\displaystyle \frac{A}{2a}}sin\left(2kia\right)sin\left(ka\right)$$

Thus, $F_i(i,i)=F_j(i,i)$ along the line $i=j$.

The potential along the line $i=j$ for the given 2D lattice potential on a square lattice is a 1D periodic potential. The forces in the i and j directions are equal along this line and depend on $sin\left(2kia\right)$. The motion of a particle starting on this line, under the defined semi-discretization with independent hopping, may or may not stay on this line.

5.3.2. Analysis of Average Drift

The average drift in the i and j directions is given by the expected change in position per time step:

$$\langle \Delta i\rangle =P_{+i}-P_{-i}=max(0,min(0.25,{\eta }_{i,n}))-max(0,min(0.25,-{\eta }_{i,n}))$$

$$\langle \Delta j\rangle =P_{+j}-P_{-j}=max(0,min(0.25,{\eta }_{j,n}))-max(0,min(0.25,-{\eta }_{j,n}))$$

where ${\eta }_{i,n}={\displaystyle \frac{p_{i,n+1}h}{m\Delta x}}$ and ${\eta }_{j,n}={\displaystyle \frac{p_{j,n+1}h}{m\Delta y}}$, and the momenta are updated by:

$$p_{i,n+1}=p_{i,n}+{\displaystyle \frac{Ahsin(k\Delta x)}{\Delta x}}sin(k{i}_n\Delta x)cos(k{j}_n\Delta y)$$

$$p_{j,n+1}=p_{j,n}+{\displaystyle \frac{Ahsin(k\Delta y)}{\Delta y}}cos(k{i}_n\Delta x)sin(k{j}_n\Delta y)$$

If the initial average momenta are zero, $\langle p_{i,0}\rangle =0$ and $\langle p_{j,0}\rangle =0$, the average force in the i-direction over the lattice is:

$$\langle F_i(i_n,j_n)\rangle ={\displaystyle \frac{Asin(k\Delta x)}{\Delta x}}\langle sin(k{i}_n\Delta x)cos(k{j}_n\Delta y)\rangle$$

If $k\Delta x$ and $k\Delta y$ are not integer multiples of $2\pi$, then $\langle F_i\rangle =0$ and $\langle F_j\rangle =0$. This implies $\langle p_{i,n+1}\rangle =0$ and $\langle p_{j,n+1}\rangle =0$, and consequently $\langle {\eta }_{i,n}\rangle =0$ and $\langle {\eta }_{j,n}\rangle =0$ on average.

Substituting these into the average drift equations gives:

$$\langle \langle \Delta i\rangle \rangle =0$$

$$\langle \langle \Delta j\rangle \rangle =0$$

Therefore, under the assumption of zero initial average momentum and that $k\Delta x$ and $k\Delta y$ are not integer multiples of $2\pi$, the average drift for this potential is zero. A non-zero average drift could occur with non-zero initial average momentum or if the average force is non-zero.

Average Drift with Non-Zero Initial Average Momentum

Suppose at the initial time $n=0$, the particle has a non-zero average momentum $\langle p_{i,0}\rangle =P_{i0}$ and $\langle p_{j,0}\rangle =P_{j0}$. Assuming the average force from the potential is zero, the average momenta at subsequent time steps remain constant:

$$\langle p_{i,n}\rangle =P_{i0}$$

$$\langle p_{j,n}\rangle =P_{j0}$$

This leads to average values of $\eta$:

$$\langle {\eta }_{i,n}\rangle ={\displaystyle \frac{P_{i0}h}{m\Delta x}}$$

$$\langle {\eta }_{j,n}\rangle ={\displaystyle \frac{P_{j0}h}{m\Delta y}}$$

The average drift in the i and j directions becomes:

$$\langle \langle \Delta i\rangle \rangle =max(0,min(0.25,{\displaystyle \frac{P_{i0}h}{m\Delta x}}))-max(0,min(0.25,-{\displaystyle \frac{P_{i0}h}{m\Delta x}}))$$

$$\langle \langle \Delta j\rangle \rangle =max(0,min(0.25,{\displaystyle \frac{P_{j0}h}{m\Delta y}}))-max(0,min(0.25,-{\displaystyle \frac{P_{j0}h}{m\Delta y}}))$$

If $P_{i0}>0$, there is a positive average drift in the i-direction. If $P_{i0}<0$, there is a negative average drift. The magnitude of the drift increases with $|P_{i0}|$ until it saturates at $0.25$ steps per time step. Similar behavior occurs for the drift in the j-direction based on the sign and magnitude of $P_{j0}$.

A non-zero initial average momentum leads to a constant average drift in the direction of the initial momentum. The speed of this drift depends on the magnitude of the initial momentum and system parameters, with a maximum drift rate of $0.25$ lattice sites per time step in each direction.

Average Drift Due to Specific Values of Potential Parameters

The forces are given by:

$$F_i(i,j)\approx {\displaystyle \frac{Asin(k\Delta x)}{\Delta x}}sin(ki\Delta x)cos(kj\Delta y)$$

$$F_j(i,j)\approx {\displaystyle \frac{Asin(k\Delta y)}{\Delta y}}cos(ki\Delta x)sin(kj\Delta y)$$

For a non-zero average force over the lattice, the averages $\langle sin(ki\Delta x)cos(kj\Delta y)\rangle$ or $\langle cos(ki\Delta x)sin(kj\Delta y)\rangle$ must be non-zero. If $k\Delta x=2\pi m_x$ (integer $m_x$), then $sin(ki\Delta x)=0$, so $F_i=0$. Similarly, if $k\Delta y=2\pi m_y$, then $sin(kj\Delta y)=0$, so $F_j=0$.

For a non-zero average of $sin(ki\Delta x)$ or $cos(ki\Delta x)$ over all integers i, $k\Delta x$ would typically need to be a multiple of $2\pi$ (making sine zero) or zero (making sine zero).

HoIver, in a finite lattice of size $L_x\times L_y$, averages like ${\displaystyle \frac{1}{L_x}}{\sum }_{i=0}^{L_x-1}sin(ki\Delta x)$ might be non-zero for specific relationships betIen $k\Delta x$ and $L_x$.

For the infinite lattice case, it is generally challenging to obtain a sustained non-zero average force with this potential through specific values of $k,\Delta x,$ or $\Delta y$ alone.

Achieving a sustained non-zero average force, and thus a constant drift, through specific choices of the potential parameters $k,\Delta x,$ or $\Delta y$ is difficult for this potential over an infinite lattice. Such a drift might be possible in finite lattices or with specific relationships betIen these parameters and the lattice dimensions.

Effects of Boundaries on Average Drift

In a finite lattice with boundaries, the probabilistic hopping rule needs to be modified when the particle reaches an edge. I consider two common types of boundaries: reflecting and absorbing.

Reflecting Boundaries

Consider a lattice bounded by $0\le i<L_x$ and $0\le j<L_y$.

• If the particle is at $i_n=0$ and ${\eta }_{i,n}<0$ (attempting a hop to $i=-1$), a reflecting boundary would prevent this. A simple rule might be that the probability $P_{-i}$ is redistributed to the probability of staying at $i=0$ or hopping in the $+i$ direction. This would create a net average drift away from the boundary at $i=0$.
• Similarly, at $i_n=L_x-1$ with ${\eta }_{i,n}>0$ (attempting a hop to $i=L_x$), the hop would be prevented, and the probability $P_{+i}$ would be redistributed, leading to an average drift away from the boundary at $i=L_x-1$.
• The same logic applies to the boundaries in the j direction at $j=0$ and $j=L_y-1$ based on the sign of ${\eta }_{j,n}$.

In general, reflecting boundaries in a finite lattice with a symmetric potential tend to cause an average drift away from the boundaries.

Absorbing Boundaries

Again, consider a lattice bounded by $0\le i<L_x$ and $0\le j<L_y$.

• If the particle reaches any of the boundaries (e.g., $i=0,i=L_x-1,j=0,j=L_y-1$), it is absorbed and removed from the system.
• If I consider an ensemble of particles, the total number of particles within the lattice will decrease over time.
• The average position of the particles that remain in the lattice might shift depending on the potential and the initial distribution. For a symmetric potential and a symmetric initial distribution, the average position of the surviving particles might stay roughly at the center, but the probability distribution will narrow as particles at the edges are absorbed.

Absorbing boundaries do not directly create a drift of the average position if the potential and initial conditions are symmetric. HoIver, they do change the overall distribution and lead to a loss of particles from the system.

Summary

Reflecting boundaries typically induce an average drift away from the boundaries due to the prevention and redistribution of hopping probabilities that would lead outside the lattice. Absorbing boundaries lead to the removal of particles upon reaching them, which changes the overall distribution but may not cause a net drift in the average position of the remaining particles in cases of symmetry.

Symmetries of the Forces

The force components are given by:

$$F_i(i,j)={\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

$$F_j(i,j)={\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

Permutation Symmetry

Under the exchange of indices $i\leftrightarrow j$ and lattice spacings $\Delta x\leftrightarrow \Delta y$, the forces transform as:

$$F_i(j,i,\Delta y,\Delta x)={\displaystyle \frac{Asin(kj\Delta y)cos(ki\Delta x)sin(k\Delta y)}{\Delta y}}=F_j(i,j,\Delta y,\Delta x)$$

$$F_j(j,i,\Delta y,\Delta x)={\displaystyle \frac{Acos(kj\Delta y)sin(ki\Delta x)sin(k\Delta x)}{\Delta x}}=F_i(i,j,\Delta x,\Delta y)$$

If $\Delta x=\Delta y=a$, then $F_i(j,i)=F_j(i,j)$ and $F_j(j,i)=F_i(i,j)$. The forces are interchanged under permutation of indices on a square lattice.

Reflection Symmetry

Reflection about the i-axis ($j\leftrightarrow -j$):

$$F_i(i,-j)={\displaystyle \frac{Asin(ki\Delta x)cos(-kj\Delta y)sin(k\Delta x)}{\Delta x}}=F_i(i,j)$$

$$F_j(i,-j)={\displaystyle \frac{Acos(ki\Delta x)sin(-kj\Delta y)sin(k\Delta y)}{\Delta y}}=-F_j(i,j)$$

Reflection about the j-axis ($i\leftrightarrow -i$):

$$F_i(-i,j)={\displaystyle \frac{Asin(-ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}=-F_i(i,j)$$

$$F_j(-i,j)={\displaystyle \frac{Acos(-ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}=F_j(i,j)$$

Rotational Symmetry (for $\Delta x=\Delta y=a$)

Rotation by 90 degrees ($(i,j)\to (-j,i)$):

$$F_i(-j,i)={\displaystyle \frac{Asin(-kja)cos\left(kia\right)sin\left(ka\right)}{a}}=-{\displaystyle \frac{Asin\left(kja\right)cos\left(kia\right)sin\left(ka\right)}{a}}=-F_j(i,j)$$

$$F_j(-j,i)={\displaystyle \frac{Acos(-kja)sin\left(kia\right)sin\left(ka\right)}{a}}={\displaystyle \frac{Acos\left(kja\right)sin\left(kia\right)sin\left(ka\right)}{a}}=F_i(j,-i)$$

5.4. Symmetries of the Forces

The force components are given by:

$$F_i(i,j)={\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

$$F_j(i,j)={\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

Permutation Symmetry

Under the exchange of indices $i\leftrightarrow j$ and lattice spacings $\Delta x\leftrightarrow \Delta y$, the forces transform as:

$$F_i(j,i,\Delta y,\Delta x)={\displaystyle \frac{Asin(kj\Delta y)cos(ki\Delta x)sin(k\Delta y)}{\Delta y}}=F_j(i,j,\Delta y,\Delta x)$$

$$F_j(j,i,\Delta y,\Delta x)={\displaystyle \frac{Acos(kj\Delta y)sin(ki\Delta x)sin(k\Delta x)}{\Delta x}}=F_i(i,j,\Delta x,\Delta y)$$

If $\Delta x=\Delta y=a$, then $F_i(j,i)=F_j(i,j)$ and $F_j(j,i)=F_i(i,j)$. The forces are interchanged under permutation of indices on a square lattice.

Reflection Symmetry

Reflection about the i-axis ($j\leftrightarrow -j$):

$$F_i(i,-j)={\displaystyle \frac{Asin(ki\Delta x)cos(-kj\Delta y)sin(k\Delta x)}{\Delta x}}=F_i(i,j)$$

$$F_j(i,-j)={\displaystyle \frac{Acos(ki\Delta x)sin(-kj\Delta y)sin(k\Delta y)}{\Delta y}}=-F_j(i,j)$$

Reflection about the j-axis ($i\leftrightarrow -i$):

$$F_i(-i,j)={\displaystyle \frac{Asin(-ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}=-F_i(i,j)$$

$$F_j(-i,j)={\displaystyle \frac{Acos(-ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}=F_j(i,j)$$

Rotational Symmetry (for $\Delta x=\Delta y=a$)

Rotation by 90 degrees ($(i,j)\to (-j,i)$):

$$F_i(-j,i)={\displaystyle \frac{Asin(-kja)cos\left(kia\right)sin\left(ka\right)}{a}}=-{\displaystyle \frac{Asin\left(kja\right)cos\left(kia\right)sin\left(ka\right)}{a}}=-F_j(i,j)$$

$$F_j(-j,i)={\displaystyle \frac{Acos(-kja)sin\left(kia\right)sin\left(ka\right)}{a}}={\displaystyle \frac{Acos\left(kja\right)sin\left(kia\right)sin\left(ka\right)}{a}}=F_i(j,-i)$$

Summary of Force Symmetries

The forces exhibit specific symmetry properties under permutation and reflection of the lattice indices. On a square lattice, the force components are interchanged upon swapping i and j. Reflections about the axes invert the sign of the force component along the other axis. Rotational symmetry is more complex and involves a relationship betIen the components at the rotated points.

5.5. Symmetries of the Motion

The motion of the particle is governed by the momentum update rules and the probabilistic hopping on the lattice. I consider how these behave under the discrete symmetries.

Reflection Symmetry

Reflection about the i-axis ($j\leftrightarrow -j$): The potential $V(i,-j)=V(i,j)$. The forces transform as $F_i(i,-j)=F_i(i,j)$ and $F_j(i,-j)=-F_j(i,j)$. The momentum update rules become:

$$p_{i,n+1}^{\prime }=p_{i,n}^{\prime }+F_i(i_n^{\prime },-j_n^{\prime })h=p_{i,n}+F_i(i_n,-j_n)h=p_{i,n}+F_i(i_n,j_n)h=p_{i,n+1}$$

$$p_{j,n+1}^{\prime }=p_{j,n}^{\prime }+F_j(i_n^{\prime },-j_n^{\prime })h=-p_{j,n}+(-F_j(i_n,j_n))h=-(p_{j,n}+F_j(i_n,j_n)h)=-p_{j,n+1}$$

The $\eta$ parameters transform as ${\eta }_{i,n}^{\prime }={\eta }_{i,n}$ and ${\eta }_{j,n}^{\prime }=-{\eta }_{j,n}$. This leads to a swap in the hopping probabilities in the j direction: $P_{+j}^{\prime }=P_{-j}$ and $P_{-j}^{\prime }=P_{+j}$, while the probabilities in the i direction remain the same. If the initial distribution of j and $p_j$ is symmetric, the average motion will respect this reflection symmetry.

Reflection about the j-axis ($i\leftrightarrow -i$): The potential $V(-i,j)=V(i,j)$. The forces transform as $F_i(-i,j)=-F_i(i,j)$ and $F_j(-i,j)=F_j(i,j)$. Similar to the previous case, if $p_{i,n}$ changes sign, then $p_{i,n+1}$ changes sign, and the hopping probabilities in the i direction are swapped. The motion respects reflection about the j-axis on average with symmetric initial conditions.

Rotational Symmetry (for $\Delta x=\Delta y=a$)

Under a 90-degree rotation $(i,j)\to (-j,i)$, the potential is invariant. The forces transform, and if the momenta also transform as $p_i^{\prime }=-p_j$ and $p_j^{\prime }=p_i$, a more complex analysis shows that the dynamics should respect the rotational symmetry on average if the initial state is also rotationally symmetric.

Summary of Motion Symmetries

The semi-discrete motion exhibits the discrete symmetries of the underlying potential on average, provided that the initial conditions of the system are also symmetric with respect to the transformations considered. Reflections about the axes lead to corresponding reflections in the average motion, and rotational symmetry of the potential and initial state implies rotational symmetry in the average dynamics.

Conclusion:

Summary of Observations and Results for Semi-Discrete Motion with Permutation Symmetry Potential

• Semi-Discretization: I applied the same semi-discretization methodology as before.
• Forces: I calculated the forces $F_i(i,j)$ and $F_j(i,j)$ using finite difference approximations.
• Momentum Update: I used Euler’s method to update the momenta based on these forces.
• Probabilistic Hopping: The position on the lattice evolves probabilistically based on the updated momenta.
• Discrete Symmetries: The potential $V(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$ exhibits translational symmetry under certain conditions on $k\Delta x$ and $k\Delta y$, permutation symmetry when $\Delta x=\Delta y$, and also possesses reflection and rotational symmetries (for a square lattice).
• Conservation of Average Energy: The average energy is generally not conserved for similar reasons as with the previous potential.
• Stationary Average Position: The average position is stationary if the average momenta are zero.
• Potential Diffusion-Like Behavior: The system has the potential for diffusion-like behavior, with the mean squared displacement growing over time due to the fluctuating momenta.
• Motion Along $i=j$: For a square lattice with lattice spacing a, the potential along this line is $A{cos}^2\left(kia\right)$, and the forces $F_i$ and $F_j$ are equal to ${\displaystyle \frac{A}{2a}}sin\left(2kia\right)sin\left(ka\right)$.
• Average Drift: The average drift is zero if the initial average momentum and average force are zero. A non-zero initial average momentum leads to a constant average drift. The average force from this potential is typically zero over an infinite lattice.
• Variance of the Position: The variance of the position grows over time, indicating diffusion, even when the average drift is zero.
• Symmetries of the Forces: The forces exhibit specific symmetries under permutation and reflection of indices.
• Symmetries of the Motion: The motion respects the symmetries of the potential and forces on average, provided the initial conditions are also symmetric.

6. Results

Thorough Comparison and Contrast of the Two Analyzed Potentials

I have analyzed two potential functions using the semi-discretization methodology:

• $V_1(i,j)=Acos(k_{x}i\Delta x+k_{y}j\Delta y+\varphi )$
• $V_2(i,j)=Acos(ki\Delta x)cos(kj\Delta y)$

I. Comparison and Contrast of Properties within Each Potential:

For both potential functions:

• Semi-Discretization Methodology: The same methodology was applied to both.
• Conservation of Average Energy: Average energy is generally not conserved for both due to the Euler method and probabilistic hopping.
• Stationary Average Position: Requires zero average momenta in both i and j directions for both.
• Potential Diffusion-Like Behavior: Both potentials show the potential for diffusion-like behavior driven by fluctuating momenta.
• Average Drift: For both, average drift is zero under zero initial average momentum and average force. Non-zero initial average momentum leads to drift.

II. Pairwise Comparison and Contrast of the Two Potentials Across Analyzed Properties:

• Discrete Symmetries:

  -

  $V_1$: Translational, rotational (if square lattice and $k_x=k_y$), and reflection symmetries.

  -

  $V_2$: Translational, permutation (if $\Delta x=\Delta y$), rotational (if square lattice), and reflection symmetries.

  -

  Contrast:$V_2$ possesses permutation symmetry on a square lattice, which is not a general symmetry of $V_1$.

• Forces:

  -

  $V_1$: Coupled forces where each component depends on both i and j through the argument $k_{x}i\Delta x+k_{y}j\Delta y+\varphi$.

  -

  $V_2$: Separable forces:

  $$F_i(i,j)\approx {\displaystyle \frac{Asin(ki\Delta x)cos(kj\Delta y)sin(k\Delta x)}{\Delta x}}$$

  $$F_j(i,j)\approx {\displaystyle \frac{Acos(ki\Delta x)sin(kj\Delta y)sin(k\Delta y)}{\Delta y}}$$

  -

  Contrast:$V_1$ produces coupled forces, while $V_2$ produces separable forces.

• Motion Along $i=j$ (Square Lattice, $\Delta x=\Delta y=a$):

  -

  $V_1$: $V_1(i,i)=Acos((k_x+k_y)ia+\varphi )$.

  -

  $V_2$: $V_2(i,i)=A{cos}^2\left(kia\right)={\displaystyle \frac{A}{2}}(1+cos\left(2kia\right))$.

  -

  Contrast: Different forms of 1D periodic potentials along the line $i=j$.

• Symmetries of the Forces:

  -

  The specific symmetry properties under permutation and reflection differed betIen the two potentials, reflecting the forms of their respective force expressions.

• Symmetries of the Motion:

  -

  Both potentials resulted in motion that respects the symmetries of the potential and forces on average, given symmetric initial conditions.

III. Group-wise Comparison and Contrast Across All Analyzed Properties:

The semi-discretization methodology established a common framework, leading to similarities in high-level behaviors such as energy non-conservation, conditions for stationary average position, and the potential for diffusion. The key differences observed in This analysis stem from the distinct mathematical structures of the two potential functions, particularly regarding their symmetries and the nature of the forces they produce (coupled vs. separable).