---
title: "Time Complexity, Inferential uncertainty and Spacekime Analytics"
subtitle: "Action Principles, Hamiltonian-Jacobi Equation, and Kime Representation"
author: "SOCR Team"
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output:
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 3
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
  word_document:
    toc: true
    toc_depth: '3'
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amssymb}
  - \usepackage{physics}
  - \usepackage{mathrsfs}
  - \usepackage{bm}
---

```{r setup, include=FALSE, warning=FALSE, message=FALSE, error=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(dplyr)
```


Let's explore the connections between mathematical and physical interconnections 
between several fundamental concepts in theoretical physics: the action principle,
functional derivatives, group generators, distribution actions, analytic continuation,
and energy conservation laws. 

# The Action Principle in Classical Mechanics

The *action* $S$ is defined as the time integral of the Lagrangian $L$ over the
path of a system from an initial time $t_1$ to a final time $t_2$

$$S[q] = \int_{t_1}^{t_2} L(q(t), \dot{q}(t), t) \, dt ,$$

where $q(t)$ represents generalized coordinates and $\dot{q}(t)$ represents their time derivatives.

For a classical mechanical system with $n$ degrees of freedom, the Lagrangian is

$$L(q_1, \ldots, q_n, \dot{q}_1, \ldots, \dot{q}_n, t) = T(\dot{q}_1, \ldots, \dot{q}_n) - V(q_1, \ldots, q_n, t)$$

where $T$ is the kinetic energy and $V$ is the potential energy of the system.

## Hamilton's Principle of Least Action

Hamilton's principle states that the actual path taken by a physical system between
two fixed configurations is the one that makes the action *stationary* (typically a minimum), i.e., $\delta S[q] = 0$,
where $\delta$ denotes a *variation* of the path while keeping the endpoints fixed.

## Euler-Lagrange Equations

The condition $\delta S = 0$ leads to the Euler-Lagrange equations, i.e.,
the equations of motion for the system

$$\frac{\partial L}{\partial q_i} - \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) = 0, \quad i = 1, \ldots, n .$$

For a particle of mass $m$ moving in a potential $V(\vec{r})$, with $\vec{r} = (x, y, z)$, the Lagrangian is

$$L = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2 + \dot{z}^2) - V(x, y, z) .$$

Applying the Euler-Lagrange equations yields Newton's second law $m\ddot{\vec{r}} = -\vec{\nabla}V(\vec{r})$.

# Functional Derivatives of the Action

The *functional derivative* $\frac{\delta S}{\delta q(t)}$ generalizes the 
concept of partial derivatives to functionals and measures how the action 
changes when the path $q(t)$ is perturbed by an infinitesimal amount $\delta q(t)$.
GIven a *functional* $F[q] = \int_{t_1}^{t_2} f(q(t), \dot{q}(t), t) \, dt$, the functional derivative is 

$$\frac{\delta F}{\delta q(t)} = \frac{\partial f}{\partial q} - \frac{d}{dt}\left(\frac{\partial f}{\partial \dot{q}}\right)$$

The Euler-Lagrange equations can be rewritten in terms of functional derivatives

$$\frac{\delta S}{\delta q_i(t)} = 0, \quad i = 1, \ldots, n ,$$

which is a compact form emphasizing that the physical path is one where the action is stationary with respect to all possible variations.

The functional derivative of the action is related to the generalized force $Q_i$ 
acting on the system $Q_i = \frac{\delta S}{\delta q_i(t)}$.
When the system follows the physical path, these generalized forces vanish, 
indicating equilibrium in the variational sense.

# Conservation Laws and Noether's Theorem

Noether's theorem establishes that each continuous symmetry of the action corresponds
to a conservation law. For time-translation symmetry, the conserved quantity is energy.
When the Lagrangian does not explicitly depend on time 
($\frac{\partial L}{\partial t} = 0$), the energy (Hamiltonian) $H$ is conserved, i.e.,

$$H = \sum_{i=1}^{n} \dot{q}_i \frac{\partial L}{\partial \dot{q}_i} - L = \text{constant} .$$

The Hamiltonian is obtained from the Lagrangian via a *Legendre transformation*

$$H(q, p, t) = \sum_{i=1}^{n} p_i \dot{q}_i - L(q, \dot{q}, t) ,$$

where the generalized momenta are $p_i = \frac{\partial L}{\partial \dot{q}_i}$.

## Hamilton's Equations of Motion

From the Hamiltonian, we can derive Hamilton's canonical equations

$$\dot{q}_i = \frac{\partial H}{\partial p_i}, \quad \dot{p}_i = -\frac{\partial H}{\partial q_i},$$
which represent an alternative formulation to the *Euler-Lagrange equations.*

## The Energy-Time Relation

For systems with time-translation symmetry, the Hamiltonian generates time evolution via

$$\frac{d}{dt}f(q, p, t) = \{f, H\} + \frac{\partial f}{\partial t},$$

where $\{f, H\}$ is the *Poisson bracket*,

$$\{f, H\} = \sum_{i=1}^{n} \left(\frac{\partial f}{\partial q_i}\frac{\partial H}{\partial p_i} - \frac{\partial f}{\partial p_i}\frac{\partial H}{\partial q_i}\right) .$$

# Group Theory and Generators in Physics

A *Lie group* $G$ is a continuous group that is also a differentiable manifold 
where the *tangent space* at the *identity element* forms a Lie algebra $\mathfrak{g}$.
For a one-parameter subgroup $g(t) \in G$, the *generator* $X \in \mathfrak{g}$ is

$$X = \left.\frac{d}{dt}g(t)\right|_{t=0} .$$

*Generators* of symmetry transformations correspond to conserved quantities in physical systems, e.g,. 

 - *Time translations*: Energy (Hamiltonian) $H$
 - *Spatial translations*: Linear momentum $\vec{P}$
 - *Rotations*: Angular momentum $\vec{L}$

Symbolically, if $U(g)$ is a unitary representation of a Lie group element $g$, and
$X$ is the corresponding generator, then $U(g) = e^{iXt}$.
In *classical mechanics*, generators form a Lie algebra under the *Poisson bracket*,
i.e., the components of angular momentum satisfy

$$\{L_i, L_j\} = \sum_{k=1}^{3} \epsilon_{ijk} L_k ,$$

where $\epsilon_{ijk}$ is the Levi-Civita symbol.
In *quantum mechanics*, the Poisson bracket corresponds to the commutator

$$[L_i, L_j] = i\hbar \sum_{k=1}^{3} \epsilon_{ijk} L_k .$$
In essence, the Poisson bracket is the realization of the classical limit of 
the quantum commutator

$$\lim_{\hbar \to 0} \frac{1}{i\hbar}[\hat{A}, \hat{B}] = \{A, B\} .$$

{\bf{Theorem}} [Noether's theorem in terms of generators] Let $G$ be a generator 
of a symmetry transformation that leaves the action invariant, then $G$ is a 
conserved quantity $\frac{dG}{dt} = 0$.

# Distribution Action and Path Integrals

In quantum mechanics, the *transition amplitude between states* is given by 
Feynman's path integral

$$\langle q_f, t_f | q_i, t_i \rangle = \int \mathcal{D}q \, e^{\frac{i}{\hbar}S[q]},$$

where $\mathcal{D}q$ represents integration over all possible paths connecting $(q_i, t_i)$ to $(q_f, t_f)$.
The probability amplitude $e^{\frac{i}{\hbar}S[q]}$ assigns a complex weight to each path, with the classical path corresponding to the stationary phase.

More rigorously, the path integral is defined as a limit

$$\int \mathcal{D}q \, e^{\frac{i}{\hbar}S[q]} = \lim_{N \to \infty} \left(\frac{m}{2\pi i \hbar \epsilon}\right)^{N/2} \int \prod_{j=1}^{N-1} dq_j \, \exp\left[\frac{i}{\hbar}\sum_{j=0}^{N-1} \epsilon L\left(\frac{q_{j+1}-q_j}{\epsilon}, \frac{q_{j+1}+q_j}{2}\right)\right] ,$$

where $\epsilon = \frac{t_f - t_i}{N}$ and $q_0 = q_i$, $q_N = q_f$.

{\bf{Definition}} [Distribution Action] The *distribution action* extends the concept
of action to encompass statistical distributions over paths. In statistical field theory, a *partition function* is

$$Z = \int \mathcal{D}\phi \, e^{-S_E[\phi]},$$

where $S_E$ is the *Euclidean action obtained via analytic continuation.*

The *distribution action* treats the action as a functional that assigns 
probabilities (or probability amplitudes) to different field configurations.

# Analytic Continuation 

Analytic continuation extends a function from a subset of its domain to a larger 
domain, preserving its analytic properties. A function $f(z)$ is analytic at a point
$z_0$ if it is differentiable in a neighborhood of $z_0$.
When a pair of analytic functions $f(z)$ and $g(z)$ agree on a set with an 
accumulation point, then $f(z) = g(z)$ throughout their domains of analyticity.

The *Wick rotation* is a specific analytic continuation from *real time* $t$ 
(event order) to *purely imaginary time* $\tau = it$. Under this transformation,
the Minkowski metric becomes Euclidean

$$ds^2 = -dt^2 + d\vec{x}^2 \to d\tau^2 + d\vec{x}^2 .$$

The (general) *action* transforms as $S = \int dt \, L \to iS_E = i\int d\tau \, L_E$, 
where $S_E$ is the *Euclidean action* and the path integral becomes

$$\int \mathcal{D}\phi \, e^{\frac{i}{\hbar}S[\phi]} \to \int \mathcal{D}\phi 
\, e^{-\frac{1}{\hbar}S_E[\phi]} .$$

This transforms *oscillatory integrals* into *exponentially decaying* integrals
that are more amenable to numerical and analytical techniques.

Under eriodic boundary conditions 
$\phi(\tau) = \phi(\tau + \beta)$, the analytic continuation $t \to i\tau$ 
connects quantum field theory to thermal field theory, where 
$\beta = \frac{1}{k_B T}$.
The *partition function* of a quantum system at temperature $T$ is

$$Z = \text{Tr}(e^{-\beta \hat{H}}) = \int \mathcal{D}\phi \, e^{-S_E[\phi]},$$
where the functional integration is performed over fields satisfying periodic boundary conditions in *imaginary time.*

Analytic continuation is useful for deriving dispersion relations that
connect the real and imaginary parts of response functions through the 
Kramers-Kronig relations

$$\text{Re}\chi(\omega) = \frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty} \frac{\text{Im}\chi(\omega')}{\omega' - \omega} d\omega' ,$$

$$\text{Im}\chi(\omega) = -\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty} \frac{\text{Re}\chi(\omega')}{\omega' - \omega} d\omega' ,$$

where $\mathcal{P}$ denotes the principal value. These relations follow from the analyticity of response functions in the upper half of the complex plane, which is a consequence of causality.

# The Action as a Central Concept

The action principle provides as a unifying framework for classical and quantum physics

 - In *classical mechanics*, it yields the equations of motion through the principle of stationary action.
 - In *quantum mechanics*, it determines probability amplitudes through the path integral.
 - In *statistical mechanics*, the Euclidean action (after analytic continuation) determines equilibrium probabilities.
 - In *spacekime*, the action of the kime-phase distribution, $\Phi$, on kime-test functions, see [TCIU Resources](https://www.socr.umich.edu/TCIU/HTMLs/Chapter6_TCIU_KimeInterpretation_of_QuantumMechanics.html). This is still *work in progress* ...

The *action of the kime phase distribution* $\Phi_t$ on *kime-test functions* $\psi$ is
defined by
$$\langle \Phi_t, \psi \rangle = \int_{-\pi}^{\pi} \psi(t e^{i\theta}) p_\Phi(\theta) \, d\theta .$$
This distribution action on test-functions generalizes measurable quantities by focusing on their functional properties rather than individual values of \( \theta \).
This analytic representation emphasizes the role of smoothing and averaging over phase, which parallels the experimental approach of repeated draws of a random variable. The former is more abstract but better suited to describe underlying kime properties, where as the latter is more useful for data-driven estimation, prediction and dynamical quantification.

The extension from *global* to *local symmetries* introduces gauge fields and leads to gauge theories. 
For a field theory with global symmetry under a transformation $\phi \to \phi + \delta\phi$, requiring local invariance under $\phi \to \phi + \delta\phi(x)$ necessitates introducing a gauge field $A_\mu(x)$ with transformation properties that compensate for the local nature of the symmetry.

The *gauge-invariant action* takes the form
$$S = \int d^4x \left[\mathcal{L}_{\text{matter}}(\phi, D_\mu\phi) - \frac{1}{4}F_{\mu\nu}F^{\mu\nu}\right] ,$$

where $D_\mu = \partial_\mu - igA_\mu$ is the *covariant derivative* and 
$F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu + ig[A_\mu, A_\nu]$ 
is the *field strength tensor.*

{\it{Geometric Interpretation}}: The action principle can be formulated in terms of differential geometry

 - The *configuration space* of a system forms a manifold $\mathcal{M}$.
 - The Lagrangian defines a *one-form* $\mathcal{L}dt$ on the tangent bundle $T\mathcal{M}$.
 - The *action* is the integral of this one-form along a path in $T\mathcal{M}$.
 - The Euler-Lagrange equations express the condition that the *path is an extremal* of this integral.

In general relativity, the *Einstein-Hilbert action* is
$S_{EH} = \frac{1}{16\pi G}\int d^4x \sqrt{-g}R ,$
where $g$ is the determinant of the metric tensor $g_{\mu\nu}$ and $R$ is the Ricci scalar.

The concepts of *action*, *functional derivatives*, *group generators*, 
*distribution actions*, analytic continuation, and energy are interconnected and their
further integration with complex-time (kime) representation may be valuable. The action principle appears as the central unifying concept to build upon.
The interplay between physical principles and mathematical structures, demonstrate
how symmetry, geometry, and variational principles come together describe the fundamental laws of nature. The main question is how to connect the statistical and data-science formulation of complex time where the repeated sampling of the 
controlled processes provides a mechanism to obtain a higher-dimensional representation
that facilitates probing the unobservable kime-phase distributions. 

# Experiments and Simulations

## Experiment 1: The Principle of Stationary Action (Harmonic Oscillator)

In this experiment, we showcase Hamilton's principle of least action using a simple harmonic oscillator. We explore how different paths connecting the same initial and
final points have different values of the action, and how the classical path 
corresponds to the path of stationary action.

To start we'll use these basic functions.

```{r functions, echo=TRUE, fig.width=7, fig.height=8, fig.align='center', warning=FALSE, message=FALSE, error=FALSE}
library(plotly)
library(purrr)

# Core functions for our simulations
calculate_lagrangian <- function(x, v, m, k) {
  # Lagrangian L = T - V = 0.5*m*v^2 - 0.5*k*x^2
  L <- 0.5*m*v^2 - 0.5*k*x^2
  return(L)
}

calculate_action <- function(x, t, m, k) {
  # Calculate velocity using central difference
  dt <- t[2] - t[1]
  v <- c(diff(x)/dt, (x[length(x)] - x[length(x)-1])/dt)
  
  # Compute Lagrangian at each point
  L <- calculate_lagrangian(x, v, m, k)
  
  # Action is integral of L over time
  S <- sum(L) * dt
  return(S)
}

# Function to generate a path with perturbations
generate_path <- function(t, omega, amplitude = 1, perturbation = NULL, 
                          pert_amplitude = 0, pert_frequency = 1) {
  x_classical <- amplitude * cos(omega * t)
  
  if (!is.null(perturbation)) {
    x <- x_classical + pert_amplitude * perturbation(pert_frequency * omega * t)
  } else {
    x <- x_classical
  }
  
  return(x)
}

# Function to calculate action for multiple paths
calculate_actions_for_paths <- function(paths_list, t, m, k) {
  # Use lapply and unlist instead of map_dbl to avoid dependency on purrr
  actions <- unlist(lapply(paths_list, function(path) calculate_action(path, t, m, k)))
  return(actions)
}
```

The first interactive simulation allows how various parameters affect the action
and observe the principle of stationary action in practice.
The simulation above demonstrates several key aspects of the action principle:

 - *Principle of Stationary Action*: The classical path (*green*) corresponds to the
path with the minimum action value among all paths connecting the same endpoints.

 - *Variational Principle*: As we move away from the classical path (with increasing perturbation amplitude), the action increases, demonstrating that the classical path makes the action stationary.
 - *Physical Meaning*: The color gradient in the perturbation paths visualizes how far each path deviates from optimality in terms of action. Paths closer to the classical one (bluer) have action values closer to the minimum.

```{r sim1, fig.width=7, fig.height=8, fig.align='center', warning=FALSE, message=FALSE, error=FALSE}
# Visualization of a simple harmonic oscillator action

# # Parameters
# m <- 1  # mass
# k <- 1  # spring constant
# omega <- sqrt(k/m)  # angular frequency
# T <- 2*pi/omega  # period
# 
# # Generate data for different paths
# t <- seq(0, T, length.out=100)
# x_classical <- cos(omega*t)  # classical path
# x_perturbed1 <- cos(omega*t) + 0.2*sin(2*omega*t)  # perturbed path 1
# x_perturbed2 <- cos(omega*t) + 0.4*sin(3*omega*t)  # perturbed path 2
# 
# # Calculate actions
# calculate_action <- function(x, t) {
#   # Approximate velocities
#   v <- c(diff(x)/diff(t), 0)
#   
#   # Lagrangian L = T - V = 0.5*m*v^2 - 0.5*k*x^2
#   L <- 0.5*m*v^2 - 0.5*k*x^2
#   
#   # Action is integral of L over time
#   S <- sum(L) * (t[2] - t[1])
#   return(S)
# }
# 
# actions <- c(
#   calculate_action(x_classical, t),
#   calculate_action(x_perturbed1, t),
#   calculate_action(x_perturbed2, t)
# )
# 
# # Create data frame for plotting
# df <- data.frame(
#   time = rep(t, 3),
#   position = c(x_classical, x_perturbed1, x_perturbed2),
#   path = rep(c("Classical", "Perturbed 1", "Perturbed 2"), each=length(t))
# )
# 
# # Create plot
# ggplot(df, aes(x=time, y=position, color=path)) +
#   geom_line() +
#   labs(title="Different Paths of a Harmonic Oscillator",
#        subtitle=paste("Classical path has stationary action S =", round(actions[1], 2)),
#        x="Time", y="Position") +
#   theme_minimal()


# Simulation parameters - user adjustable
m <- 1.0  # mass
k <- 2.0  # spring constant
omega <- sqrt(k/m)  # angular frequency
T <- 2*pi/omega  # period
num_points <- 200  # resolution of our simulation
num_paths <- 10  # number of perturbed paths to generate

# Time points for simulation
t <- seq(0, T, length.out = num_points)

# Generate classical path
x_classical <- generate_path(t, omega)

# Generate perturbed paths with varying amplitudes
perturbation_amplitudes <- seq(0.05, 0.5, length.out = num_paths)
paths_list <- list()
paths_list[[1]] <- x_classical  # First path is classical

for (i in 1:num_paths) {
  # Add sine perturbation with varying amplitude and frequency
  paths_list[[i+1]] <- generate_path(
    t, omega, 
    perturbation = sin,
    pert_amplitude = perturbation_amplitudes[i],
    pert_frequency = 2
  )
}

# Calculate action for each path
path_names <- c("Classical", paste("Perturbed", 1:num_paths))
actions <- calculate_actions_for_paths(paths_list, t, m, k)

# Create data frame for plotting trajectories
trajectory_df <- map_dfr(1:(num_paths+1), function(i) {
  data.frame(
    time = t,
    position = paths_list[[i]],
    path = path_names[i],
    action = actions[i]
  )
})

# Create the trajectory plot
trajectory_plot <- plot_ly(height = 500) %>%
  add_trace(data = trajectory_df %>% filter(path == "Classical"),
            x = ~time, y = ~position, type = 'scatter', mode = 'lines',
            name = 'Classical Path',
            line = list(color = 'green', width = 3)) %>%
  layout(title = "Different Paths of a Harmonic Oscillator",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Position"),
         hovermode = "closest")

# Add perturbed paths with color gradient based on action value
perturbed_only <- trajectory_df %>% filter(path != "Classical")
action_range <- range(perturbed_only$action)
norm_actions <- (perturbed_only$action - action_range[1]) / (action_range[2] - action_range[1])

for (i in 1:num_paths) {
  path_data <- perturbed_only %>% filter(path == paste("Perturbed", i))
  # Color from blue (close to classical) to red (far from classical)
  color_val <- rgb(norm_actions[i*num_points], 0, 1-norm_actions[i*num_points])
  
  trajectory_plot <- trajectory_plot %>%
    add_trace(data = path_data, x = ~time, y = ~position, 
              type = 'scatter', mode = 'lines',
              name = paste0(path_data$path[1], " (S = ", round(path_data$action[1], 2), ")"),
              line = list(color = color_val, width = 1.5))
}

# Plot the action values for each path
action_df <- data.frame(
  path = factor(path_names, levels = path_names),
  action = actions
)

action_plot <- plot_ly(data = action_df, x = ~path, y = ~action, type = 'bar',
                       marker = list(color = c('green', colorRampPalette(c('blue', 'red'))(num_paths)))) %>%
  layout(title = "Action Value for Each Path",
         xaxis = list(title = "Path"),
         yaxis = list(title = "Action Value"),
         showlegend = FALSE)

# Display the plots
subplot(trajectory_plot, action_plot, nrows = 2, heights = c(0.7, 0.3)) %>%
  layout(title = "Principle of Stationary Action Demonstration",
         annotations = list(
           list(
             x = 0.5,
             y = 1.05,
             text = "The classical path minimizes the action functional",
             showarrow = FALSE,
             xref = "paper",
             yref = "paper"
           )
         ))
```

## Experiment 2: Group Generators and Conservation Laws

This experiment shows Noether's theorem, connecting symmetries to conservation laws
and shows how time-translation symmetry leads to energy conservation, and how spatial
symmetries relate to momentum conservation.

 - Harmonic Oscillator (Time-Translation Symmetry):
   - The Lagrangian has no explicit time dependence
   - The generator of time translations is the Hamiltonian/energy
   - Energy is perfectly conserved as shown by the flat red line


 - Free Particle (Spatial-Translation Symmetry):
   - The system is invariant under spatial translations
   - The generator of spatial translations is momentum
   - Momentum is conserved (constant purple line)


 - Particle in Gravity (Broken Time-Translation Symmetry):
   - The potential energy depends on position, which changes with time
   - Time-translation symmetry is broken
   - Total energy is not conserved, as seen in the varying red line.


```{r sim2, warning=FALSE, message=FALSE, error=FALSE}
# Function to simulate a physical system with symmetries
simulate_system_with_symmetries <- function(t, initial_conditions, 
                                           params = list(m = 1, k = 1, g = 9.8),
                                           system_type = "harmonic") {
  
  # Unpack initial conditions
  x0 <- initial_conditions$x0
  v0 <- initial_conditions$v0
  
  # Unpack parameters
  m <- params$m
  k <- params$k
  g <- params$g
  
  # Define different system types
  if (system_type == "harmonic") {
    # Harmonic oscillator (time-translation symmetry)
    x <- x0 * cos(sqrt(k/m) * t) + v0/sqrt(k/m) * sin(sqrt(k/m) * t)
    v <- -x0 * sqrt(k/m) * sin(sqrt(k/m) * t) + v0 * cos(sqrt(k/m) * t)
    
    # Potential and kinetic energy
    V <- 0.5 * k * x^2
    T <- 0.5 * m * v^2
    
    # Total energy (should be conserved due to time symmetry)
    E <- T + V
    
    return(data.frame(t = t, x = x, v = v, T = T, V = V, E = E))
    
  } else if (system_type == "free_particle") {
    # Free particle (spatial translation symmetry)
    x <- x0 + v0 * t
    v <- rep(v0, length(t))
    
    # Momentum (should be conserved due to spatial symmetry)
    p <- m * v
    
    # Kinetic energy
    T <- 0.5 * m * v^2
    
    return(data.frame(t = t, x = x, v = v, p = p, T = T))
    
  } else if (system_type == "gravity") {
    # Particle in gravity (broken time symmetry)
    x <- x0 + v0 * t - 0.5 * g * t^2
    v <- v0 - g * t
    
    # Potential and kinetic energy
    V <- m * g * x
    T <- 0.5 * m * v^2
    
    # Total energy (not conserved due to time-dependent potential)
    E <- T + V
    
    return(data.frame(t = t, x = x, v = v, T = T, V = V, E = E))
  }
}

# Time points for simulation
t <- seq(0, 10, length.out = 200)

# Initial conditions
initial_conditions <- list(x0 = 1, v0 = 0)

# Simulate three different systems
harmonic_data <- simulate_system_with_symmetries(t, initial_conditions, 
                                                params = list(m = 1, k = 2),
                                                system_type = "harmonic")

free_particle_data <- simulate_system_with_symmetries(t, initial_conditions, 
                                                     params = list(m = 1),
                                                     system_type = "free_particle")

gravity_data <- simulate_system_with_symmetries(t, initial_conditions, 
                                               params = list(m = 1, g = 0.5),
                                               system_type = "gravity")

# Interactive plots for each system
harmonic_plot <- plot_ly(harmonic_data, height = 400) %>%
  add_trace(x = ~t, y = ~E, type = 'scatter', mode = 'lines', name = 'Total Energy',
            line = list(color = 'red', width = 3)) %>%
  add_trace(x = ~t, y = ~T, type = 'scatter', mode = 'lines', name = 'Kinetic Energy',
            line = list(color = 'blue')) %>%
  add_trace(x = ~t, y = ~V, type = 'scatter', mode = 'lines', name = 'Potential Energy',
            line = list(color = 'green')) %>%
  layout(title = "Harmonic Oscillator (Time-Translation Symmetry)",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Energy"),
         annotations = list(
           list(
             x = 5,
             y = max(harmonic_data$E) * 1.1,
             text = "Energy is conserved due to time-translation symmetry",
             showarrow = FALSE
           )
         ))

free_particle_plot <- plot_ly(free_particle_data, height = 400) %>%
  add_trace(x = ~t, y = ~p, type = 'scatter', mode = 'lines', name = 'Momentum',
            line = list(color = 'purple', width = 3)) %>%
  add_trace(x = ~t, y = ~T, type = 'scatter', mode = 'lines', name = 'Kinetic Energy',
            line = list(color = 'blue')) %>%
  layout(title = "Free Particle (Spatial-Translation Symmetry)",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Momentum/Energy"),
         annotations = list(
           list(
             x = 5,
             y = max(free_particle_data$p) * 1.1,
             text = "Momentum is conserved due to spatial-translation symmetry",
             showarrow = FALSE
           )
         ))

gravity_plot <- plot_ly(gravity_data, height = 400) %>%
  add_trace(x = ~t, y = ~E, type = 'scatter', mode = 'lines', name = 'Total Energy',
            line = list(color = 'red', width = 3)) %>%
  add_trace(x = ~t, y = ~T, type = 'scatter', mode = 'lines', name = 'Kinetic Energy',
            line = list(color = 'blue')) %>%
  add_trace(x = ~t, y = ~V, type = 'scatter', mode = 'lines', name = 'Potential Energy',
            line = list(color = 'green')) %>%
  layout(title = "Particle in Gravity (Broken Time-Translation Symmetry)",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Energy"),
         annotations = list(
           list(
             x = 5,
             y = min(gravity_data$E) * 0.9,
             text = "Energy is not conserved due to broken time-translation symmetry",
             showarrow = FALSE
           )
         ))

# Display the plots
subplot(harmonic_plot, free_particle_plot, gravity_plot, nrows = 3) %>%
  layout(title = "Noether's Theorem: Symmetries and Conservation Laws")
```


## Experiment 3: Path Integrals and Quantum Tunneling

This experiment visualizes the quantum path integral approach by simulating quantum
tunneling through a potential barrier and shows the classical path as just one
contribution to the quantum amplitude, and how multiple paths contribute to quantum
phenomena.

The path integral formulation of quantum mechanics shows:

 - *Multiple Paths*: Unlike classical mechanics where only one path contributes, quantum mechanics considers all possible paths between initial and final states.
 - *Complex Amplitudes*: Each path contributes with a complex amplitude $e^{iS/\hbar}$. The distribution of these amplitudes in the complex plane (bottom plot) shows how paths interfere.
 - *Quantum Tunneling*: Paths that cross the potential barrier (middle plot) have non-zero probability amplitudes, demonstrating quantum tunneling - a phenomenon impossible in classical mechanics.
 - Role of $\hbar$: The Planck's constant, $\hbar$, controls the "quantumness" of the system. Smaller values make the system more classical, larger values enhance quantum effects.


```{r sim3, warning=FALSE, message=FALSE, error=FALSE}
# Function to calculate potential
calculate_potential <- function(x, barrier_height = 1, barrier_width = 1) {
  # Double well potential with barrier
  potential <- numeric(length(x))
  
  for (i in 1:length(x)) {
    if (abs(x[i]) < barrier_width/2) {
      potential[i] <- barrier_height
    } else {
      potential[i] <- 0.1 * (x[i]^2 - barrier_width^2)^2
    }
  }
  
  return(potential)
}

# Function to generate random paths for path integral
generate_random_paths <- function(x_start, x_end, n_steps, n_paths, max_deviation = 0.5) {
  paths <- matrix(0, nrow = n_paths, ncol = n_steps)
  
  # Linear path from start to end
  straight_path <- seq(x_start, x_end, length.out = n_steps)
  
  for (i in 1:n_paths) {
    # Start with straight path
    current_path <- straight_path
    
    # Add random deviations, keeping endpoints fixed
    deviations <- c(0, rnorm(n_steps-2, 0, max_deviation), 0)
    paths[i,] <- current_path + deviations
  }
  
  return(paths)
}

# Calculate quantum action (with ħ term)
calculate_quantum_action <- function(path, dt, m = 1, potential_fn, hbar = 1) {
  n <- length(path)
  dx <- diff(path)
  v <- dx/dt
  
  # Kinetic energy
  T <- 0.5 * m * v^2
  
  # Potential energy (average of adjacent points)
  V <- potential_fn(path[-1]) # potential at midpoints
  
  # Quantum action includes i*ħ
  S <- sum((T - V) * dt)
  
  return(S)
}

# Calculate quantum amplitude (complex)
calculate_amplitude <- function(action, hbar = 1) {
  return(exp(1i * action / hbar))
}

# Simulate quantum tunneling
simulate_tunneling <- function(x_start = -2, x_end = 2, n_steps = 50, n_paths = 100, 
                              barrier_params = list(height = 1.5, width = 1.5),
                              hbar = 0.5) {
  # Time step
  dt <- 0.1
  
  # Generate random paths
  paths <- generate_random_paths(x_start, x_end, n_steps, n_paths)
  
  # Create potential function with barrier
  potential_fn <- function(x) calculate_potential(x, barrier_params$height, barrier_params$width)
  
  # Position grid for potential plot
  x_grid <- seq(-3, 3, length.out = 200)
  potential_grid <- potential_fn(x_grid)
  
  # Calculate action and amplitude for each path
  results <- data.frame(
    path_id = 1:n_paths,
    action = numeric(n_paths),
    amplitude_real = numeric(n_paths),
    amplitude_imag = numeric(n_paths),
    amplitude_abs = numeric(n_paths)
  )
  
  path_data <- list()
  
  for (i in 1:n_paths) {
    # Calculate action for this path
    action_val <- calculate_quantum_action(paths[i,], dt, potential_fn = potential_fn, hbar = hbar)
    amplitude <- calculate_amplitude(action_val, hbar)
    
    results$action[i] <- action_val
    results$amplitude_real[i] <- Re(amplitude)
    results$amplitude_imag[i] <- Im(amplitude)
    results$amplitude_abs[i] <- Mod(amplitude)
    
    # Store path data for plotting
    path_data[[i]] <- data.frame(
      x = paths[i,],
      t = seq(0, (n_steps-1)*dt, by = dt),
      path_id = i,
      action = action_val,
      amplitude_abs = Mod(amplitude)
    )
  }
  
  # Combine path data
  all_paths <- do.call(rbind, path_data)
  
  # Create potential plot
  potential_plot <- plot_ly(height = 300) %>%
    add_trace(x = x_grid, y = potential_grid, type = 'scatter', mode = 'lines',
              name = 'Potential', line = list(color = 'black', width = 2)) %>%
    layout(title = "Potential Barrier",
           xaxis = list(title = "Position"),
           yaxis = list(title = "Potential Energy"))
  
  # Plot a subset of paths colored by amplitude
  paths_to_plot <- sample(1:n_paths, min(30, n_paths))
  
  path_plot <- plot_ly(height = 300) %>%
    layout(title = "Sample Quantum Paths",
           xaxis = list(title = "Time"),
           yaxis = list(title = "Position"))
  
  for (i in paths_to_plot) {
    path_subset <- all_paths %>% filter(path_id == i)
    
    # Normalize amplitude for color
    norm_amplitude <- (path_subset$amplitude_abs[1] - min(results$amplitude_abs)) / 
                      (max(results$amplitude_abs) - min(results$amplitude_abs))
    
    # Color from blue (low amplitude) to red (high amplitude)
    color_val <- rgb(norm_amplitude, 0, 1-norm_amplitude)
    
    path_plot <- path_plot %>%
      add_trace(data = path_subset, x = ~t, y = ~x, type = 'scatter', mode = 'lines',
                line = list(color = color_val),
                showlegend = FALSE)
  }
  
  # Create amplitude distribution plot
  amplitude_plot <- plot_ly(data = results, height = 300) %>%
    add_trace(x = ~amplitude_real, y = ~amplitude_imag, type = 'scatter', mode = 'markers',
              marker = list(size = 5, color = ~amplitude_abs, 
                            colorscale = 'Viridis', showscale = TRUE)) %>%
    layout(title = "Distribution of Path Amplitudes in Complex Plane",
           xaxis = list(title = "Re(amplitude)"),
           yaxis = list(title = "Im(amplitude)"))
  
  # Return all plots and data
  return(list(
    potential_plot = potential_plot,
    path_plot = path_plot,
    amplitude_plot = amplitude_plot,
    results = results,
    all_paths = all_paths
  ))
}

# Run the quantum tunneling simulation
tunneling_sim <- simulate_tunneling(hbar = 0.5)

# Display the plots
subplot(tunneling_sim$potential_plot, 
        tunneling_sim$path_plot, 
        tunneling_sim$amplitude_plot, 
        nrows = 3) %>%
  layout(title = "Quantum Tunneling via Path Integral Approach",
         annotations = list(
           list(
             x = 0.5,
             y = 1.05,
             text = paste("Quantum paths contribute with complex amplitudes e^(iS/ħ),",
                          "allowing tunneling through classically forbidden regions"),
             showarrow = FALSE,
             xref = "paper",
             yref = "paper"
           )
         ))
```


## Experiment 4: Analytic Continuation and Euclidean Path Integral

This experiment demonstrates how analytic continuation transforms oscillatory 
path integrals into more manageable Euclidean path integrals, facilitating both
numerical calculations and connections to statistical mechanics.

```{r sim4, warning=FALSE, message=FALSE, error=FALSE}
# Function to simulate Minkowski and Euclidean path integrals
simulate_analytic_continuation <- function(n_paths = 1000, n_steps = 20, hbar = 1) {
  
  # Time parameters
  t_max <- 2
  t <- seq(0, t_max, length.out = n_steps)
  dt <- t[2] - t[1]
  
  # Generate paths for simple harmonic oscillator
  m <- 1
  k <- 1
  omega <- sqrt(k/m)
  
  # Classical path
  x_classical <- cos(omega * t)
  
  # Generate random paths with fixed endpoints
  x_start <- x_classical[1]
  x_end <- x_classical[n_steps]
  
  paths <- matrix(0, nrow = n_paths, ncol = n_steps)
  paths[,1] <- x_start
  paths[,n_steps] <- x_end
  
  # Random fluctuations around straight line
  for (i in 1:n_paths) {
    # Generate random intermediate points
    for (j in 2:(n_steps-1)) {
      # Linear interpolation plus random deviation
      paths[i,j] <- x_start + (j-1)/(n_steps-1) * (x_end - x_start) + 
                     rnorm(1, 0, 0.3 * sqrt(dt))
    }
  }
  
  # Calculate actions and amplitudes for both Minkowski and Euclidean cases
  results <- data.frame(
    path_id = 1:n_paths,
    minkowski_action = numeric(n_paths),
    minkowski_amplitude_re = numeric(n_paths),
    minkowski_amplitude_im = numeric(n_paths),
    euclidean_action = numeric(n_paths),
    euclidean_weight = numeric(n_paths)
  )
  
  for (i in 1:n_paths) {
    # Current path
    path <- paths[i,]
    
    # Calculate velocities
    v <- c(diff(path)/dt, 0)
    
    # Lagrangian L = T - V
    L <- 0.5 * m * v[1:(n_steps-1)]^2 - 0.5 * k * path[1:(n_steps-1)]^2
    
    # Minkowski action (real time)
    S_M <- sum(L * dt)
    results$minkowski_action[i] <- S_M
    
    # Minkowski amplitude e^(iS/ħ)
    amplitude <- exp(1i * S_M / hbar)
    results$minkowski_amplitude_re[i] <- Re(amplitude)
    results$minkowski_amplitude_im[i] <- Im(amplitude)
    
    # Euclidean action (imaginary time)
    # Mathematically: t → iτ, S → iS_E
    # Lagrangian becomes L_E = T + V (sign change)
    L_E <- 0.5 * m * v[1:(n_steps-1)]^2 + 0.5 * k * path[1:(n_steps-1)]^2
    S_E <- sum(L_E * dt)
    results$euclidean_action[i] <- S_E
    
    # Euclidean weight e^(-S_E/ħ)
    results$euclidean_weight[i] <- exp(-S_E / hbar)
  }
  
  # Normalize Euclidean weights for visualization
  results$euclidean_weight_norm <- results$euclidean_weight / max(results$euclidean_weight)
  
  # Prepare path data for plotting
  path_data <- list()
  selected_paths <- sample(1:n_paths, 30)  # Select random paths to plot
  
  for (i in selected_paths) {
    path_data[[length(path_data) + 1]] <- data.frame(
      t = t,
      x = paths[i,],
      path_id = i,
      minkowski_action = results$minkowski_action[i],
      euclidean_weight = results$euclidean_weight[i],
      euclidean_weight_norm = results$euclidean_weight_norm[i]
    )
  }
  
  # Add classical path
  path_data[[length(path_data) + 1]] <- data.frame(
    t = t,
    x = x_classical,
    path_id = 0,
    minkowski_action = NA,
    euclidean_weight = NA,
    euclidean_weight_norm = NA
  )
  
  all_paths <- bind_rows(path_data)
  
  # Return all the computed data
  return(list(
    paths = paths,
    results = results,
    path_data = all_paths,
    classical_path = x_classical,
    t = t
  ))
}

# Load required libraries
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(viridis)
library(grid)
library(gridExtra)

# Run the simulation
set.seed(42)
sim_results <- simulate_analytic_continuation(n_paths = 1000, n_steps = 50, hbar = 0.5)

# Plot a selection of paths with color representing Euclidean weight
p1 <- ggplot(filter(sim_results$path_data, path_id != 0)) +
  geom_line(aes(x = t, y = x, group = path_id, color = euclidean_weight_norm), alpha = 0.7) +
  geom_line(data = filter(sim_results$path_data, path_id == 0), 
            aes(x = t, y = x), color = "red", size = 1, linetype = "dashed") +
  scale_color_viridis_c(name = "Euclidean\nWeight", direction = -1) +
  labs(title = "Quantum Paths for Harmonic Oscillator",
       subtitle = "Colored by Euclidean weight (higher = darker)",
       x = "Time",
       y = "Position") +
  theme_minimal() +
  theme(legend.position = "right")

# Distribution of Minkowski actions and amplitudes
p2 <- ggplot(sim_results$results) +
  geom_point(aes(x = minkowski_amplitude_re, y = minkowski_amplitude_im, color = minkowski_action),
             alpha = 0.5) +
  scale_color_viridis_c(name = "Minkowski\nAction") +
  coord_fixed() +
  labs(title = "Minkowski Path Integral",
       subtitle = "Complex amplitude distribution",
       x = "Re(e^(iS/ħ))",
       y = "Im(e^(iS/ħ))") +
  theme_minimal() +
  theme(legend.position = "right")

# Distribution of Euclidean actions and weights
p3 <- ggplot(sim_results$results) +
  geom_histogram(aes(x = euclidean_action, weight = euclidean_weight), bins = 30, fill = "steelblue") +
  labs(title = "Euclidean Path Integral",
       subtitle = "Action distribution weighted by e^(-S/ħ)",
       x = "Euclidean Action",
       y = "Weighted Frequency") +
  theme_minimal()

# Relationship between Minkowski and Euclidean actions
p4 <- ggplot(sim_results$results) +
  geom_point(aes(x = euclidean_action, y = minkowski_action, color = euclidean_weight_norm),
             alpha = 0.5) +
  scale_color_viridis_c(name = "Euclidean\nWeight", direction = -1) +
  labs(title = "Relationship Between Actions",
       subtitle = "Minkowski vs. Euclidean",
       x = "Euclidean Action",
       y = "Minkowski Action") +
  theme_minimal() +
  theme(legend.position = "right")

# Combine plots
combined_plot <- (p1 + p2) / (p3 + p4) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "right")

print(combined_plot)

# Additional analysis: Compute expectation values
compute_expectation_values <- function(sim_results) {
  # Extract relevant data
  paths <- sim_results$paths
  weights <- sim_results$results$euclidean_weight
  t <- sim_results$t
  n_steps <- length(t)
  
  # Normalize weights
  weights_norm <- weights / sum(weights)
  
  # Compute position expectation value at each time
  pos_expectation <- numeric(n_steps)
  pos2_expectation <- numeric(n_steps)
  
  for (j in 1:n_steps) {
    pos_expectation[j] <- sum(paths[, j] * weights_norm)
    pos2_expectation[j] <- sum(paths[, j]^2 * weights_norm)
  }
  
  # Compute variance
  pos_variance <- pos2_expectation - pos_expectation^2
  
  # Return expectation values
  return(data.frame(
    t = t,
    position = pos_expectation,
    position_squared = pos2_expectation,
    position_variance = pos_variance
  ))
}

# Compute expectation values
expectation_values <- compute_expectation_values(sim_results)

# Plot expectation values
p5 <- ggplot() +
  geom_line(data = expectation_values, aes(x = t, y = position), color = "blue", size = 1) +
  geom_ribbon(data = expectation_values, 
              aes(x = t, ymin = position - sqrt(position_variance), 
                  ymax = position + sqrt(position_variance)), 
              alpha = 0.3, fill = "blue") +
  geom_line(data = data.frame(t = sim_results$t, x = sim_results$classical_path),
            aes(x = t, y = x), color = "red", linetype = "dashed") +
  labs(title = "Quantum vs. Classical Position",
       subtitle = "Blue = quantum expectation with uncertainty, Red = classical path",
       x = "Time",
       y = "Position") +
  theme_minimal()

print(p5)

# Create interpretive visualization comparing Minkowski and Euclidean formulations
# First, create demonstration data
t_values <- seq(0, 2*pi, length.out = 100)
minkowski_wave <- data.frame(
  t = t_values,
  real = cos(t_values),
  imag = sin(t_values),
  type = "Minkowski"
)

euclidean_decay <- data.frame(
  t = t_values,
  real = exp(-t_values/2),
  imag = 0,
  type = "Euclidean"
)

combined_data <- bind_rows(minkowski_wave, euclidean_decay)

p6 <- ggplot(combined_data, aes(x = t, color = type)) +
  geom_line(aes(y = real, linetype = "Real Part"), size = 1) +
  geom_line(data = filter(combined_data, type == "Minkowski"),
            aes(y = imag, linetype = "Imaginary Part"), size = 1) +
  scale_color_manual(values = c("Minkowski" = "blue", "Euclidean" = "darkred")) +
  scale_linetype_manual(values = c("Real Part" = "solid", "Imaginary Part" = "dashed")) +
  labs(title = "Analytic Continuation: Minkowski vs. Euclidean Formulations",
       subtitle = "Oscillatory behavior transforms to exponential decay",
       x = "Time (t) / Euclidean time (τ)",
       y = "Amplitude / Weight") +
  theme_minimal() +
  theme(legend.title = element_blank())

print(p6)

# Final combined analysis plot
final_plot <- p5 / p6 + plot_layout(heights = c(3, 2))
print(final_plot)

# Summary of findings
cat("Summary of Analytic Continuation Experiment:\n\n")
cat("1. Minkowski path integral: Oscillatory phase factors lead to highly fluctuating sums\n")
cat("2. Euclidean path integral: Exponential damping factors lead to well-behaved sums\n")
cat("3. Classical path corresponds closely to the path of maximum weight in Euclidean formulation\n")
cat("4. Quantum uncertainty is visible in the spread around expectation values\n")
cat("5. The Wick rotation (t → iτ) transforms oscillatory behavior to exponential decay\n")
```

This final experiment demonstrates the power of analytic continuation in quantum
mechanics, particularly when working with path integrals. The key insights include:

 1. *Oscillatory vs. Damping Behavior*: In the Minkowski formulation, the phase factor $e^{-S_E/\hbar}$ oscillates rapidly, making numerical computation challenging. The Euclidean formulation, with its weight factor $e^{-S_E/\hbar}$, provides exponential damping that emphasizes paths near the classical solution.

 2. *Statistical Connection*: The Euclidean path integral has a direct interpretation as a statistical mechanical partition function. This connection allows techniques from statistical mechanics to be applied to quantum field theory.

 3. *Numerical Advantages*: The Euclidean approach vastly improves numerical stability and convergence in path integral calculations, as shown in the weighted histogram of actions.

 4. *Visualization of Quantum Behavior*: The plots illustrate how quantum paths distribute around the classical solution, with the Euclidean weighting providing a natural importance sampling for the most relevant paths.

 5. *Uncertainty Principle*: The spread in position values around the expectation demonstrates the quantum uncertainty inherent in the system.

The technique of analytic continuation via *Wick rotation* ($t \to i\tau$) is a cornerstone of modern quantum field theory and provides one of the most important connections between quantum mechanics and statistical physics.

# Hamiltonian-Jacobi Equation for the Action

Next we'll consider the *Lagrangian approach* based on the Poincaré one-form,
the *action*, $\delta S=0$ and the *Hamiltonian approach* based on the Poincaré two-form yielding the weak (local) Noether invariants.

The Hamilton-Jacobi (HJ) equation provides a reformulation of classical mechanics
that bridges the gap between Hamiltonian mechanics and quantum mechanics. 
HJ equations focuse on the action $S$ as the fundamental quantity.
Starting from Hamilton's principal function $S(q,t)$, defined as the action 
along a classical trajectory ending at position $q$ at time $t$, the HJ equation
emerges from a canonical transformation that trivializes the equations of motion.
Symbolically, the Hamilton-Jacobi equation is

$$\frac{\partial S}{\partial t} + H\left(q, \frac{\partial S}{\partial q}, t\right) =0.$$

where $H$ is the *Hamiltonian* of the system, $q$ represents *generalized coordinates*,
and $\frac{\partial S}{\partial q}$ replaces the *canonical momentum* $p$.

This equation arises by defining a *type-2 generating function* $F_2(q,P,t) = S(q,t)$
for a canonical transformation $(q,p) \rightarrow (Q,P)$ that makes the new 
Hamiltonian $K(Q,P,t) = 0$, resulting in trivial dynamics where $Q$ and $P$ are constants.

The Hamilton-Jacobi equation has several properties

 1. *Complete Integral*: A complete solution $S(q, \alpha, t)$, where $\alpha$ are integration constants, allows us to find all possible trajectories of the system.

 2. *Action as a Generator*: The *action* $S$ serves as a generating function for canonical transformations that simplify the equations of motion.

 3. *Connection to Wave Mechanics*: The HJ equation is the classical limit of the *Schrödinger equation* when $\hbar \rightarrow 0$. If we substitute $\Psi = e^{iS/\hbar}$ into the Schrödinger equation, the leading order term yields the HJ equation as $\hbar \rightarrow 0$.

 4. *Eikonal Equation*: In optics, the HJ equation reduces to the eikonal equation, connecting mechanics with geometric optics $\left(\nabla S\right)^2 = n^2(q)$, where $n$ is the refractive index.

Below is a simulation showing the Hamilton-Jacobi solutions.

```{r hj-equation, echo=TRUE, warning=FALSE, message=FALSE, error=FALSE}
library(plotly)

# Hamilton-Jacobi solution for a free particle
generate_HJ_free_particle <- function(x_range, t_range, resolution=50, m=1) {
  # Create grid
  x <- seq(from=x_range[1], to=x_range[2], length.out=resolution)
  t <- seq(from=t_range[1], to=t_range[2], length.out=resolution)
  grid <- expand.grid(x=x, t=t)
  
  # Calculate action S(x,t) for multiple trajectories
  S_values <- matrix(0, nrow=resolution, ncol=resolution)
  
  # Different momenta (integration constants)
  momenta <- c(0.5, 1, 1.5, 2, 2.5)
  colors <- c('blue', 'red', 'green', 'purple', 'orange')
  
  # Calculate trajectories and action surfaces
  trajectory_data <- list()
  
  for(i in seq_along(momenta)) {
    p <- momenta[i]
    
    # For each momentum p, compute action S(x,t) = px - p²t/(2m)
    for(ix in 1:resolution) {
      for(it in 1:resolution) {
        x_val <- x[ix]
        t_val <- t[it]
        S_values[ix, it] <- p*x_val - (p^2)/(2*m)*t_val
      }
    }
    
    # Store action surface
    surface_data <- expand.grid(x=x, t=t)
    surface_data$S <- c(S_values)
    surface_data$momentum <- p
    
    # Generate trajectory for this momentum
    traj <- data.frame(
      t = t,
      x = p*t/m,  # Free particle trajectory x(t) = p*t/m + x₀ (setting x₀=0)
      p = p
    )
    
    trajectory_data[[i]] <- list(
      surface = surface_data,
      trajectory = traj,
      color = colors[i]
    )
  }
  
  return(list(
    x = x,
    t = t,
    trajectories = trajectory_data
  ))
}

# Generate data
hj_data <- generate_HJ_free_particle(
  x_range = c(0, 10), 
  t_range = c(0, 5),
  resolution = 40
)

# Create 3D surface plot of action
action_plot <- plot_ly(height=600)

# Add action surfaces for different momenta
for(i in seq_along(hj_data$trajectories)) {
  traj_data <- hj_data$trajectories[[i]]
  
  # Convert to matrix form for surface plot
  S_matrix <- matrix(NA, nrow=length(hj_data$x), ncol=length(hj_data$t))
  
  # Fill in the matrix
  for(ix in 1:length(hj_data$x)) {
    for(it in 1:length(hj_data$t)) {
      x_val <- hj_data$x[ix]
      t_val <- hj_data$t[it]
      p <- unique(traj_data$surface$momentum)
      S_matrix[ix, it] <- p*x_val - (p^2)/(2)*t_val
    }
  }
  
  # Add surface
  action_plot <- action_plot %>%
    add_surface(
      z = S_matrix,
      x = hj_data$t,
      y = hj_data$x,
      opacity = 0.7,
      colorscale = list(c(0, 1), c(traj_data$color, traj_data$color)),
      showscale = FALSE,
      name = paste("p =", unique(traj_data$surface$momentum))
    )
}

# Add trajectories on the action surface
for(i in seq_along(hj_data$trajectories)) {
  traj_data <- hj_data$trajectories[[i]]
  traj <- traj_data$trajectory
  p <- unique(traj$p)
  
  # Calculate S values along trajectory
  S_traj <- p*traj$x - (p^2)/(2)*traj$t
  
  action_plot <- action_plot %>%
    add_trace(
      x = traj$t,
      y = traj$x,
      z = S_traj,
      type = 'scatter3d',
      mode = 'lines',
      line = list(color='black', width=5),
      name = paste("Trajectory p =", p)
    )
}

# Final layout adjustments
action_plot %>%
  layout(
    title = "Hamilton-Jacobi Action Surfaces for Free Particle",
    scene = list(
      xaxis = list(title = "Time"),
      yaxis = list(title = "Position"),
      zaxis = list(title = "Action S(x,t)"),
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1.2))
    )
  )
```

This Hamilton-Jacobi equation simulation shows

 1. Colored surfaces representing solutions to the HJ equation for a free particle with a *different momentum* (integration constant).

 2. The black lines show the *classical trajectories*, which follow paths of constant momentum along the action surfaces.

 3. The *gradient of the action* $\nabla S$ gives the momentum at each point, while the time derivative $\partial S/\partial t$ relates to the energy.

 4. The complete integral of the HJ equation gives us a family of solutions parameterized by integration constants (different momenta in this case).

## Lagrangian Approach and the Poincaré One-Form

The Lagrangian formulation of mechanics is elegantly recast in terms of 
differential geometry through the Poincaré one-form, providing a deep geometric
understanding of the principle of least action.

The *Poincaré one-form* $\theta$ on the cotangent bundle $T^*Q$ of the 
configuration manifold $Q$ is given in local coordinates by $\theta = p_i dq^i$,
where $p_i$ are the *canonical momenta* and $q^i$ are the *generalized coordinates.*
The *Lagrangian* $L(q, \dot{q}, t)$ defines a map from the tangent bundle $TQ$ 
to the reals. The *action* $S$ is the integral of the Lagrangian over a path
$\gamma$ in the extended configuration space
$S[\gamma] = \int_{\gamma} L(q, \dot{q}, t) \, dt$.

In terms of the Poincaré one-form, the action can be written as
$$S[\gamma] = \int_{\gamma} \left( \theta - H \, dt \right) ,$$
where $H$ is the Hamiltonian.

The *principle of stationary action* states that the physical path $\gamma$ makes 
the action $S[\gamma]$ stationary with respect to variations that keep the endpoints
fixed $\delta S[\gamma] = 0$. Geometrically, this variational principle reflects
paths along which the integral of the Poincaré one-form is stationary. 
The resulting Euler-Lagrange equations are

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}^i}\right) - \frac{\partial L}{\partial q^i} = 0 .$$

The geometric interpretation of this approach includes

 1. The Poincaré one-form $\theta$ provides a natural pairing between velocities and momenta.

 2. The exterior derivative of $\theta$ gives the symplectic two-form $\omega = d\theta$, which is invariant under Hamiltonian flow.

 3. The action principle selects those paths for which the integral of $\theta - H dt$ is stationary.

The following simulation demonstrated the *action principle*.

```{r action-principle, echo=TRUE, warning=FALSE, message=FALSE, error=FALSE}
# Function to calculate the Lagrangian for a simple pendulum
pendulum_lagrangian <- function(theta, theta_dot, length=1, mass=1, g=9.8) {
  T = 0.5 * mass * (length * theta_dot)^2  # Kinetic energy
  V = mass * g * length * (1 - cos(theta))  # Potential energy
  return(T - V)  # Lagrangian L = T - V
}

# Function to calculate action for a given path
calculate_action <- function(path, dt, length=1, mass=1, g=9.8) {
  # Extract theta and compute theta_dot
  theta <- path
  theta_dot <- c(diff(theta)/dt, diff(theta)[length(theta)-1]/dt)
  
  # Calculate Lagrangian at each point
  L <- pendulum_lagrangian(theta, theta_dot, length, mass, g)
  
  # Action is integral of L over time
  S <- sum(L) * dt
  return(S)
}

# Generate a classical path and perturbed paths
generate_pendulum_paths <- function(t_max, n_steps, n_paths=10, perturbation_max=0.5) {
  t <- seq(0, t_max, length.out=n_steps)
  dt <- t[2] - t[1]
  
  # Initial conditions
  theta_0 <- pi/6  # Initial angle
  omega <- sqrt(9.8)  # Natural frequency for length=1
  
  # Classical solution for small oscillations
  theta_classical <- theta_0 * cos(omega * t)
  
  # Store all paths
  paths <- list()
  paths[[1]] <- theta_classical
  
  # Generate perturbed paths with fixed endpoints
  for(i in 2:n_paths) {
    # Create random perturbation (zero at endpoints)
    perturbation_amplitude <- perturbation_max * (i-1)/(n_paths-1)
    perturbation <- perturbation_amplitude * sin(pi * t/t_max) * sin(2*pi * t/t_max)
    
    # Perturbed path
    paths[[i]] <- theta_classical + perturbation
  }
  
  return(list(
    t = t,
    dt = dt,
    paths = paths
  ))
}

# Generate paths for pendulum
pendulum_data <- generate_pendulum_paths(
  t_max = 4,
  n_steps = 100,
  n_paths = 20
)

# Calculate action for each path
actions <- numeric(length(pendulum_data$paths))
for(i in seq_along(pendulum_data$paths)) {
  actions[i] <- calculate_action(pendulum_data$paths[[i]], pendulum_data$dt)
}

# Create data frame for plotting
plot_data <- data.frame()
path_names <- c("Classical", paste("Perturbed", 1:(length(pendulum_data$paths)-1)))

for(i in seq_along(pendulum_data$paths)) {
  path_df <- data.frame(
    time = pendulum_data$t,
    theta = pendulum_data$paths[[i]],
    path = path_names[i],
    action = actions[i]
  )
  plot_data <- rbind(plot_data, path_df)
}

# Normalize action values for color mapping
action_range <- range(actions[-1])  # Exclude classical path
normalized_actions <- (actions - action_range[1]) / diff(action_range)

# Plot the paths and their actions
trajectory_plot <- plot_ly(height=450)

# Add classical path first
classical_data <- subset(plot_data, path == "Classical")
trajectory_plot <- trajectory_plot %>%
  add_trace(
    data = classical_data,
    x = ~time,
    y = ~theta,
    type = 'scatter',
    mode = 'lines',
    name = 'Classical Path',
    line = list(color = 'black', width = 3)
  )

# Add perturbed paths
for(i in 2:length(pendulum_data$paths)) {
  path_subset <- subset(plot_data, path == path_names[i])
  # Color gradient based on action (blue=low, red=high)
  norm_action <- normalized_actions[i]
  color_val <- rgb(norm_action, 0, 1-norm_action)
  
  trajectory_plot <- trajectory_plot %>%
    add_trace(
      data = path_subset,
      x = ~time,
      y = ~theta,
      type = 'scatter',
      mode = 'lines',
      name = sprintf("%s (S = %.2f)", path_names[i], actions[i]),
      line = list(color = color_val, width = 1.5)
    )
}

# Plot actions as bars
action_df <- data.frame(
  path = factor(path_names, levels = path_names),
  action = actions
)

action_plot <- plot_ly(
  data = action_df,
  x = ~path,
  y = ~action,
  type = 'bar',
  marker = list(
    color = c('black', colorRamp(c('blue', 'red'))(normalized_actions[-1])),
    line = list(color = 'black', width = 1)
  )
) %>%
  layout(
    title = "Action Values for Different Paths",
    xaxis = list(title = "Path"),
    yaxis = list(title = "Action Value")
  )

# Combine plots
subplot(
  trajectory_plot %>% layout(title = "Pendulum Paths (Classical vs Perturbed)"),
  action_plot,
  nrows = 2,
  heights = c(0.7, 0.3)
) %>%
  layout(
    title = "Principle of Stationary Action for a Pendulum",
    showlegend = FALSE
  )
```

In this principle of stationary action simulation, we use a pendulum experiment to
show 

 1. The *black lines* representing the classical paths corresponding to the action stationary.

 2. The *colored lines* represent perturbed paths with the same endpoints, with colors indicating the action value (blue for lower, red for higher).

 3. The *bar chart* shows clearly that the classical path minimizes the action compared to the perturbed paths.

 4. The principle $\delta S = 0$ is illustrated as small perturbations away from the classical path increase the action value.

## Hamiltonian Approach and the Poincaré Two-Form

The Hamiltonian formulation builds upon the Poincaré one-form by taking its exterior derivative to obtain the symplectic two-form, which provides the foundation for understanding phase space, canonical transformations, and Noether's theorem.

The Poincaré two-form $\omega$ is obtained by taking the *exterior derivative*
of the Poincaré one-form $\theta$, i.e., $\omega = d\theta = dp_i \wedge dq^i$.
This symplectic two-form has several properties

 1. It is *non-degenerate*, allowing us to define a unique Hamiltonian vector field $X_H$ corresponding to any function $H$ on phase space.

 2. It is *closed* ($d\omega = 0$), which is a necessary condition for Liouville's theorem on the preservation of phase space volume.

 3. It defines a *Poisson bracket structure* on phase space:

$$\{F, G\} = \omega(X_F, X_G) = \frac{\partial F}{\partial q^i}\frac{\partial G}{\partial p_i} - \frac{\partial F}{\partial p_i}\frac{\partial G}{\partial q^i}.$$

The Hamiltonian vector field $X_H$ defines a flow on phase space that preserves the symplectic structure $i_{X_H}\omega = dH$, where $i_{X_H}$ denotes the interior product with $X_H$.
A function $F$ on phase space is a conserved quantity if it is invariant along the Hamiltonian flow $\frac{dF}{dt} = \{F, H\} = 0$.
Noether's theorem establishes a deep connection between symmetries and conserved quantities:

 1. *Strong Version*: If the Hamiltonian is invariant under a continuous symmetry transformation, there exists a corresponding conserved quantity.

 2. *Weak (Local) Version*: Even when the symmetry is only local or the Hamiltonian is not strictly invariant, we can still identify quantities that are conserved to first order or within certain regions of phase space.

These "weak Noether invariants" are crucial in systems with approximate symmetries or gauge theories, where they lead to constraints and generate gauge transformations.

This demonstration shows symplectic structures and conservation laws.

```{r symplectic, echo=TRUE, warning=FALSE, message=FALSE, error=FALSE}
# Function to simulate Hamiltonian dynamics
simulate_hamiltonian_system <- function(
  H_func,       # Hamiltonian function
  grad_H_func,  # Gradient of Hamiltonian
  initial_state,# Initial conditions [q1, q2, p1, p2]
  t_max,        # Maximum time
  dt,           # Time step
  invariant_func = NULL  # Optional invariant function
) {
  # Number of steps
  n_steps <- ceiling(t_max / dt)
  
  # Allocate arrays for results
  states <- matrix(0, nrow=n_steps, ncol=length(initial_state))
  times <- seq(0, t_max, length.out=n_steps)
  energy <- numeric(n_steps)
  invariant <- NULL
  if(!is.null(invariant_func)) {
    invariant <- numeric(n_steps)
  }
  
  # Initialize
  states[1,] <- initial_state
  energy[1] <- H_func(initial_state)
  if(!is.null(invariant_func)) {
    invariant[1] <- invariant_func(initial_state)
  }
  
  # Symplectic integration (4th order Runge-Kutta)
  for(i in 2:n_steps) {
    # Current state
    state <- states[i-1,]
    
    # RK4 integration
    k1 <- dt * hamiltonian_flow(state, grad_H_func)
    k2 <- dt * hamiltonian_flow(state + 0.5*k1, grad_H_func)
    k3 <- dt * hamiltonian_flow(state + 0.5*k2, grad_H_func)
    k4 <- dt * hamiltonian_flow(state + k3, grad_H_func)
    
    # Update state
    states[i,] <- state + (k1 + 2*k2 + 2*k3 + k4)/6
    
    # Calculate energy and invariant
    energy[i] <- H_func(states[i,])
    if(!is.null(invariant_func)) {
      invariant[i] <- invariant_func(states[i,])
    }
  }
  
  # Return results
  result <- list(
    times = times,
    states = states,
    energy = energy
  )
  
  if(!is.null(invariant_func)) {
    result$invariant <- invariant
  }
  
  return(result)
}

# Hamiltonian flow vector field
hamiltonian_flow <- function(state, grad_H_func) {
  # Extract position and momentum
  q <- state[1:2]
  p <- state[3:4]
  
  # Calculate gradient of Hamiltonian
  grad_H <- grad_H_func(state)
  
  # Hamilton's equations: ẋ = J∇H where J is the symplectic matrix
  flow <- c(
    grad_H[3:4],   # dq/dt = ∂H/∂p
    -grad_H[1:2]   # dp/dt = -∂H/∂q
  )
  
  return(flow)
}

# Hamiltonian for a 2D system with central force
central_force_hamiltonian <- function(state) {
  # Extract position and momentum
  q1 <- state[1]
  q2 <- state[2]
  p1 <- state[3]
  p2 <- state[4]
  
  # Kinetic energy T = p²/2m (m=1)
  T <- 0.5 * (p1^2 + p2^2)
  
  # Potential energy for central force V = k/r
  r <- sqrt(q1^2 + q2^2)
  V <- 1/r
  
  # Hamiltonian H = T + V
  return(T + V)
}

# Gradient of Hamiltonian
central_force_gradient <- function(state) {
  # Extract position and momentum
  q1 <- state[1]
  q2 <- state[2]
  p1 <- state[3]
  p2 <- state[4]
  
  # Calculate r and its cube
  r <- sqrt(q1^2 + q2^2)
  r3 <- r^3
  
  # Gradient components
  dH_dq1 <- -q1 / r3
  dH_dq2 <- -q2 / r3
  dH_dp1 <- p1
  dH_dp2 <- p2
  
  return(c(dH_dq1, dH_dq2, dH_dp1, dH_dp2))
}

# Angular momentum (Noether invariant for rotational symmetry)
angular_momentum <- function(state) {
  # L = q1*p2 - q2*p1
  return(state[1]*state[4] - state[2]*state[3])
}

# Simulate the system for different initial conditions
simulate_central_force <- function() {
  # Initial conditions (varying angular momentum)
  initial_states <- list(
    c(1.0, 0.0, 0.0, 1.0),    # Circular orbit
    c(1.0, 0.0, 0.0, 1.2),    # Elliptical orbit 1
    c(1.0, 0.0, 0.0, 0.8)     # Elliptical orbit 2
  )
  
  # Simulation parameters
  t_max <- 20
  dt <- 0.01
  
  # Run simulations
  results <- list()
  for(i in seq_along(initial_states)) {
    results[[i]] <- simulate_hamiltonian_system(
      H_func = central_force_hamiltonian,
      grad_H_func = central_force_gradient,
      initial_state = initial_states[[i]],
      t_max = t_max,
      dt = dt,
      invariant_func = angular_momentum
    )
  }
  
  return(results)
}

# Run simulations
central_force_results <- simulate_central_force()

# Prepare data for plotting
trajectory_data <- data.frame()
conservation_data <- data.frame()

orbit_names <- c("Circular", "Elliptical 1", "Elliptical 2")
colors <- c('#1f77b4', '#ff7f0e', '#2ca02c')

for(i in seq_along(central_force_results)) {
  result <- central_force_results[[i]]
  
  # Trajectory data
  traj <- data.frame(
    orbit = orbit_names[i],
    time = result$times,
    q1 = result$states[,1],
    q2 = result$states[,2],
    p1 = result$states[,3],
    p2 = result$states[,4],
    energy = result$energy,
    angular_momentum = result$invariant
  )
  
  trajectory_data <- rbind(trajectory_data, traj)
  
  # Conservation data (sampled)
  sample_idx <- seq(1, length(result$times), by=10)
  conservation <- data.frame(
    orbit = orbit_names[i],
    time = result$times[sample_idx],
    energy = result$energy[sample_idx],
    angular_momentum = result$invariant[sample_idx],
    normalized_energy = (result$energy[sample_idx] - min(result$energy))/
                        (max(result$energy) - min(result$energy)),
    normalized_momentum = (result$invariant[sample_idx] - min(result$invariant))/
                          (max(result$invariant) - min(result$invariant))
  )
  
  conservation_data <- rbind(conservation_data, conservation)
}

# Create plots
# 1. Phase space trajectories
phase_space_plot <- plot_ly(height=400)

for(i in 1:length(orbit_names)) {
  orbit_data <- subset(trajectory_data, orbit == orbit_names[i])
  
  phase_space_plot <- phase_space_plot %>%
    add_trace(
      data = orbit_data,
      x = ~q1,
      y = ~q2,
      type = 'scatter',
      mode = 'lines',
      line = list(color = colors[i], width = 2),
      name = orbit_names[i]
    )
}

phase_space_plot <- phase_space_plot %>%
  layout(
    title = "Configuration Space Trajectories",
    xaxis = list(title = "q₁"),
    yaxis = list(title = "q₂")
  )

# 2. Conservation plots
energy_plot <- plot_ly(height=300)
momentum_plot <- plot_ly(height=300)

for(i in 1:length(orbit_names)) {
  orbit_data <- subset(conservation_data, orbit == orbit_names[i])
  
  energy_plot <- energy_plot %>%
    add_trace(
      data = orbit_data,
      x = ~time,
      y = ~normalized_energy,
      type = 'scatter',
      mode = 'lines',
      line = list(color = colors[i], width = 2),
      name = orbit_names[i]
    )
  
  momentum_plot <- momentum_plot %>%
    add_trace(
      data = orbit_data,
      x = ~time,
      y = ~normalized_momentum,
      type = 'scatter',
      mode = 'lines',
      line = list(color = colors[i], width = 2),
      name = orbit_names[i]
    )
}

energy_plot <- energy_plot %>%
  layout(
    title = "Energy Conservation (Normalized)",
    xaxis = list(title = "Time"),
    yaxis = list(title = "Normalized Energy")
  )

momentum_plot <- momentum_plot %>%
  layout(
    title = "Angular Momentum Conservation (Normalized)",
    xaxis = list(title = "Time"),
    yaxis = list(title = "Normalized Angular Momentum")
  )

# 3. Poincaré section visualization (2D section of phase space)
poincare_data <- data.frame()

for(i in 1:length(orbit_names)) {
  orbit_data <- subset(trajectory_data, orbit == orbit_names[i])
  
  # Find points where trajectory crosses q2=0 plane with p2>0
  # Modified approach to ensure we get some crossings
  crossing_indices <- c()
  for(j in 1:(nrow(orbit_data)-1)) {
    # Check for crossing from negative to positive q2
    if(orbit_data$q2[j] <= 0 && orbit_data$q2[j+1] > 0) {
      crossing_indices <- c(crossing_indices, j)
    }
    # Also check for crossing from positive to negative q2
    # to ensure we capture some points for each orbit
    else if(orbit_data$q2[j] >= 0 && orbit_data$q2[j+1] < 0) {
      crossing_indices <- c(crossing_indices, j)
    }
  }
  
  # If we still don't have crossings, sample some points
  if(length(crossing_indices) < 3) {
    # Take points closest to q2=0
    dist_to_plane <- abs(orbit_data$q2)
    ordered_indices <- order(dist_to_plane)
    crossing_indices <- ordered_indices[1:10]  # Take 10 closest points
  }
  
  # Create Poincaré section data
  if(length(crossing_indices) > 0) {
    crossings <- orbit_data[crossing_indices,]
    crossings$section <- "Near q₂ = 0"
    poincare_data <- rbind(poincare_data, crossings)
  }
}

poincare_plot <- plot_ly(height=400)

for(i in 1:length(orbit_names)) {
  orbit_data <- subset(poincare_data, orbit == orbit_names[i])
  
  if(nrow(orbit_data) > 0) {
    poincare_plot <- poincare_plot %>%
      add_trace(
        data = orbit_data,
        x = ~q1,
        y = ~p1,
        type = 'scatter',
        mode = 'markers',
        marker = list(color = colors[i], size = 10),
        name = orbit_names[i]
      )
  }
}

poincare_plot <- poincare_plot %>%
  layout(
    title = "Poincaré Section (q₂ = 0, p₂ > 0)",
    xaxis = list(title = "q₁"),
    yaxis = list(title = "p₁")
  )

# Combine all plots
subplot(
  phase_space_plot,
  subplot(energy_plot, momentum_plot, nrows=1),
  poincare_plot,
  nrows = 3,
  heights = c(0.4, 0.3, 0.3)
) %>%
  layout(
    title = "Symplectic Structure and Noether Invariants",
    showlegend = TRUE
  )
```

The Hamiltonian simulation with Noether invariants shows

 1. *Phase Space Structure*: The *top plot* shows trajectories in configuration space for a central force system (like a gravitational or Coulomb potential). The circular and elliptical orbits arise from the same Hamiltonian with different initial conditions.

 2. *Conservation Laws*: The *middle plots* demonstrate the conservation of energy and angular momentum. Angular momentum is the Noether invariant corresponding to rotational symmetry of the central force potential.

 3. *Symplectic Structure*: The *bottom plot* shows a Poincaré section, which provides a means to visualize the phase space structure by recording points where trajectories intersect a chosen plane (in this case, when $q_2= 0$ with $p_2 > 0$). The invariant curves in the Poincaré section reflect the preservation of the symplectic form and phase space volume (Liouville's theorem).

The preservation of these invariants, even in this discrete numerical simulation,
demonstrates how the symplectic structure fundamentally constrains the dynamics
in phase space.


## Complex Time (Kime) and Analytic Continuation

Complex time, or kime ($\kappa$), is defined as $\kappa = t e^{i\theta}$, where 
$t$ is the standard time parameter and $\theta$ is a random kime-phase distributed
according to a phase distribution $\Phi(\theta)$ supported on $[-\pi, \pi)$. 
This formulation extends time into the complex plane, providing a mechanism for
connecting quantum and classical mechanics through regularizing singular potentials.

Complex time appears in several contexts in theoretical physics

 1. *Kime-Phase Distribution*: The phase $\theta$ is sampled from a distribution $\Phi(\theta)$ supported on $[-\pi, \pi)$, characterizing the stochastic nature of the complex time extension.

 2. *Wick Rotation*: Traditionally viewed as the transformation $t \to i\tau$ (rotating time into the imaginary axis), this can be understood through kime as a special case where $\theta = \pi/2$. This transforms between Minkowski and Euclidean space-time, and the action transforms as $S_M \to iS_E$, where $S_M$ is the Minkowski action and $S_E$ is the Euclidean action.

 3. *Path Integrals*: In the path integral formulation, the quantum amplitude is $\langle q_f, t_f | q_i, t_i \rangle = \int \mathcal{D}q \, e^{iS[q]/\hbar}$. When examined through the kime formulation with appropriate phase distribution $\Phi(\theta)$, this allows for a more general analysis beyond the standard Wick rotation.

 4. *Analytic Continuation of Classical Trajectories*: Classical trajectories can be extended into the complex plane via kime, allowing for the description of classically forbidden processes like tunneling, with the phase distribution $\Phi(\theta)$ determining the probability of various tunneling pathways.

The Hamilton-Jacobi equation can be analytically continued into complex time

$$\frac{\partial S}{\partial \kappa} + H\left(q, \frac{\partial S}{\partial q}, \kappa\right) = 0,$$

where $\kappa = t e^{i\theta}$ is the complex time variable with $\theta \sim \Phi(\theta)$.

Complex-time trajectories with stochastic phase distributions provide solutions that connect classical turning points, enabling

 1. *Stochastic WKB Approximation*: Improved accuracy in semiclassical approximations by incorporating phase distribution $\Phi(\theta)$
 2. *Distributed Instantons*: Nonperturbative tunneling solutions with phase-dependent pathways
 3. *Extended Resurgence Theory*: Connection between perturbative and nonperturbative effects through phase-distribution-weighted contributions

Here is a quantum tunneling simulation using complex time.

```{r complex-time, echo=TRUE, warning=FALSE, message=FALSE, error=FALSE}
# Function to generate potential
double_well_potential <- function(x, a=1, b=0.25) {
  return(a*x^4 - b*x^2)
}

# Calculate complex trajectories for tunneling
calculate_complex_trajectories <- function(E=0.1, a=1, b=0.25, resolution=50) {
  # Find classical turning points by solving V(x) = E
  # For double-well potential, we need numerical solution
  find_turning_points <- function(E) {
    # Starting points for numerical search
    x_guess <- c(-1, -0.1, 0.1, 1)
    turning_points <- numeric(0)
    
    for(x_start in x_guess) {
      # Find root where V(x) = E
      result <- uniroot(function(x) double_well_potential(x, a, b) - E, 
                         c(x_start - 0.5, x_start + 0.5), 
                         tol=1e-10, extendInt="yes")
      turning_points <- c(turning_points, result$root)
    }
    
    # Remove duplicates and sort
    turning_points <- unique(turning_points)
    turning_points <- turning_points[order(turning_points)]
    
    return(turning_points)
  }
  
  turning_points <- find_turning_points(E)
  
  # Define integration paths in complex time plane
  # We'll use 3 segments:
  # 1. Real time from left turning point inwards
  # 2. Imaginary time from left to right (tunneling)
  # 3. Real time from right turning point outwards
  
  # Setup grid
  x_range <- range(turning_points) * 1.5
  x_grid <- seq(x_range[1], x_range[2], length.out=resolution)
  
  # Calculate potential on grid
  V_grid <- double_well_potential(x_grid, a, b)
  
  # Path 1: Real time trajectory (left well)
  x_left <- seq(turning_points[1], -0.1, length.out=resolution/3)
  t_left <- cumsum(c(0, sqrt(2/a) * diff(x_left) / 
                    sqrt(abs(E - double_well_potential(x_left[-length(x_left)], a, b)))))
  
  # Path 2: Imaginary time trajectory (tunneling)
  x_tunnel <- seq(-0.1, 0.1, length.out=resolution/3)
  tau_tunnel <- cumsum(c(0, sqrt(2/a) * diff(x_tunnel) / 
                      sqrt(abs(double_well_potential(x_tunnel[-length(x_tunnel)], a, b) - E))))
  
  # Path 3: Real time trajectory (right well)
  x_right <- seq(0.1, turning_points[4], length.out=resolution/3)
  t_right <- cumsum(c(0, sqrt(2/a) * diff(x_right) / 
                     sqrt(abs(E - double_well_potential(x_right[-length(x_right)], a, b)))))
  
  # Combine paths in complex time plane
  kappa_real <- c(t_left, rep(max(t_left), length(tau_tunnel)), max(t_left) + t_right)
  kappa_imag <- c(rep(0, length(t_left)), tau_tunnel, rep(max(tau_tunnel), length(t_right)))
  x_combined <- c(x_left, x_tunnel, x_right)
  
  # Calculate tunneling amplitude
  S_tunnel <- sum(sqrt(2 * a) * diff(x_tunnel) * 
                 sqrt(abs(double_well_potential(x_tunnel[-length(x_tunnel)], a, b) - E)))
  tunneling_probability <- exp(-2 * S_tunnel / sqrt(a))
  
  return(list(
    turning_points = turning_points,
    x_grid = x_grid,
    V_grid = V_grid,
    x_combined = x_combined,
    kappa_real = kappa_real,
    kappa_imag = kappa_imag,
    tunneling_probability = tunneling_probability,
    energy = E
  ))
}

# Generate tunneling data
tunnel_data <- calculate_complex_trajectories(E=0.05)

# Create 3D Visualization of Complex Time Tunneling
# 1. Potential Plot
potential_plot <- plot_ly(height=300) %>%
  add_trace(
    x = tunnel_data$x_grid,
    y = tunnel_data$V_grid,
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'black', width = 3),
    name = 'Double-Well Potential'
  ) %>%
  add_trace(
    x = c(min(tunnel_data$x_grid), max(tunnel_data$x_grid)),
    y = c(tunnel_data$energy, tunnel_data$energy),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'red', dash = 'dash', width = 2),
    name = 'Energy Level'
  ) %>%
  add_trace(
    x = tunnel_data$turning_points,
    y = rep(tunnel_data$energy, length(tunnel_data$turning_points)),
    type = 'scatter',
    mode = 'markers',
    marker = list(color = 'blue', size = 10, symbol = 'circle'),
    name = 'Turning Points'
  ) %>%
  layout(
    title = paste0("Double-Well Potential with E = ", 
                  tunnel_data$energy,
                  " (Tunneling Probability ≈ ", 
                  format(tunnel_data$tunneling_probability, digits=4),
                  ")"),
    xaxis = list(title = "Position x"),
    yaxis = list(title = "Potential V(x)")
  )

# 2. Complex Time Trajectory
complex_time_plot <- plot_ly(height=350) %>%
  add_trace(
    x = tunnel_data$kappa_real,
    y = tunnel_data$kappa_imag,
    z = tunnel_data$x_combined,
    type = 'scatter3d',
    mode = 'lines',
    line = list(color = 'blue', width = 5),
    name = 'Complex Time Trajectory'
  ) %>%
  layout(
    scene = list(
      xaxis = list(title = "Re(κ)"),
      yaxis = list(title = "Im(κ)"),
      zaxis = list(title = "Position x")
    ),
    title = "Quantum Tunneling in Complex Time"
  )

# 3. Projection plots
real_time_plot <- plot_ly(height=250) %>%
  add_trace(
    x = tunnel_data$kappa_real,
    y = tunnel_data$x_combined,
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'green', width = 3),
    name = 'Real Time Projection'
  ) %>%
  layout(
    title = "Trajectory Projection on Real Time Plane",
    xaxis = list(title = "Re(κ)"),
    yaxis = list(title = "Position x")
  )

imag_time_plot <- plot_ly(height=250) %>%
  add_trace(
    x = tunnel_data$kappa_imag,
    y = tunnel_data$x_combined,
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'purple', width = 3),
    name = 'Imaginary Time Projection'
  ) %>%
  layout(
    title = "Trajectory Projection on Imaginary Time Plane",
    xaxis = list(title = "Im(κ)"),
    yaxis = list(title = "Position x")
  )

# Combine all plots
subplot(
  potential_plot,
  complex_time_plot,
  subplot(real_time_plot, imag_time_plot, nrows=1),
  nrows = 3,
  heights = c(0.3, 0.4, 0.3)
) %>%
  layout(
    title = "Complex Time (Kime) Analysis of Quantum Tunneling",
    showlegend = FALSE
  )
```

This simulation illustrates quantum tunneling through a double-well potential using complex time analysis. In the context of the kime representation ($\kappa = t e^{i\theta}$), we can interpret this as follows

 1. *Potential Landscape*: The *top plot* shows the double-well potential with an energy level (red dashed line) below the barrier height. Classical motion is confined to either the left or right well, bounded by turning points (blue dots). This represents the system configuration in purely real time.

 2. *Complex Time Trajectory*: The *middle plot* shows a complex-time trajectory connecting the two classically allowed regions. In the kime framework, this represents a path where the phase $\theta$ varies from 0 (real time) to values approaching $\pi/2$ (imaginary time) for tunneling, and back to 0. This trajectory is a particular realization from the phase distribution $\Phi(\theta)$.

3. *Projections*: The *bottom plots* show projections of the trajectory onto the real and imaginary time planes, highlighting how the kime-phase $\theta$ modulates between classical motion ($\theta \approx 0$) and quantum tunneling ($\theta \approx \pi/2$).

The tunneling probability, calculated from the imaginary action along the complex path, represents an integral over the possible phase paths weighted by the distribution $\Phi(\theta)$. In more general treatments, the phase distribution $\Phi(\theta)$ would determine an ensemble of possible tunneling pathways, each with its own probability contribution.

## Connections and Unifying Framework

These four mathematical frameworks—Hamilton-Jacobi theory, Lagrangian formalism with the Poincaré one-form, Hamiltonian mechanics with the symplectic two-form, and complex time analysis—provide a cohesive picture of classical and quantum mechanics. Here's how they interconnect

The Hamilton-Jacobi equation
$\frac{\partial S}{\partial t} + H\left(q, \frac{\partial S}{\partial q}, t\right) = 0$
becomes the quantum Hamilton-Jacobi equation when we write the wavefunction as $\Psi = e^{iS/\hbar}$ and substitute into the Schrödinger equation

$$\frac{\partial S}{\partial t} + H\left(q, \frac{\partial S}{\partial q}, t\right) = \frac{\hbar^2}{2m} \frac{1}{\sqrt{\rho}} \nabla^2 \sqrt{\rho},$$

where $\rho = |\Psi|^2$ is the probability density. The right-hand term, which vanishes as $\hbar \to 0$, is known as the "quantum potential" and accounts for quantum effects.

The Poincaré one-form $\theta = p_i dq^i$ generates the symplectic two-form $\omega = d\theta = dp_i \wedge dq^i$ through exterior differentiation. The action principle involving the one-form, 
$\delta \int \left( \theta - H \, dt \right) = 0$, leads to Hamilton's equations, 
which preserve the symplectic two-form

$$\frac{dq^i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial q^i}.$$

Eextending time into the complex plane, the action becomes a complex-valued function. The paths of stationary action in complex time can describe both classical dynamics and quantum tunneling. The transformation between Minkowski and Euclidean formulations (Wick rotation) connects quantum mechanics to statistical mechanics through $Z = \text{Tr}(e^{-\beta H}) = \int \mathcal{D}q \, e^{-S_E[q]/\hbar}$, 
where $Z$ is the partition function and $\beta = 1/k_BT$ is the inverse temperature.

Below is a unified phase space simulation with quantum corrections.

```{r unified, echo=TRUE, warning=FALSE, message=FALSE, error=FALSE}
# Generate unified phase space visualization with quantum corrections
generate_unified_phase_space <- function(hbar_values = c(0, 0.1, 0.5), resolution = 50) {
  # Define phase space grid
  q_range <- c(-2, 2)
  p_range <- c(-2, 2)
  
  q_grid <- seq(q_range[1], q_range[2], length.out = resolution)
  p_grid <- seq(p_range[1], p_range[2], length.out = resolution)
  
  grid_points <- expand.grid(q = q_grid, p = p_grid)
  
  # Classical Hamiltonian (harmonic oscillator)
  H_classical <- function(q, p, m = 1, k = 1) {
    return(p^2/(2*m) + 0.5*k*q^2)
  }
  
  # Quantum corrected Hamiltonian with "quantum potential" term
  H_quantum <- function(q, p, hbar, m = 1, k = 1) {
    # Classical term
    H_cl <- p^2/(2*m) + 0.5*k*q^2
    
    # Quantum correction (simplified form of quantum potential)
    # In full quantum mechanics, this would be more complex
    H_q <- hbar^2/(8*m) * (k/m) / (p^2/(2*m) + 0.5*k*q^2 + 0.001)
    
    return(H_cl + H_q)
  }
  
  # Calculate Hamiltonians for each hbar value
  results <- list()
  
  for(i in seq_along(hbar_values)) {
    hbar <- hbar_values[i]
    
    # Calculate energy at each grid point
    if(hbar == 0) {
      # Classical case
      energy <- apply(grid_points, 1, function(point) {
        H_classical(point[1], point[2])
      })
    } else {
      # Quantum case
      energy <- apply(grid_points, 1, function(point) {
        H_quantum(point[1], point[2], hbar)
      })
    }
    
    # Vector field for phase space flow
    flow_q <- grid_points$p  # dq/dt = ∂H/∂p ≈ p
    
    # Classical flow for p
    flow_p_classical <- -grid_points$q  # dp/dt = -∂H/∂q ≈ -q
    
    # Quantum correction to flow (simplified)
    if(hbar == 0) {
      flow_p <- flow_p_classical
    } else {
      # Add quantum correction
      quantum_correction <- -hbar^2/(4) * grid_points$q /
                           (grid_points$p^2 + grid_points$q^2 + 0.001)^2
      flow_p <- flow_p_classical + quantum_correction
    }
    
    # Normalize flow for visualization
    flow_magnitude <- sqrt(flow_q^2 + flow_p^2)
    flow_q_norm <- flow_q / (flow_magnitude + 0.001)
    flow_p_norm <- flow_p / (flow_magnitude + 0.001)
    
    # Store results
    results[[i]] <- list(
      hbar = hbar,
      q = grid_points$q,
      p = grid_points$p,
      energy = energy,
      flow_q = flow_q_norm,
      flow_p = flow_p_norm
    )
  }
  
  return(list(
    q_grid = q_grid,
    p_grid = p_grid,
    results = results
  ))
}

# Generate data
phase_space_data <- generate_unified_phase_space(
  hbar_values = c(0, 0.2, 0.5),
  resolution = 30
)

# Create visualization
phase_space_plots <- list()
hbar_labels <- c("Classical (ħ = 0)", "Semiclassical (ħ = 0.2)", "Quantum (ħ = 0.5)")
plot_colors <- c("#1f77b4", "#ff7f0e", "#2ca02c")

for(i in seq_along(phase_space_data$results)) {
  result <- phase_space_data$results[[i]]
  
  # Reshape energy for contour plot
  energy_matrix <- matrix(result$energy, 
                         nrow = length(phase_space_data$q_grid),
                         ncol = length(phase_space_data$p_grid))
  
  # Create contour plot
  plot_data <- expand.grid(
    q = phase_space_data$q_grid,
    p = phase_space_data$p_grid
  )
  plot_data$energy <- c(energy_matrix)
  
  # Subsample flow vectors for clearer visualization
  subsample <- seq(1, length(result$q), by = 5)
  
  # Create plot
  plot <- plot_ly(height = 400) %>%
    add_contour(
      z = energy_matrix,
      x = phase_space_data$q_grid,
      y = phase_space_data$p_grid,
      contours = list(
        start = 0,
        end = 3,
        size = 0.25
      ),
      colorscale = 'Viridis',
      showscale = FALSE
    ) %>%
    add_trace(
      x = result$q[subsample],
      y = result$p[subsample],
      u = result$flow_q[subsample],
      v = result$flow_p[subsample],
      type = 'scatter',
      mode = 'markers+lines',
      marker = list(size = 2, color = plot_colors[i]),
      line = list(width = 1, color = plot_colors[i]),
      name = 'Flow',
      showlegend = FALSE
    ) %>%
    layout(
      title = hbar_labels[i],
      xaxis = list(title = "Position (q)"),
      yaxis = list(title = "Momentum (p)")
    )
  
  phase_space_plots[[i]] <- plot
}

# Integration illustration
integration_plot <- plot_ly(height = 300) %>%
  add_trace(
    x = seq(-2, 2, length.out = 100),
    y = cos(pi * seq(-2, 2, length.out = 100)),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'blue', width = 3),
    name = 'Classical Path'
  ) %>%
  add_trace(
    x = seq(-2, 2, length.out = 100),
    y = cos(pi * seq(-2, 2, length.out = 100)) + 
       0.3 * sin(3 * pi * seq(-2, 2, length.out = 100)),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'red', width = 3),
    name = 'Quantum Path'
  ) %>%
  layout(
    title = "Path Integration: Classical vs Quantum",
    xaxis = list(title = "Time"),
    yaxis = list(title = "Position")
  )

# Combine all plots
subplot(
  subplot(phase_space_plots[[1]], phase_space_plots[[2]], phase_space_plots[[3]], 
          nrows = 1, shareY = TRUE, shareX = TRUE),
  integration_plot,
  nrows = 2,
  heights = c(0.7, 0.3)
) %>%
  layout(
    title = "Unified Framework: From Classical to Quantum Mechanics",
    showlegend = TRUE
  )
```

the graph illustrates the unified framework connecting classical and quantum mechanics

 1. *Phase Space Evolution*: The *top row* shows phase space representations of a harmonic oscillator with increasing values of $\hbar$: (1) classical ($\hbar = 0$), perfect elliptical contours representing energy conservation; (2) semiclassical ($\hbar = 0.2$), slightly modified contours due to quantum corrections; and (3) quantum ($\hbar = 0.5$), more substantial deformation of phase space structure.

 2. *Path Integration*: The *bottom plot *compares classical and quantum paths. While classical mechanics selects a single stationary path, quantum mechanics involves integration over neighboring paths, allowing phenomena like tunneling and interference.

The progressive deformation of phase space as $\hbar$ increases demonstrates how 
quantum effects modify, but don't completely abolish, the underlying symplectic
structure of classical mechanics.

To summarize

 1. *Hamilton-Jacobi Theory* gives us a powerful tool for solving mechanical problems by finding complete integrals, connecting classical mechanics to the WKB approximation in quantum mechanics.

 2. *The Lagrangian Approach with the Poincaré One-Form* provides a geometric interpretation of the principle of stationary action, emphasizing the importance of the path itself.

 3. *The Hamiltonian Approach with the Symplectic Two-Form* shifts focus to phase space and its invariant structure, revealing how symmetries lead to conservation laws through Noether's theorem.

 4. *Complex Time Analysis* bridges quantum and classical regimes by extending dynamics into the complex plane, explaining tunneling phenomena and connecting quantum mechanics to statistical mechanics.


# References

 1. Goldstein, H., Poole, C., & Safko, J. (2002). *Classical Mechanics* (3rd ed.). Addison Wesley.

 2. Peskin, M. E., & Schroeder, D. V. (1995). *An Introduction to Quantum Field Theory*. Westview Press.

 3. Nakahara, M. (2003). *Geometry, Topology and Physics* (2nd ed.). Institute of Physics Publishing.

 4. Weinberg, S. (1995). *The Quantum Theory of Fields, Volume 1: Foundations*. Cambridge University Press.

 5. Zinn-Justin, J. (2021). *Path Integrals in Quantum Mechanics*. Oxford University Press.

 6. Arnold, V. I. (1989). *Mathematical Methods of Classical Mechanics* (2nd ed.). Springer-Verlag.
 
 7. Dinov, ID and Velev, MV (2021) *Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics*, De Gruyter (STEM Series), Berlin/Boston, ISBN 9783110697803 / 3110697807, DOI: 10.1515/9783110697827. 
 
 8. [Time Complexity and Inferencial Uncertainty (TCIU)](https://www.socr.umich.edu/TCIU/).
 
 <!--html_preserve-->
<div>
    	<footer><center>
			<a href="https://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="https://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="https://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->