Let's run some tests and see what the dynamics of the defect density looks like.

First, start and stay in the ferromagnetic phase.

The defect density is indeed zero; similarly the system always stays in the instantaenous ground state, since the Hamiltonian is constant.

Now let's try constant fields in the paramagnetic phase, for which $g>g_c=1$.

The defect density stays at a constant value, which is determined by the strength of $g$. Higher and higher values of $g$, taking the system deeper into the paramagnetic phase, should cause the defect density to saturate at 0.5 (although it doesn't really make sense to talk about topological defects in the paramagnetic phase, as the symmetry is not broken; maybe kink density is better).

Looks good.

Now, we will use the linear quench. Let's fix the chain size at $N=10$ and see if we can access different regimes for the dynamics by changing the rate at which we quench.

The quench occurs from time $t=0$ to time $t=T$, from $g=2$ to $g=0$. So $g(t)=2-\frac{2}{T}t$. This means that the quench time is $\tau_Q=T/2$.

From Dziarmaga: using Landau-Zener one can show that if $\tau_Q<\tau_Q^{ad}$, then the quench becomes nonadiabatic, where $\tau_Q^{ad}=\frac{\hbar N^2}{2\pi^3 J}\approx 0.016 N^2$. So by varying $T$ on either side of this threshhold, we should observe adiabatic behavior, where the final defect density is close to 0 and the fidelity remains high (indicating that the system remains close to its instantaneous ground state), and nonadiabatic (diabatic?) behavior, where the defect density does something else, and the fidelity drops. This nonadiabatic regime is what we are interested in, since we want to control the defect density.

For the largest quench times, $\tau_Q=25,10$ or $T=50,20$, $D$ decreases all the way to 0 in the ferromagnetic phase. Correspondingly, the fidelities for these two quench times are very close to 1.0 at all times, except for a slight dip for the $\tau_Q=10$ curve when $g$ is close to its critical value.

Clearly, as the quench time is reduced, the dynamics becomes less adiabatic. The fidelities for smaller $\tau_Q$ never recover once the critical point is passed at $2t/\tau_Q=0.5$, and $D(T)$ starts to rise. Once $\tau_Q<\tau_Q^{ad}$, markedly different behavior appears, in that $f(T)$ becomes much lower and $D(T)$ becomes much higher. Indeed, for $\tau_Q=0.05$, the defect density appears to be constant, indicating that the system is "frozen out", and that it cannot respond to a change in the magnetic field. This means that we probably would not be able to control the system on such a short time scale either, but we can still try.

Of course, the above plots are for $N=10$. Different sizes will yield different behaviors of $D$ and $f$ for these values of $\tau_Q$, but if the quench times are scaled according to $\tau_Q \sim N^2$, then the same behavior should be observed. Let's check this. Double $N$, and quadruple all the quench times above.

The agreement is not perfect, but somewhat close. Perhaps this is because although the analysis in deriving $\tau_Q$ makes sense for finite $N$, it also assumes that $N$ is large. Anyway, this suggests that larger systems are more responsive to changes in $g$ and therefore more amenable to being controlled over a greater variety of time scales. We will try to simulate as large systems as possible, since the computational expense of the numerics scales linearly with $N$ (but also linearly with $T$, which we may also increase as the system size gets bigger).

If we want to increase the number of defects at the final time, we can try quenching from a higher initial value of $g$, so that the system starts deeper in the paramagnetic phase. We still have to fix $g(T)=0$ ( $g<-1$ again belongs to the paramagnetic phase), so if we do this we will have to adjust $T$ or $\tau_Q$ accordingly. Let's check to see how starting at a higher value of $g(0)$ affects the defect density. We will fix $N=10$ and $T=1$, so that the dynamics is nonadiabatic.

We see that the defect density gets higher, and that the fidelities change at later and later times; this is because the critical point is being crossed later. We could try changing the pulses we optimize in this way later, but this change may be identical to just making $\tau_Q$ shorter.

Now we will check that the code reproduces the KZ scaling, which together with Landau-Zener states that the density of defects should scale with the quench time as $D=\frac{1}{2\pi}\sqrt{\frac{\hbar}{2J\tau_Q}}$ (Zurek and Dziarmaga). In numerical simulations, the exponent of $\tau_Q$ was found to decrease with $N$, down to 0.58 for $N=100$. We will use smaller $N$ for now, so if we obtain an exponent higher than 0.5 we should not be surprised. We will fix $N=50$, $g(0)=5$ (to match them), and $g(T)=0$. Zurek et al. plot the defect density against the dimensionless quench rate $\tau_0/\tau_Q$, where $\tau_0=\hbar/2J$. For us, then, $\tau_0=1/2$ and $\tau_Q=T/5$, so $\tau_0/\tau_Q=5/2T$. We will calculate $D(T)$ for $\tau_0/\tau_Q=0.01,$ $0.1,$ $0.2,$ $0.3,$ and $0.4$. To compare with the following figure from Zurek et al.:

With these few points, we can reproduce Zurek et al.'s results. Our fited exponent is 0.59, and our coefficient is 0.19, both of which are in excellent agreement with their results of 0.58 and 0.16 using larger $N$ and more points.

Now, let's see what D-MORPH can do by looking at a (partially) optimized pulse. For this test simulation we set $N=20$, $T=1$, and used an an initial control a linear quench such that $g(0)=2$ and $g(T)=0$. The optimization lasts until $s$, the DMORPH parameter, reaches a pre-specified value; we chose $s_{max}=1$ for now. First, we will examine the initial and final control fields $g(t)$:

Now we will examine the resulting defect densities and fidelities as functions of time.

We have managed to lower the final defect density by a little. Letting DMORPH run longer should produce a better result.