Trouble with a custom EquationOfMotion class

Geometry, Fields and Transportation
1

We are running a solid-state physics simulation, including charge carriers with an electric field bias. We have written our own equation-of-motion class (G4CMPEqEMField) as a subclass of G4EqMagElectricField.

I’m having trouble reconciling the sparse guidance in the Toolkit Guide with the code in G4EqMagElectricField::EvaluateRhsGivenB(), which we used as a guide for our own code. The documentation says

the derivative vector dydx is a vector with six components. The first three are the derivative of the position with respect to the curve length. Thus they should set equal to the normalised velocity, which is components 3, 4 and 5 of y.
(dydx[0], dydx[1], dydx[2]) = (y[3], y[4], y[5])

But the code assumes that y[3,4,5] is the actual momentum (in Geant4 energy units), and computes the output components above as

   G4double pSquared = y[3]*y[3] + y[4]*y[4] + y[5]*y[5] ;
   G4double pModuleInverse  = 1.0/std::sqrt(pSquared) ;
   dydx[0] = y[3]*pModuleInverse ;                         
   dydx[1] = y[4]*pModuleInverse ;                         
   dydx[2] = y[5]*pModuleInverse ;

So dydx[0,1,2] does come out as a unit (direction) vector. The documentation also doesn’t mention the dydx[6]=0. and dydx[7]=inverse_velocity; components.

Besides the documentation, I’m also having a physics problem. I would expect the field to change the particle’s kinetic energy by an amount equal to the voltage drop (field * path length). This is what happens with G4EqMagElectricField, but not with our code. Instead, I seem to get systematically too much energy (by about 25% or so).

Unfortunately, the “stepper” code is sufficiently abstract that I really can’t tell how it is that the particle’s momentum and energy get changed by the output of EvaluateRhsGivenB(). Is there somewhere (slides from a talk, maybe?) that would provide a physically motivated connection that I could follow?

1 Like
2

@japost John, is this something you could help with, or point me toward old talks?

The specific issue I’m fighting with is we have non-trivial kinematics. In particular, anisotropic electron transport and differences between momentum and velocity vectors. I have confirmed that we are using correct coordinate transforms to rotate the electric field vector so that it applies the correct anisotropic force on the electron, but we are losing energy conservation in the process. The electron’s change in kinetic energy (after transport) does not match the voltage drop between the pre- and post-step position.

1 Like
3

Hi Mike.

In the original version of the integration of motion for field propagation, it was the (unit) direction which was integrated in components [3], [4] and [5]. However as colleague(s) at FNAL correctly pointed out, it is numerically more stable to integrate the momentum - so we adopted this. If you find remnants of the old convention in comments or data member names , I/we will need to correct them. ( There should be no remnants in the code. )

My first suggestion is to control the error of integration - set Epsilon Min/Max to the same small value, say 1.0e-6.

The higher order integration methods are a bit difficult to parse upon the first ‘season’ of exposure to them. But you could use a lower order method ( e.g. Euler, SimpleRunge 2nd order, Simple Heum 3rd order) if you want to gain a feeling how a lower order method would fare, and intuition.

And you could check whether a different method (e.g. G4CashKarpRKF45 or G4BogackiShampine) might give slightly different results.

But my guess is that there is a problem elsewhere.

From your description I am having difficulty understanding how your equation of motion differs from G4EqMagElectricField.

I would suggest also that you/we make (experimental) changes in the methods that load the integration vector, to make y[6] = kineticEnergy y[8] = Voltage and make corresponding changes to the equation of motion. That way we should be able to confirm whether it is the integration that is causing any problem, or something else.

4

One key item is that G4EqMagElectricField calculated the velocity of the particle from it’s momentum and applies

d/ds P-vector = 1/velocity * d/dt P-vector

Is this where your equation differs ?

5

Thanks for the replies, John.

No. That’s what I calculate, specifically

force = Efield * fCharge*vinv*c_light;

where “G4threeVector force” goes into the dP/dx components dydz[3,4,5].

What makes our equation of motion different is that we transform Efield (initially Efield = G4ThreeVector(field[3],field[4],field[5])) into an anisotropic coordinate system corresponding to the electron charge carrier’s Herring-Vogt frame. That is a coordinate frame in the Brillouin zone’s valley in which gives the electron an isotropic potential energy surface.

  fGlobalToLocal.ApplyAxisTransform(Efield);
  fLocalToLattice.ApplyAxisTransform(Efield);
  fLatticeToValley[i].ApplyAxisTransform[Efield);
  Efield *= theLattice->GetSqrtInvTensor();

followed by the inverse transforms to get back to world coordinates. The “SqrtInvTensor()” is a 3x3 matrix diagonal matrix

Diag(1/sqrt(m(xx)/mHV), 1/sqrt(m(yy)/mHV), 1/sqrt(m(zz))/mHV)

normalized to unity by mHV = 3/[1/m(xx) + 1/m(yy) + 1/m(zz)]. The three masses correspond to the electron’s anisotropic effective mass along (m(xx)) and across the valley direction. This is an approximation, but it’s one we use consistently throughout the package.

In the Geant4 kinematics (G4DynamicParticle), we set the kinetic energy and momentum direction of the electron, and set the mass equal to the electron’s effective mass along the momentum direction (m(eff) = p^2 / 2Ekin). This is done by way of a step limiter process, to ensure that the kinematics stay accurate.

In your first reply, you wrote,

Okay, so y[3,4,5] is the momentum, and we should set dydx[0,1,2] = p.unit(), and dydx[7] as the inverse velocity (m/p with the factors of c and c^2 taken out). That’s what we are doing now, so that is reassuring.

I’d like to try that, but I’m having trouble finding where to set those. In our field manager, we set stepperLength = 1e-9*mm (which may be too small). I’ve looked at the .hh files for all the transport classes (chord finder, driver and stepper) but I don’t see where to set the integration error bounds.

6

I think that a step length of 1.0e-9 mm is too small. Also it is not obvious that something like this is needed.

You can find the recommendations for revising the integration accuracy (and other ways to control the step) in the Applications developer guide section on adjusting integration accuracy for external EM fields

The epsilon max value is the most important, as it defines the relative accuracy for the smallest steps.

G4double minEps= 1.0e-6;  //   Minimum & value for largest steps
G4double maxEps= 1.0e-5;  //   Maximum & value for smallest steps

fieldManager->SetMinimumEpsilonStep( minEps );
fieldManager->SetMaximumEpsilonStep( maxEps );

You can choose to set the same value - this will be the acceptable relative accuracy for all steps.

Otherwise you can provide a length scale for the biggest acceptable error in a single step:

globalFieldManager->SetDeltaOneStep( 0.1 * CLHEP::micron );

The relative errors must still remain within the bounds given, so this value will be used only when it results in a relative error within those bounds.

7

Thank you! I had been looking in all the individual classes; I should have looked at G4FieldManager itself. I tried setting the limits as you described, but they don’t change my results.

I have done two things which are pointing toward some solution, but I don’t yet understand it.

  1. I set the charge carrier mass matrix to be the identity matrix (i.e., the same value for m(xx), m(yy), m(zz). That means the electron behaves like a normal particle with a scalar mass, and our custom equation of motion follows that. All the transforms still happen, but there’s no alteration of the electric field components. We end up with the “force” vector (qEc/v) pointed along the field direction. In this configuration, there’s no discrepancy between the voltage drop and the electron’s energy gain.

  2. I inserted a rotation to an off-axis force (qEc/v doesn’t point along the E field), for normal particles. No unusual kinematics, just point the force (dydx[3,4,5]) in a different direction than the input electric field (field[3,4,5]).

The second case resulted in the same discrepancy between voltage drop and energy gain. The momentum change seems correct (at least, the direction of momentum change was along the force direction), but the corresponding kinetic energy change is not. In some steps, even the sign is wrong (the track lost energy while gaining momentum).

Does the stepper assume that the force and field directions must align? I wouldn’t think so, because that’s not true for magnetic fields, nor for combined E and B fields.

8

All general steppers are agnostic about what is being integrated. These include the older ones G4ClassicalRK4 and G4CashKarpRKF45, and the newer ones such as G4DormandPrince745.

9

Okay, thanks, John!

To summarize my situation, I find that if I have a pure electric field and I compute the force in EquationOfMotion such that it points in a different direction, then G4Transportation no longer conserves energy (the energy of the track after each step doesn’t correspond to the voltage drop along that step).

I think this makes sense, since the rotated force means the position change has a component that doesn’t lie along the field direction. So the computed voltage drop is different from the energy gained from the force.

I have tried rescaling the force vector to get either average energy conservation (using a hardcoded number), or even step-by-step energy conservation, but it’s not working right.

10

@japost Maybe I should open a new forum topic about this; I want to reframe my original question slightly differently, as I’ve thought of a better way to think about what our needs are.

How do I implement a “composite” field and equation of motion class? My situation is such that I have two effective force components, one component which changes the particle’s energy (the electric field) and a second component which affects the particle’s direction of motion, but doesn’t change its energy.

How should I implement this in the EqOfMotion class? I looked at G4EqMagElectricField, which does the qE + qvxB, but I don’t see where the kinetic energy change comes in. I also don’t see what do with the field class.

11

Okay, I dove into the G4Transportation, G4FieldTrack, G4PropagatorInField maze and I think I understand better. It appears that the propagator directly integrates the dydx[] steps supplied by EqOfMotion::EvaluateRHSGivenB(), and it then passes the integral result as an array to G4FieldTrack::LoadFromArray(). In that function, valArr[3,4,5] is the final momentum, and that function computes fKineticEnergy = momentum (+) mass.

I think I can answer my own question above. There are two separate forces in play: the electric field force should do what it does, and change the electron’s momentum (and consequently energy) appropriately. There’s a separate force which merely acts to keep the electron moving along the valley, essentially rotating it’s momentum direction without chaing the magnitude.

I think the proper way to treat this in Geant4 is to define two different fields, one a normal electric field, and the other a “magnetic like” field that know how to rotate the momentum correctly without changing its direction. That’s true if you can have multiple fields registered for a volume.

As an alternative, I can rescale the rotated force in the code I have now, such that its total magnitude corresponds to the portion of the rotated magnitude along the true E-field direction. Then it’ll give the electron a kick along the rotated direction, but it’ll change the electron’s momentum magnitude (and hence energy) by an amount corresponding to the true voltage drop.

12

I agree that you need to create a (single) custom equation of motion which combines the ‘external’ electric field and the force which the particle feels from the differential effective mass.

The ‘external’ E-field can be a parameter to this equation of motion, which it can consult or integrated into it (less general but maybe faster.)