为什么 sx，sy，sz（电子自旋）会影响我的代码中的 x，y，z？（虽然不应该）

问：

我试图绘制电子在尾流场中的捕获，当电子在 t=0 和 v（t=0） 时位于气泡的边缘时，我想知道电子在每个位置的自旋。此代码应返回不同时间点的位置、速度和旋转。但是位置和旋转不应该相互干扰，但它以某种方式会相互干扰。有人可以解释一下这个问题吗？

这是我的代码。在我之前的代码中，我测试了位置和速度的常微分方程的实现，而没有自旋的常微分方程，反之亦然，它有效，但现在组合代码不起作用。版本 1 显示了预期的正确图，版本 2 是有问题的图。

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import solve_ivp,odeint, quad
import matplotlib.pyplot as plt

c = 1.
m = 1.
e = 1.
q = 1.
a = 1e-3

def tbmt(t,X,Vb):#X includes x,y,z,dx,dy,dz,sx,sy,sz
    x=X[0]
    y=X[1]
    z=X[2]
    p_x=X[3]*m
    p_y=X[4]*m
    p_z=X[5]*m
    
    v=np.array([X[3],X[4],X[5]])
    L_fak = np.sqrt(1 + p_x**2 + p_y**2)
    
    dx = p_x / L_fak - Vb
    dy = p_y / L_fak
    dz = p_z / L_fak
    dpx = -(1 + Vb) * x / 4 + (p_z * z) / (4 * L_fak)
    dpy = -(1 + p_x / L_fak) * y / 4
    dpz = -(1 + p_x / L_fak) * z / 4
    
    E_x=0.5*x
    E_y=0.25*y
    B_z=-E_y
    B_x=0.
    E_z=0.25*z
    B_y=E_z
    
    E=np.array([E_x,E_y,E_z])
    B=np.array([B_x,B_y,B_z])
    
    s=np.array([1,0,0]) #Version 1 with correct plot for the electron trapping
    #s=np.array([X[6],X[7],X[8]]) #Version 2 with the wrong plot for the electron trapping
    
    gamma = 1 / np.sqrt(1 - np.dot(v/c, v/c))
    OmegaB = a + 1. / gamma
    OmegaV = a * gamma / (gamma + 1.)
    OmegaE = a + 1. / (1. + gamma)

    Omega = q * e / m / c * (OmegaB * B - OmegaV * np.dot(v/c, B) * v / c - OmegaE * np.cross(v/c, E))
    dsdt = - np.cross(Omega, s)
    
    dsx=dsdt[0]
    dsy=dsdt[1]
    dsz=dsdt[2]
    
    return np.array([dx, dy,dz, dpx, dpy, dpz]) #Version 1 with correct plot for the electron trapping
    #return np.array([dx, dy,dz, dpx, dpy, dpz,dsx,dsy,dsz]) #Version 2 with the wrong plot for the electron trapping

r_b=10 #radius of the bubble
rho=9 #z(t=0)
vx0=0 # v_x(t=0)
vy0=0 # v_y(t=0)
vz0=0 # v_z(t=0)
gamma0=4 

x0=np.sqrt(r_b**2-rho**2)
y0=np.sqrt(r_b**2-rho**2)
z0=rho

V=np.sqrt(1-(1/gamma0**2))
params=(V,)

s0=np.array([x0,y0,z0,vx0,vy0,vz0]) #Version 1 with correct plot for the electron trapping
#s0=np.array([x0,y0,z0,vx0,vy0,vz0,1.,0.,0.]) #Version 2 with the wrong plot for the electron trapping


Tend = 300  
N = 1000    

t_span = (0.0, Tend)               
tau = np.linspace(0, Tend, N+1) 

sol=solve_ivp(tbmt, t_span, s0, t_eval = tau, args = params)

x = sol.y[0]    
y = sol.y[1]    
z = sol.y[2] 
vx = sol.y[3]   
vy = sol.y[4]   
vz = sol.y[5] 
t = sol.t  

#plot
plt.figure(dpi = 200)
plt.plot(x, z, color = 'r')
plt.gca().add_patch(plt.Circle((0, 0), r_b, ec = 'black', fc = 'none', linewidth = 2))
plt.axis('equal')
plt.xlim(-r_b, r_b)
plt.ylim(-r_b, r_b)
plt.xlabel('$\zeta$')
plt.ylabel('y')
#plt.grid()
plt.show()

版本 1：我应该从代码中期待什么 版本 2：轨迹错误