Relation between initial spin(nors) and kinetic

[jesusmontera]

Hello, in your sample rel_particle2d.py the initial spinors are created based on KX ans KY

init_spinor = np.array([ones*(mc*px - 1.0j*mc*py + 
                                  (px - 1.0j*py)*omega)/(p*den),
                        zeros, zeros, ones*p/den], np.complex128)

is this a restriction to real phisycs?, or the spin can be anyone you like independent from Kinetic, like:

 sv = Statevector.from_label(statelabel) #z='0' or '1' x='+' or '-'  y='r' or  'l'
 statesqubit=sv.data 
init_spinor = np.array([ones * statesqubit[0],  zeros, zeros,
                            ones*statesqubit[1]],   np.complex128)

[marl0ny]

My reasoning for choosing that form for init_spinor is that it corresponds to a positive energy eigenstate for the free particle Hamiltonian. I think you can set the four components of the wave function to whatever you want, so long as boundary conditions are preserved, the initial wave function is a valid state in your Hilbert space, etc.

For the forms of the positive (and negative) energy eigenstates, I listed and used them here, which I computed using this script, using Sympy. They are also explained in Wikipedia, although they may have expressed it differently due to the degeneracies in the positive/negative energies, and they may have used a different choice for the gamma matrices.

[jesusmontera]

Thanks, I understood that because you initialized psi0 in splitstep with a specific Kinectic on it's spinors, then each K reflects a specific spin, and if this spin is disturbed (for example a magnetic field), then K must be disturbed too. And the operations to init psi with specific spinors by Kinetic are: psi[0] pos_eig1[1] (spin + up) psi[3] pos_eig1[2] (spin + down) But if i want to initialize psi also with negative I suppose it will be: psi[2] neg_eig1[1] ( spin - up) psi[1] neg_eig1[2] ( spin - down) That will be:

init_spinor = np.array( [pos_eig1[1], neg_eig1[2],
                        neg_eig1[1], pos_eig1[2]], np.complex128) * ones

psi0 *=  init_spinor

No? I tested with positive spinors in 3D diracsplitstep with a variant magnetic field in Z B and Kinetic Y by doing a np.dot after U step between the initial spin (in 3d bloch vector) and B. The spin changes and the particle positon looks ok like Stern-Gerlach. But if I add negative spinors, I will have 2 spins? or I have to add up + with neg + Note the spins are extracted:

up = np.sum(psi[0]) + np.sum(psi[2])
down = np.sum(psi[3]) + np.sum(psi[1])     

[marl0ny]

To add some more explanation to why I initialized psi0 the way I did, a wave packet consisting of positive energy eigenstates have a group velocity that's in the same direction for my initial choice of momentum. If I did this with negative energy eigenstates instead, the wave packet would propagate in the opposite direction. In general psi0 can consist of both a combination of positive and negative energy eigenstates, where in this case an initial wave packet may split in two with each piece propagating in opposite directions.

As mentioned in these notes, states of definite spin are in general not energy eigenstates. However, for positive energy solutions and for speeds much less than the speed of light, the first two components are significantly larger than the last two, and in this limit the first two components can be well approximated by the familiar two component spinor as described by the Pauli equation. I think this works as well for the negative energy solutions, but with the bottom two components dominating over the top two.

Perhaps you could try is something like (pos_eig1 @ np.array([spin[0], spin[1], 0.0, 0.0]))*np.conj(pos_eig1) + pos_eig2 @ np.array([spin[0], spin[1], 0.0, 0.0])*np.conj(pos_eig2), or (neg_eig1 @ np.array([0.0, 0.0, spin[0], spin[1]]))*np.conj(neg_eig1) + neg_eig2 @ np.array([0.0, 0.0, spin[0], spin[1]])*np.conj(neg_eig2), then normalize it (I should have mentioned that the the energy eigenstates are actually In bra form and should be conjugated to get them in ket form). This is basically just projecting spin states to energy eigenstates of a given momentum.

[jesusmontera]

Thanks, ok, so since u pass any spin to make the initial spinnors, that mean's that not always the spin is consequence of momentum, but in that special case plotting it in Bloch sphere points to the moving direction. But in the sample rel_particle2d.py (first message) u initialize spinnors [0] and [3], which I suppose to be positive up down, and in the last message u use [0] and [1] for positive, so my doubt is witch expect the dirac_splitstep.py for positive, [0][3] or [0][1]?, because there is multiplitcation psi0 * initspinors before starting the U steps process.

[marl0ny]

In my comment above, I could have written it also as c1 = pos_eig1 @ np.array([spin[0], spin[1], 0.0, 0.0]) # inner product c2 = pos_eig2 @ np.array([spin[0], spin[1], 0.0, 0.0]) init_4spinor = c1 np.conj(pos_eig1) + c2 np.conj(pos_eig2) where the third and fourth components are not necessarily zero.

I can take any linear combination of pos_eig1 and pos_eig2 to get a positive energy solution since they are degenerate, where in general each of the components can be non zero.

[jesusmontera]

Thanks, now giving + or - spinnors is ok, but for extracting spin from spinnors (or psi at each step) like I said does not work:

up = np.sum(psi[0]) + np.sum(psi[2])
down = np.sum(psi[3]) + np.sum(psi[1])   

so I put spin up is s[0] and spin down is s[3] when spinnor(s) is initialized positive, and up is s[2] and down s[1] when spinnor is init negative.

[marl0ny]

I don't know if this resolves your issue, but the term containing σ•B only appears in the Pauli equation once you take the non-relativistic limit of the Dirac equation for positive energy solutions, where the Pauli equation only acts on the first two components of the Dirac spinor. This doesn't apply to the last two components, which work differently (if χ is the spinor containing the first two components and φ is for the bottom two, then φ ~ σ•p χ, where p is its momentum.

I think the magnetic fields can only be implemented properly through appropriately setting the vector potential (it's the vector_potential parameter in the initializer), although you may find that those vector potential terms that do not contribute to the magnetic field tend to dominate.

[jesusmontera]

I took the idea from this simulation with B here, so I did the same for the 4 Dirac psi's and the particle moves ok. The dot is between B[x][y][z] and the initial spin and add that complex contribution after U step. I tried with B and each spìn at [i][x][y][z] , but was before fixing spinnors with bad result. Now I'm modeling the magnets with magpylib. . My goal is to study spins under entanglement and B, so I modified your two_particles2d.py to DiracSplitStep, to apply a B field and extract spins for each particle, but I still working in how to separate (trace) the 4D psi to extract th particles individual spins. I suppose they must have simetricall positions (diagonal) for the trace to work.

[marl0ny]

I skimmed through the pdf in the linked repo, and they are indeed using the Pauli equation, but without the vector potential A(r) in the momentum term. If you just want to look at spin under a magnetic field then I think just using the Pauli equation will work just fine. Approximating it with the split operator method, and without the vector potential in the momentum term, then Ψ(r, Δt) = exp(- i H Δt / ħ) Ψ(r, 0) ≈ exp(i q σ• B(r) Δt/2) exp(- i p² Δt / (2 m ħ) ) exp(i q σ• B(r) Δt/2) Ψ(r, 0). I think you can use this (equation 2) to find the form of exp(i q σ• B(r) Δt/2).

[jesusmontera]

Thanks, I will take look at Pauli equation to compare results with your Dirac equation. Now I'm rendering it in OpenGL cloud points.

[jesusmontera]

I still have a doubt.Does spin vector rotates? In this 2 videos in OpenGL, both spinnors are initalized from a kinetic vector kx=19, ky= 0,kz, so in both cases the spin extracted from psi points between Z and X. (bigger kx makes to point more to X), but in one the spin rotates and in the other doesn't. So witch is the correct?. There's no magnetic or potential acting on them.

[marl0ny]

After glancing at your code, I think I understand now how you want to visualize spin. I don't think summing over the components in the manner that's currently being done will work. I think instead you should compute the total spin expectation value s, where its x component can be computed as sₓ = <ψ| σₓ |ψ> / <ψ|ψ>, where σₓ is the Pauli matrix in the x direction and |ψ> is the wave function. In Python this can be transcribed as sx = np.sum(np.conj(psi)*np.einsum('ij,j...->i...', sigma_x, psi))/np.sum(np.abs(psi)**2). The y and z components are similarly computed. Note that since the wave function has four components the sigma matrices are of course 4 x 4 matrices, where I think you can construct them by putting their 2 by 2 equivalents on the diagonal and zero everywhere else. In Python this may be written as sigma2_x = np.array([[0.0, 1.0], [1.0, 0.0]]); sigma4_x = np.block([[sigma2_x, np.zeros([2, 2])],[np.zeros([2, 2]), sigma2_x]]).

Also I didn’t intend kx or ky to mean kinetic. It’s just a variable name to represent the wave number in the x or y direction (i.e. the spatial frequency with respect to spatial extent of the simulation). I should have named them something like nx or ny instead to avoid confusion.

[jesusmontera]

Thanks, I've changed the code and now the specific spin works fine normalizing bloch [sx,sy,sz], and it does not rotate. But the angle in the np.dot between B and spin vectors influence how big the deviation is, and that does not happens on Stern Gerlach, so I suppose I must set spin to up or down assuming that the magnetic field change the initial spin aligning it. So a spin that have a little up will always go up(if points to B North) and the same for down.

[marl0ny]

For the Dirac solver, I think the only way to get magnetic fields to work properly is by setting the vector potential terms. In the Dirac equation the wave function is explicitly coupled with the vector potential A and not B, and you can just use the relation B = curl A to find what the magnetic field is in terms of the vector potential.

On the other hand for the two component Schrodinger or Pauli equation using the split step method, referring to a previous comment the term exp(i q σ• B(r) Δt) is where the magnetic field interaction occurs. Then it is just applying the matrix exp(i q σ• B(r) Δt) to the wave function at each step.

[jesusmontera]

Thanks, I've applied that B operation in Dirac and Schrod. The particle moves ok, but the spin does not change from it's initial value. And don't know what B = curl A means.
Now I'm trying two Dirac particles in 4D by making two 2d gaussians with specific spin and merging then to 4D with psi4d=np.einsum('jq,ki->ijkq', psi2dA,psi2dB]). Then they are evolved with your DiracSplit class U. But while this works on Schrod in Dirac not at all after U step (the spinors?).
Also I think that without potential(coulomb) the particles will not interact and get entangled, and my purpose is to see what happens on entaglement when I apply a magnetic field only to one. Will the other particle change it's position position too?

[marl0ny]

B is the magnetic field, curl is curl, and A is the magnetic vector potential, sorry if I wasn’t clear. The vector potential can be set through the vector_potential parameter in the initializer, although I think there are currently division by zero errors if there are any zero values in this input (just change these values to something non-zero but tiny). I think working with the vector potential terms is the only way to properly get magnetic fields working when using the Dirac split step method. I only intended the Dirac solver to work for single particles: it will not work properly for multiple particles.

For the expression exp(i q σ• B(r) Δt) mentioned previously, and using equation (2) from here, this can be expressed in Python as exp_magnetic = np.array([[np.cos(a) + 1.0j*nz*np.sin(a), 1.0j*(nx - 1.0j*ny)*np.sin(a)], [1.0j*(nx + 1.0j*ny)*np.sin(a), np.cos(a) - 1.0j*nz*np.sin(a)]]), where the variables a = q*np.sqrt(B[0]**2 + B[1]**2 + B[2]**2)*dt and nx, ny, nz = B/np.sqrt(B[0]**2 + B[1]**2 + B[2]**2) were previously defined (note that B is used in a denominator, so watch out for division by zero errors). Applying the magnetic field term to the wave function is just wave_func = np.einsum('ij...,j...->i...', exp_magnetic, wave_func).

[jesusmontera]

Thanks, I've think I've got how to make spinnors for 2 particles in 4d from your equation to simulate a magnetic field. Is just apply Pauli 4x4 pauli spinnors rotations. And each of the 4 psi's are evolved with your splitstep, witch gives the same result.