# Vortex motion equation for flux-flow instability

I am solving the motion equation of an isolated vortex oscillation in a superconductor. I am assuming that the driving force is generated by an RF current oscillating at 1.3 GHz. The main point of this calculation is to see what's the effect of flux-flow instability under RF drive. The code I wrote is reported below.

f = 1.3 10^9; (*Hz*)
\[Omega] = 2. Pi f;
Trf = 1./f;

Tc = 9.25; (*K*)
T = 1.5; (*K*)

\[Mu]0 = 4. Pi 10^-7; (*H/m*)
\[Phi]0 = 2.07 10^-15; (*Wb*)
n = 5.56 10^28; (*m^-3*)
e = 1.6 10^-19; (*C*)
m = 9.1 10^-31; (*kg*)

\[Tau] = 0.5 10^-9; (*s*)
vf = 1.37 10^6; (*m/s*)
Bc2 = 410. 10^-3; (*T*)
l = 100 10^-9; (*m*)
\[Lambda] = 49.5 10^-9; (*m*)
\[Xi] = 28 10^-9; (*m*)
\[Kappa] = \[Lambda]/\[Xi];

\[Gamma] = (\[Phi]0)/(\[Mu]0 \[Lambda] );

v0 = Sqrt[(l vf Sqrt[14 Zeta (1 - T/Tc)])/(3 Pi \[Tau])];

\[Rho] = (m vf)/(n e^2 l);
\[Eta] = (\[Phi]0 Bc2)/\[Rho];

g = Log[\[Kappa]]+0.5+Exp[-0.4-0.8 Log[\[Kappa]]-0.1(Log[\[Kappa]])^2];
\[Epsilon] = (\[Phi]0^2 g)/(4. Pi \[Mu]0 \[Lambda]^2);

Tmax = 3. Trf;
Zmax = 1. 10^-6; (*m*)

B = 100. 10^-3; (*T*)

sol = NDSolve[{
(v0^2 \[Eta])/(v0^2 + D[x[t, z], t]^2)D[x[t, z],t] == \[Epsilon] D[x[t, z], z, z] + \[Gamma] B Cos[\[Omega] t] Exp[-(z/\[Lambda])],
x[0., z] == 0.,
x[t, Zmax] == 0.,
(D[x[t, z], z] /. z -> 0) == 0.
}, x, {t, 0., Tmax}, {z, 0., Zmax},
AccuracyGoal -> 13,
PrecisionGoal -> 2,
MaxSteps -> Infinity];
u[t_, z_] = Evaluate[x[t, z] /. sol][];
du[t_, z_] = D[u[t, z],t];

phaseSpace =
ParallelTable[{Re[u[t, 0.]] 10^6, Re[du[t, 0.]] 10^-3}, {t, Trf,
Tmax, Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All]


I am able to get a almost correct result out of this calculation only if I fix the solution domain to 1 micron (Zmax = 1. 10^-6). This is however not correct because the vortex should be able to oscillate more than that inside the bulk, therefore I am literally chopping my solution using a too small domain.

If I fix Zmax to a larger number (e.g. 5 um) I need to fix MaxStepSize to a small number in order to avoid errors and the calculation would take more than 2 days on my PC (I actually stopped it after 2 days because the PC froze up!). Does anybody know how to solve this problem?

Also, I was wondering if there is any way to solve the equation when the vortex velocity is higher than v0 (the velocity onset for flux-flow instability). If I fix B to let's say 200 mT, then the vortex velocity would be larger than v0 and the calculation does not return a correct solution, just a series of spikes. Is there any Method I can use to solve this problem?

Mattia

• I posted an answer, but it is probably a numerical artifact. I noticed that all of your BCs and ICs are zero. Is that correct? Jun 13, 2019 at 10:55
• Yes it is correct. In the bulk the vortex do not move and it is perpendicular to the surface. However, the initial condition might be set differently I believe. The best condition would be to put the IC for a vortex with zero velocity, but I do not know the analytical form of it, that's why I put zero as IC (static vortex). I figured that if the solution is calculated for a large enough number of oscillations, it converges to a steady state solution anyway. Jun 17, 2019 at 15:46

When working with small spatial and temporal scales, it is often conducive to scale the variables by characteristic lengths and times (or reciprocal frequencies). The natural choices for length is $$\lambda$$ and reciprocal time is $$\omega$$. If we create dimensionless $$x$$, $$z$$, and $$t$$, and do the appropriate substitutions, we should be able to derive a dimensionless differential equation that depends on several dimensionless groups.

My non-dimensionalized version of your dimensioned equation is

$$\frac{{\frac{{\partial x}}{{\partial t}}}}{{{\nu ^2}{{\left( {\frac{{\partial x}}{{\partial t}}} \right)}^2} + 1}} = \delta \left( {\rho \frac{{{\partial ^2}x}}{{\partial {z^2}}} + Cos(t)Exp( - z)} \right)$$

If we use your initial parameters:

f = 1.3 10^9;(*Hz*)ω = 2. Pi f;
Trf = 1./f;

Tc = 9.25;(*K*)
T = 1.5;(*K*)
μ0 = 4. Pi 10^-7;(*H/m*)
ϕ0 = 2.07 10^-15;(*Wb*)
n = 5.56 10^28;(*m^-3*)
e = 1.6 10^-19;(*C*)
m = 9.1 10^-31;(*kg*)
τ = 0.5 10^-9;(*s*)
vf = 1.37 10^6;(*m/s*)
Bc2 = 410. 10^-3;(*T*)
l = 100 10^-9;(*m*)
λ = 49.5 10^-9;(*m*)
ξ = 28 10^-9;(*m*)
κ = λ/ξ;

γ = (ϕ0)/(μ0 λ);

v0 = Sqrt[(l vf Sqrt[14 Zeta (1 - T/Tc)])/(3 Pi τ)];

ρ = (m vf)/(n e^2 l);
η = (ϕ0 Bc2)/ρ;

g = Log[κ] + 0.5 +
Exp[-0.4 - 0.8 Log[κ] - 0.1 (Log[κ])^2];
ϵ = (ϕ0^2 g)/(4. Pi μ0 λ^2);

Tmax = 3. Trf;
Zmax = 1. 10^-6;(*m*)
B = 100. 10^-3(*T*)


Here is my modified non-dimensionalized version of your code.

Zmax = 100;
Tmax = 20;
vc = λ ω;
nu = vc/v0;
rho = ϵ/(λ B γ);
delta = (B γ)/(η vc);
sol = NDSolve[{D[x[t, z], t]/(1 + nu^2 D[x[t, z], t]^2) ==
delta ( rho D[x[t, z], z, z] + Cos[ t] Exp[-z]), x[0., z] == 0.,
x[t, Zmax] == 0., (D[x[t, z], z] /. z -> 0) == 0.},
x, {t, 0., Tmax}, {z, 0., Zmax}];
u[t_, z_] = Evaluate[x[t, z] /. sol][];
du[t_, z_] = D[u[t, z], t];
phaseSpace =
Table[{Re[u[t, 0.]] , Re[du[t, 0.]] }, {t, Trf ω, Tmax,
Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All] It took about 2.5 seconds to run on my machine. I won't have time to double check my work for a couple days, but I hope this gets you started.

• I received two messages: NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 2.888363945453023*^-45. NDSolve::eerr: Warning: scaled local spatial error estimate of 576279.2946597305 at t = 2.888363945453023*^-45 in the direction of independent variable z is much greater than the prescribed error tolerance.  Jun 13, 2019 at 4:04
• You are correct. I did not hear the beeps last night, but I hear them now and I see the same message. More investigation is needed. Jun 13, 2019 at 10:10
• @AlexTrounev I noticed that all the BCs and ICs are zero, which seems like it would be difficult/impossible to evolve. I'll ask the OP for additional guidance. Jun 13, 2019 at 10:51
• Thank you for the suggestion to make everything dimensionless, I believe it helped a lot on speeding up the calculation. I get the same error Alex Trounev get, but the code still spits out the solution. Jun 17, 2019 at 16:07

A small code change allows you to solve both problems. Figure 1 shows solutions for increasing B and Zmax.

f =
1.3 10^9;(*Hz*)\[Omega] = 2. Pi f;
Trf = 1./f;

Tc = 9.25;(*K*)T = 1.5;(*K*)\[Mu]0 = 4. Pi 10^-7;(*H/m*)\[Phi]0 =
2.07 10^-15;(*Wb*)n = 5.56 10^28;(*m^-3*)e = 1.6 10^-19;(*C*)m =
9.1 10^-31;(*kg*)\[Tau] = 0.5 10^-9;(*s*)vf = 1.37 10^6;(*m/s*)Bc2 =
410. 10^-3;(*T*)l = 100 10^-9;(*m*)\[Lambda] =
49.5 10^-9;(*m*)\[Xi] = 28 10^-9;(*m*)\[Kappa] = \[Lambda]/\[Xi];

\[Gamma] = (\[Phi]0)/(\[Mu]0 \[Lambda]);

v0 = Sqrt[(l vf Sqrt[14 Zeta (1 - T/Tc)])/(3 Pi \[Tau])];

\[Rho] = (m vf)/(n e^2 l);
\[Eta] = (\[Phi]0 Bc2)/\[Rho];

g = Log[\[Kappa]] + 0.5 +
Exp[-0.4 - 0.8 Log[\[Kappa]] - 0.1 (Log[\[Kappa]])^2];
\[Epsilon] = (\[Phi]0^2 g)/(4. Pi \[Mu]0 \[Lambda]^2);

Tmax = 3. Trf;
Zmax = 5. 10^-6;(*m*)B = 300. 10^-3;(*T*)sol =
NDSolveValue[{ \[Eta]/(1 + D[x[t, z], t]^2/v0^2) D[x[t, z],
t] == \[Epsilon] D[x[t, z], z,
z] + \[Gamma] B Cos[\[Omega] t] Exp[-(z/\[Lambda])],
x[0., z] == 0.,
x[t, Zmax] == 0., (D[x[t, z], z] /. z -> 0) ==
0.}, {Re[x[t, 0]]*10^6, Re[D[x[t, 0], t]]*10^-3}, {t, 0.,
Tmax}, {z, 0., Zmax}, AccuracyGoal -> 13, PrecisionGoal -> 2,
MaxSteps -> Infinity] // Quiet;

phaseSpace = ParallelTable[sol, {t, Trf, Tmax, Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All,
PlotLabel -> Grid[{{"Zmax = ", Zmax}, {"B = ", B}}]] • Thank you Alex. There is something strange in the solutions you get, as a function of B both speed (y-axis) and displacement (x-axis) should increase as a function of the rf field amplitude B. I noticed as well that if Zmax is increased (keeping the same B) the code tends to return a solution that qualitatively is ok (an ellipse) but it is not accurate and generally underestimated. To partially solve this I used the "SingleStep" method specifying MaxStepSize and StartingStepSize. The solution in this way is more accurate when Zmax is increased. Jun 17, 2019 at 16:11
• I meant "FixedStep", sorry Jun 17, 2019 at 16:17
• @Mattia OK! You probably understand what the solution should be. I just indicated how to choose a branch of solution using  NDSolveValue instead of  NDSolve`. Jun 17, 2019 at 16:36