We develop a theoretical framework that allows us to explore the coupled motion of neutron-superfluid vortices and proton-superconductor flux tubes in a gravitationally collapsed condensate, which describe neutron stars that form pulsars. Our framework uses the 3D Gross-Pitaevskii-Poisson-Equation (GPPE) for neutron Cooper pairs, the Real-Time-Ginzburg-Landau equation (RTGLE) for proton Cooper pairs, the Maxwell equations for the vector potential ${\bf A}$, and Newtonian gravity and interactions, both direct and induced by the Poisson equation, between the neutron and proton subsystems. For a pulsar we include a crust potential, characterized by an angle $\theta$, and frictional drag. By carrying out extensive direct numerical simulations of this model, we obtain a variety of interesting results. We show that a rotating proton superconductor generates a uniform London magnetic field, which changes the field distribution inside flux tubes. In the absence of any direct interaction between the two species, they interact through the gravitational Poisson equation. The presence of attractive (repulsive) density-density interaction leads to the attraction (repulsion) between neutron vortices and proton flux tubes. The inclusion of the current-current interaction and the complete Maxwell equations allows us to quantify the entrainment effect that leads to induced magnetization of neutron vortices. We demonstrate that, with a strong external magnetic field ${\bf B}_{\rm ext}$, proton flux tubes are anchored to the crust, whereas neutron vortices leave the condensate and lead to abrupt changes of the crust angular momentum ${\rm J}_c$. The frictional term in the dynamical equation for $\theta$ yields stick-slip dynamics that leads, in turn, to glitches in the time series of ${\rm J}_c$. By calculating various statistical properties of this time series, we demonstrate that they display self-organised criticality (SOC) that has been found experimentally for several pulsars. We compare our results with those of earlier explorations of pulsar-glitch statistics in GPE-based minimal models for pulsars.

Neutron-superfluid vortices and proton-superconductor flux tubes: Development of a minimal model for pulsar glitches
Sanjay Shukla
Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India
Marc E. Brachet
Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université Université de Paris, 24 Rue Lhomond, 75005 Paris, France

Rahul Pandit ^1

§ INTRODUCTION

Recent advances in the Gross-Pitaevskii-Poisson (GPP) modelling of bosonic <cit.> and axionic <cit.> stars have led to an elegant, minimal model for pulsars <cit.>, which includes a crust potential, and leads naturally to pulsar glitches <cit.>. These GPP models have, so far, accounted only for bosons, e.g., the neutron superfluid in a neutron star <cit.>. Within the outer core of a neutron star, characterized by a density ranging from $\rho \simeq10^{13}$ to $10^{14} \ \rm g \ cm^{-3}$, lies a region consisting predominantly of neutrons ($95\%$ of the total mass) and some protons ($5\%$ of the total mass), and sufficient electrons to maintain charge neutrality; this region exhibits extraordinary properties: the neutron Cooper pairs form a superfluid <cit.> and the proton Cooper pairs a superconductor <cit.>.
Pulsars are rapidly rotating and highly magnetized neutron stars <cit.> with magnetic fields $\simeq 10^{12} \ \rm G$. Our goal is to generalise the GPP modelling framework for pulsars <cit.> to include protons, which are in a superconducting state that is affected strongly by the magnetic field.

Neutrons in a neutron-star interiors are in a superfluid state, so when the star rotates with a sufficiently large angular velocity, quantized vortices are formed; each vortex has an angular momentum that is an integer multiple of the quantum of circulation $\mathcal{K}=\frac{h}{m_n^*}$, where $m_n^*$, the mass of a neutron Cooper pair,  is twice the mass of a neutron. By contrast, the protons in a neutron star form an Abrikosov phase of a Type II superconductor <cit.>, in which the external magnetic field leads to an array of flux tubes, each carrying a magnetic flux quantum $\Phi_0=\frac{hc}{q}$, with $q=2e$ being the charge of a proton Cooper pair. Early studies by Ginzburg and Kirzhnits <cit.>, Wolf <cit.>, and Baym, Pethick, and Pines <cit.> laid the foundations for our understanding of neutron superfluidity and proton superconductivity in neutron stars.

Pulsars exhibit sudden increases, known as glitches, in their rotational frequencies <cit.>. The interaction of the pulsar crust with the neutron superfluid may provide an explanation for these glitches, as first suggested by Baym, Pethick, and Pines <cit.>, and as explored
recently by our group in Ref. <cit.>. Pulsar-glitch observations <cit.> suggest that there is a connection of glitches with the transfer of angular momentum, stored in the quantized vortices of the neutron superfluid, to the solid crust of a pulsar. Various models, such as those based on vortex avalanches <cit.> or superfluid vortex-crust interaction <cit.>, have been proposed to study the glitching phenomenon.

As we have noted above, a neutron superfluid dominates the outer core of a pulsars.  Therefore, Refs. <cit.> have utilized a simple model for a pulsar in which the Gross-Pitaevskii equation (GPE) is used for the neutron superfluid together with a pinning potential for the crust and a rotating container that is defined by a quadratic confining potential. In Ref. <cit.>, our group has removed the confining potential but introduced Newtonian gravity, which leads to a gravitationally collapsed bosonic pulsar that displays glitches whose statistical properties are akin to those seen in several pulsars.

Even though protons constitute only about $\simeq 5 \%$ of the mass of a pulsar, they play a crucial role in its rich dynamics because the strong magnetic field leads to the formation of an Abrikosov flux lattice <cit.>. The flux tubes in this lattice can interact with the vortices in the neutron superfluid <cit.>. References <cit.> and <cit.> have examined such interactions, both in equilibrium and out-of-equilibrium conditions, in the 2D and 3D GPE systems, but with (a) a static Ansatz for the proton flux tubes and (b) harmonic confinement. We go beyond approximation (a) and replace the harmonic potential in (b) by Newtonian gravity that leads to a gravitationally collapsed pulsar. In particular, we develop a theoretical framework by combining the Maxwell equations for the electromagnetic fields with the 3D Gross-Pitaevskii-Poisson-Equation (GPPE) for the neutron superfluid and the Real-Time-Ginzburg-Landau equation (RTGLE) for the proton-superconductor system. Moreover, we include (i) density-density and (ii) current-current direct interaction between neutron and proton Cooper pairs.

Before going into the details of our calculations, we present the principal results of our study of the coupled GPPE, RTGLE, and Maxwell systems:

* We show that the rotation of the proton superconductor leads to a London moment <cit.>, i.e., inside the superconductor there is a uniform magnetic field, whose magnitude depends on the rotation frequency.

* We evolve the magnetic field of the proton flux tubes by using the Maxwell equations; this leads to a precise characterization of the entrainment of proton-superconductor flux tubes around neutron-superfluid vortices, by vitrue of which these vortices also become magnetised.

* We show that, even if there is no direct interaction between the neutron and proton Cooper pairs, they interact through the gravitational Poisson equation, which leads to an attraction between neutron-superfluid vortices and proton-superconductor flux tubes.

*  Direct attractive (repulsive) density-density interactions cause the cores of neutron-superfluid vortices and proton-superconductor flux tubes to attract (repel). Furthermore, the current-current interaction, which represents the interaction between the momenta of neutron and proton Cooper pairs, leads to an entrainment effect by incorporating the entrained proton mass current into the Maxwell equation for the vector potential. This is the first demonstration of such entrainment in the GPPE context.

* If $\Theta > 0$ is the initial angle between the rotation axis and the external magnetic field, then, eventually, the proton-superconductor flux tubes tend to align themselves along the rotation axis. We demonstrate this alignment by calculating the magnetic moment of the proton Cooper pairs.

* We follow the real-time dynamics of the GPPE and RTGLE together with the crust potential, for the illustrative case  $\Theta=0$ and with the neutron and proton Cooper pairs interacting only via the gravitational potential. This gives rise to a collapsed pulsar, with a crust angular momentum that
displays glitches with signatures of self-organized criticality (SOC) <cit.>.

The remainder of this paper is organised as follows: In Section <ref>, we provide a comprehensive description of the GPPE and RTGLE models. Section <ref> outlines the units and dimensionless forms of GPPE and RTGLE, accompanied by an elucidation of the pseudospectral method employed for our study. Our results are presented in Section <ref>, followed by a discussion of conclusions in Section <ref>.

§ THE MODEL

The total Lagrangian $\mathcal{L}$ governing the dynamics of neutron and proton Cooper pairs within the system is composed of distinct Lagrangians. Section <ref> delves into $\mathcal{L}_{\rm n}$, which encapsulates the dynamics of neutron Cooper pairs. Similarly, Section <ref> focuses on $\mathcal{L}_{\rm p}$ describing the proton Cooper pairs, while the electromagnetic field is described by $\mathcal{L}_{\rm EM}$ (Section <ref>). The interactions between neutron and proton Cooper pairs are addressed through $\mathcal{L}_{\rm np}$ in Section <ref>. Finally, we present the governing equations of motion in Section <ref> by using the total Lagrangian $\mathcal{L}$.

§.§ Neutron Superfluid

In a pulsar, neutrons form Cooper pairs that lead to a superfluid <cit.>. At temperatures below the transition temperature $T_{\lambda}$, these Cooper pairs lead to a Bose-Einstein condensate (BEC), characterized by a macroscopic complex wavefunction $\psi_n$. The Lagrangian of a weakly interacting rotating BEC in a self-gravitating potential $\Phi$ is given by <cit.>:

\begin{eqnarray}
\mathcal{L}_{\rm n} &=& \frac{i\hbar}{2} \left( \psi_n^* \frac{\partial \psi_n}{\partial t}- \psi_n \frac{\partial \psi_n^*}{\partial t}\right) -\frac{\hbar^2}{2m_n}|\nabla \psi_n|^2  \nonumber\\
&-&\frac{g}{2} \bigg(|\psi_n|^2-\frac{\mu_n}{g}\bigg)^2 -m_n\Phi|\psi_n|^2- \frac{1}{8\pi G} (\nabla \Phi)^2 \nonumber \\
&-& V_{\theta} |\psi_n|^2+\frac{i\hbar}{2}({\bf\Omega} \times {\bf r}) \cdot (\psi_n\nabla \psi_n^* - \psi_n^*\nabla \psi_n)\,,
\label{eq:Lag_n}
\end{eqnarray}

where $g=4\pi a \hbar^2/m_n$ is the interaction strength between neutron Cooper pairs, $a$ is the s-wave scattering length, $m_n$ and $\mu_n$ are, respectively,
the mass and chemical potential of neutron Cooper pairs, $G$ is Newton's gravitational constant, $V_{\theta}$ is the crust potential, and $\bf\Omega$ is the rotational velocity.

§.§ Proton Superconductor

Proton Cooper pairs, which also form in a pulsar, yield a Type II superconductor with an Abrikosov flux lattice <cit.>. This superconducting system can be described by the complex wavefunction $\psi_p$, coupled to a vector potential $\bf A$ and self-gravitating potential $\Phi$, resulting in the following Lagrangian, in which we include the rotational velocity $\bf\Omega$:

\begin{equation}
\begin{aligned}
\mathcal{L}_{\rm p} &= \frac{i\hbar}{2} \left( \psi_p^* \frac{\partial \psi_p}{\partial t}- \psi_p \frac{\partial \psi_p^*}{\partial t}\right) -\frac{1}{2m_p}|D_{\bf A}\psi_p|^2 \\
&- \frac{\alpha_s}{2}\bigg(|\psi_p|^2-\frac{\mu_p}{\alpha}\bigg)^2 -q \phi |\psi_p|^2 -m_p\Phi|\psi_p|^2 - V_{\theta} |\psi_p|^2\\
&+\frac{1}{2}({\bf \Omega} \times {\bf r}) \cdot (\psi_pD_{\bf A} \psi_p^* + \psi_p^*D_{\bf A} \psi_p)\,,
\end{aligned}
\label{eq:Lag_p}
\end{equation}

where $D_{\bf A} \equiv [\frac{\hbar}{i}\nabla -q {\bf A}$] is the magnetic gradient operator, $m_p$ is the mass of a proton Cooper pair, $\alpha_s$ is the interaction strength between the proton Cooper pairs, $\mu_p$ is the proton chemical potential, $q=2e$ is the charge of a proton Cooper pair, and $\phi$ is the electric scalar potential.

The evolution of the vector potential ${\bf A}$ follows the Maxwell equations, which can be obtained from the electromagnetic Lagrangian:

\begin{equation}
\begin{aligned}
\mathcal{L}_{\rm EM} &= \epsilon_0 \bigg[ -\frac{1}{2} [{\bf E}^2 - c^2({\bf B}-{\bf B}_{ext})^2] \\
+ {\bf E} &\cdot \bigg(-\nabla \phi - \frac{\partial {\bf A}}{\partial t }\bigg) -c^2({\bf B}-{\bf B}_{\rm ext}) \cdot (\nabla \times {\bf A})\bigg]\,,
\end{aligned}
\end{equation}

where $\bf E$ and $\bf B$ are the electric and magnetic fields, respectively, ${\bf B}_{\text{ext}}$ is the external magnetic field, and $c$ is the speed of light.

§.§ Interaction between Neutron and Proton Cooper Pairs

We consider (i) the density-density and (ii) the current-current interactions between the neutron and proton subsystems <cit.> and use the Lagrangian

\begin{equation}
\begin{aligned}
\mathcal{L}_{\rm np} &= \gamma \bigg\{ g_{np} |\psi_n|^2 |\psi_p|^2\\
&-\frac{\hbar}{4i} ( \psi_n \nabla \psi_n^* - \psi_n^* \nabla \psi_n)\cdot \left[ \psi_p D_{\bf A} \psi_p^* + \psi_p^* D_{\bf A} \psi_p\right] \bigg\}\,,
\label{eq:interaction_lag}
\end{aligned}
\end{equation}

where $\gamma$ is the overall interaction strength, $g_{np}$ is the density-density coupling constant <cit.>. The current-current interaction, given in the second row, causes neutron-superfluid vortices to drag proton-superconductor flux tubes.

§.§ Total Lagrangian: equations of motion

By combining the Lagrangians (<ref>)-(<ref>) we obtain the total Lagrangian

\begin{eqnarray}
\mathcal{L} = \mathcal{L_{\rm n}}+\mathcal{L_{\rm p}}+\mathcal{L_{\rm EM}}+\mathcal{L_{\rm np}}\,.
\end{eqnarray}

The Euler-Lagrange equations for $\mathcal{L}$ yield the following:

* the Gross-Pitaevskii-Poisson equation (GPPE) for neutron Cooper pairs (variation with respect to $\psi_n^*$):

\begin{eqnarray}
i\hbar\frac{\partial \psi_n}{\partial t} &=& -\frac{\hbar^2}{2m_n}\nabla^2 \psi_n-\mu_n \psi_n  + g|\psi_n|^2 \psi_n+m_n\Phi\psi_n\nonumber\\
&+& i\hbar ({\bf \Omega} \times {\bf r}) \cdot \nabla\psi_n+ V_{\theta} \psi_n+\gamma g_{np}|\psi_p|^2 \psi_n\nonumber\\
&-& \frac{\gamma\hbar}{2i} \left[ \nabla \psi_n \cdot {\bf J}_p + \nabla \cdot \left( \psi_n {\bf J}_p\right)\right]\,;
\label{eq:GPE_neutron}
\end{eqnarray}

here, ${\bf J}_p$ is the proton current density:

\begin{eqnarray}
{\bf J}_p &=& \frac{\hbar}{2i}(\psi_p^*\nabla\psi_p - \psi_p\nabla \psi_p^*) - q {\bf A}_{\rm eff}|\psi_p|^2\,;
\label{eq:JnJp}
\end{eqnarray}

* the real-time Ginzburg-Landau-Poisson equation (RTGLPE) for proton Cooper pairs (variation with respect to $\psi_p^*$):

\begin{equation}
\begin{aligned}
i\hbar \left( \frac{\partial }{\partial t} + \frac{i}{\hbar}q\phi_{\rm eff} \right) \psi_p &= \frac{1}{2m_p}\bigg(\frac{\hbar}{i}\nabla- q{\bf A}_{\rm eff}\bigg)^2 \psi_p-\mu_p\psi_p \\
&+\alpha_s |\psi_p|^2\psi_p+m_p \Phi \psi_p+ V_{\theta}\psi_p \\
&+ \gamma g_{np} |\psi_n|^2 \psi_p+\gamma q({\bf J}_n\cdot {\bf A})\psi_p \\
&- \frac{\gamma\hbar}{2i}[{\bf J}_n.\nabla \psi_p + \nabla \cdot({\bf J}_n \psi_p)]\,;
\label{eq:GLE_proton}
\end{aligned}
\end{equation}

\begin{eqnarray}
\phi_{\rm eff} &=& \phi-\frac{m_p}{2q}\Omega^2 r^2\,;\nonumber \\
{\bf A}_{\rm eff} &=& {\bf A}+\frac{m_p}{q}({\bf \Omega} \times {\bf r})\,; \nonumber \\
{\bf J}_n &=& \frac{\hbar}{2i}(\psi_n^*\nabla\psi_n - \psi_n\nabla \psi_n^*)\,;
\label{eq:eff}
\end{eqnarray}

* the Maxwell equation for the vector potential
${\bf A}$ and Poisson equations for the gravitational potential $\Phi$ and the electric scalar potential $\phi$:

\begin{equation}
\begin{aligned}
\frac{1}{c^2}\frac{\partial^2 {\bf A}}{\partial t^2} -\nabla^2 {\bf A}-\nabla \times {\bf B}_{\rm ext}&= \mathbb{P}\bigg[\frac{q}{m_pc^2\epsilon_0} {\bf J}_{p}\\
&-\frac{\gamma q}{c^2\epsilon_0} {\bf J}_n|\psi_p|^2\bigg]\,;
\end{aligned}
\label{eq:vector_pot}
\end{equation}

\begin{equation}
\begin{aligned}
\nabla^2 \Phi &= 4\pi G \bigg(m_n|\psi_n|^2 +m_p|\psi_p|^2 -\rho_{\rm bg}\bigg)\,;
\end{aligned}
\label{eq:grav_pot}
\end{equation}

\begin{eqnarray}
\rho_{bg} = m_n\langle|\psi_n|^2\rangle+m_p\langle|\psi_p|^2\rangle\,;
\label{eq:back}
\end{eqnarray}

\begin{eqnarray}
\nabla^2\phi = -\frac{1}{\epsilon_0} (q|\psi_p|^2-qn_p)\,.
\label{eq:scalar_pot}
\end{eqnarray}

With $\gamma=0$ and without crust potential ($V_{\theta}=0$), Eq. (<ref>) has been used extensively in Refs. <cit.> to study self-gravitating BECs at temperature $T=0$ (by using a Gaussian Ansatz for $|\psi_n|^2$). References <cit.> have performed numerical simulations without rotation (${\bf \Omega}=0$). Furthermore, Ref. <cit.> has included rotation in the GPPE to study the dynamical properties of BEC dark matter (but still with $\gamma=0$ and $V_{\theta}=0$). In our previous studies, we have used the GPPE to study the formation of compact bosonic objects at finite temperatures <cit.> and their axionic counterparts, by including a quintic nonlinearity <cit.>. The imaginary time ($t\to -it$) version of Eq. (<ref>) with $\gamma=0$, $V_{\theta}=0$, and $\Phi=0$ is the well-known time-dependent Ginzburg-Landau equation <cit.> whose solutions give Type I and Type II superconductors.

In writing Eq. (<ref>), we have used the Coulomb gauge $\nabla \cdot {\bf A}=0$. The Helmholtz projector, which has components $\mathbb{P}_{ij}:= \delta_{ij}-\mathcal{F}^{-1} \frac{k_ik_j}{k^2}\mathcal{F}$, with $\mathcal{F}$ the Fourier-transform operator, projects a field onto its divergence-free part; its application in Eq. (<ref>) maintains the Coulomb gauge. In Eq. (<ref>) for $\Phi$, the gravitational potential, the subtraction of the background mean density $\rho_{\rm bg}$ [often called the Jeans Swindle <cit.>]  can be justified either by defining a Newtonian cosmological constant <cit.> or by accounting for cosmological expansion <cit.>. Furthermore, in Eq. (<ref>) for the scalar potential $\phi$, we subtract the mean charge density, $qn_p$, to maintain charge neutrality in the system. In the context of a neutron star, $qn_p$ corresponds to the background charge coming from electrons.

The important term considering the interaction between neutron and proton Cooper pairs is the last term in Eq. (<ref>), which is of the form $\frac{\gamma q}{c^2\epsilon_0} {\bf J}_n|\psi_p|^2$. This term, not considered hitherto in the GPPE and RTGLPE, causes neutron-superfluid vortices to drag proton-superconductor flux tubes, generate an entrained-proton current, because of which the neutron-superfluid vortices become magnetized, as we show below.

A schematic diagram [cf. Refs. <cit.>] of the interior of a pulsar (magnetized neutron star). The light-blue luminous central region represents the inner core, characterized by ultra-dense matter where neutrons and protons break down into quarks and gluons. Surrounding this core lies the outer core (shaded blue-white), composed of a neutron superfluid and proton superconductor, with neutron-superfluid vortices and proton-superconductor flux tubes, respectively [magnified view in the bottom-left inset]. The dark-blue crust has a crystalline lattice structure (not shown) and consists of heavy atomic nuclei and free neutrons and free electrons. The glowing white envelope (often called the outer crust) comprises atomic nuclei and free electrons. The white conical regions show radiation beams emerging from the poles of the pulsar.

§.§ Crust Potential ($V_{\theta}$)

In the absence of the crust potential $V_{\theta}$, Eqs. (<ref>)-(<ref>) govern the interplay between neutron-superfluid vortices and proton-superconductor flux tubes in the outer-core region shown in the schematic diagram of Fig. <ref>. The potential $V_{\theta}$ comes into play in the crust, which lies between the envelope and the outer core of a neutron star and provides pinning sites for vortices and flux tubes (see Fig. <ref> and Ref. <cit.>). At the level of a minimal model for pulsars <cit.> the dynamics of this crust is characterised by a single polar angle $\theta$ <cit.> that evolves as follows:

\begin{equation}
\begin{aligned}
I_c\frac{d^2\theta}{dt^2} &= \frac{1}{N_n}\bigg(\int d^3x\partial_{\theta}V_{\theta}|\psi_n|^2 +\frac{n_n}{n_p}\int d^3x\partial_{\theta}V_{\theta}|\psi_p|^2\bigg)\\
&-\alpha \frac{d\theta}{dt}\,;\\
V_{\theta} ({\bf r}_p) &= V_0\exp \left[-\frac{(|{\bf r}_p|-r_{\rm crust})^2}{(\Delta r_{\rm crust})^2}\right] \tilde{V}(x_{\theta},y _{\theta})\,;
\end{aligned}
\label{eq:crust_pot_eq}
\end{equation}

here, $I_c$ is the moment of inertia of the crust, $n_n/n_p$ is the ratio of the number densities of neutrons and protons, and the slowing down of the crust is controlled by the friction coefficient $\alpha$. We choose $\tilde{V} (x_{\theta},y_{\theta})=3+\cos(n_{\rm crust}x_{\theta})+\cos(n_{\rm crust}y_{\theta})$, with $x_{\theta}=\cos(\theta)x_p+\sin(\theta)y_p$ and $y_{\theta}=-\sin(\theta)x_p+\cos(\theta)y_p$ as in Ref. <cit.>; here, $n_{\rm crust}$ determines the number of pinning sites in the crust, $r_{\rm crust}$ is the radius at which $V_{\theta}$ assumes its maximum value, and $\Delta r_{\rm crust}$ is the thickness of the crust. We use a $2\pi$-periodic version of the coordinates, namely, ${\bf r}_p= (x_p,y_p,z_p)$, which is $\pi$-centered,  with $z_p=\pi$,

\begin{eqnarray}
x_p&=&-\sum_{n=1}^{10}\exp(-\frac{16}{100} n^2)(-1)^n \frac{\sin(n(x-\pi))}{n}\,,\nonumber\\
{\rm{and}}\,\,y_p&=&-\sum_{n=1}^{10}\exp(-\frac{16}{100} n^2)(-1)^n \frac{\sin(n(y-\pi))}{n}\,. \nonumber\\
\label{eq:x_y_periodic}
\end{eqnarray}

For the proton superconductor, in the absence of rotation, we have ${\bf A}_{\rm eff}\to{\bf A}$ and $\phi_{\rm eff}\to \phi$ [Eq. <ref>]. A superconductor that is subjected to rotation (${\bf \Omega} \neq 0$) and which is in a uniform external magnetic field displays a captivating interplay of quantum phenomena. Consider first a non-rotating Type-II superconductor in an external magnetic field; it can display a vortex-lattice phase in which flux tubes are arranged in the form of an Abrikosov lattice <cit.>; the quantized magnetic flux $\Phi_B=\int {\bf A} \cdot d{\bf l} $ passes through each vortex. Each of these magnetised vortices contributes a discrete quantum of magnetic flux to the net magnetic field inside the superconductor, which is zero outside the vortices. Next consider a rotating superconductor without an external magnetic field; this displays a uniform magnetic field, known as the London moment <cit.>. [This field is uniform away from the boundary, for distances larger than $\lambda_p$, the London penetration depth.] The London moment follows from  ${\bf A}_{\rm eff}$ [Eq. <ref>] because the RTGLPE [Eq. <ref>] has an additional term with a vector potential (the subscript $L$ stands for London)

\begin{equation}
{\bf A}_{\rm L} = -\frac{m_p}{q}({\bf \Omega }\times {\bf r}) \,,
\label{eq:london_vecpot}
\end{equation}

so that, in the absence of flux tubes,  ${\bf A}_{\rm eff}=0$ inside the superconductor. For a rotating Type II superconductor <cit.>, in a uniform magnetic field, there is a critical rotational speed ${\rm \Omega}_c^p$ beyond which vortices (here, proton-superconductor flux tubes) enter the system. The critical ${\rm \Omega}_c^p$, which follows by minimizing the energy in the rotating frame $E'\equiv E - {\bf \Omega} \cdot {\bf L}_z$, with ${\bf L}_z$ the angular momentum along the rotation axis, is

\begin{eqnarray}
{\rm \Omega}_c^p = \frac{\hbar }{m_p\lambda_p^2}\ln \bigg(\frac{\lambda_p}{\xi_p}\bigg)\,,
\label{eq:omega_c_super}
\end{eqnarray}

where $\xi_p$, and $\lambda_p$ are, respectively, the superconducting coherence length and the London penetration depth.
In the Abrikosov-lattice phase, a London moment ($\propto {\bf \Omega }$) is present inside the superconducting together with flux-tube lattice.

§ UNITS AND NUMERICAL METHOD

We use the dimensionless forms of Eqs. (<ref>)-(<ref>), which we obtain by using the general reference length $L_{\rm ref}$ and speed $V_{\rm ref}$. The scaled position ${\bf x}$, time $t$, vector potential ${\bf A}$, and scalar potential $\phi$ are

\begin{eqnarray}
{\bf x} &=& L_{\rm ref}{\bf x}'\,, \nonumber\\
t &=& \frac{L_{\rm ref}}{V_{\rm ref}} t'\,, \nonumber\\
{\bf \Omega} &=& \frac{L_{\rm ref}}{V_{\rm ref}} {\bf \Omega}'\,,\nonumber\\
{\bf A} &=&  \frac{H_{\rm c2}L_{\rm ref}}{\kappa} {\bf A}^{'}\nonumber\,, \\
{\rm {and}} \,\, \phi &=& \frac{ L^2_{\rm ref}}{\tau} \frac{H_{\rm c2}}{\kappa}  \phi^{'}\,,
\label{eq:non_dimensional}
\end{eqnarray}

where $H_{\rm c2}$ is the (zero-temperature) upper critical magnetic field of the superconductor, $\kappa = \frac{\lambda}{\xi_p}$ is the London ratio, and $\tau = \frac{L_{\rm ref}}{V_{\rm ref}}$. The wavefunctions are normalized as $\psi_n=\sqrt{n_n}\psi_n'$ and $\psi_p=\sqrt{n_p}\psi_p'$, so the non-dimensionalised neutron GPPE, proton RTGLPE, and vector, gravitational, and scalar potential equations are, respectively, (for notational simplicity we now drop the
primes that come from non-dimensionalization):

\begin{equation}
\begin{aligned}
i\frac{\partial \psi_n}{\partial t} &= -\alpha \nabla^2 \psi_n+ \beta (|\psi_n|^2-1) \psi_n  + \mathfrak{G} \Phi\psi_n  \\
&+i ({\bf \Omega} \times {\bf r}) \cdot \nabla\psi_n+ V_{\theta} \psi_n+ \gamma_p \mathfrak{g} |\psi_p|^2 \psi_n \\
&-\frac{\gamma_p\alpha}{i}  \left[ \nabla \psi_n \cdot {\bf J}_p + \nabla \cdot \left( \psi_n {\bf J}_p\right)\right]\,;
\end{aligned}
\label{eq:GPE_neutron_ndim}
\end{equation}

\begin{equation}
\begin{aligned}
i \left( \frac{\partial }{\partial t} + i\frac{L_{ref}^2}{\kappa \xi_p^2} \phi_{\rm eff}\right) \psi_p &= \alpha \bigg(\frac{\nabla}{i}-\frac{L_{\rm ref}^2}{\xi_p^2\kappa} {\bf A}_{\rm eff}\bigg) \psi_p +\mathfrak{G}\Phi \psi_p\\
&+\beta \frac{\xi_n^2}{\xi_p^2}(|\psi_p|^2 -1)\psi_p +V_{\theta}\psi_p  \\
&+ \gamma_n \mathfrak{g} |\psi_n|^2 \psi_p\\
&-  \gamma_n \alpha (2{\bf J}_n \cdot D_{\bf A}+i  \psi_p \nabla \cdot {\bf J}_n)\,;
\end{aligned}
\label{eq:GLE_proton_ndim}
\end{equation}

\begin{equation}
\begin{aligned}
\frac{V_{ref}^2}{c^2}\frac{\partial^2 {\bf A}}{\partial t^2}-\nabla^2 {\bf A} -\nabla \times {\bf B}_{\rm ext}&=& \mathbb{P}\bigg[\frac{1}{\kappa} {\bf J}_{p}-\frac{\gamma_n}{\kappa} {\bf J}_n|\psi_p|^2\bigg]\,;
\label{eq:vector_pot_ndim}
\end{aligned}
\end{equation}

\begin{eqnarray}
\nabla^2 \Phi &=&|\psi_n|^2 +\frac{n_p}{n_n}|\psi_p|^2 -n_{\rm bg}\,;
\label{eq:gravity_ndim}
\end{eqnarray}

\begin{eqnarray}
\nabla^2\phi &=& -\frac{\beta}{\kappa} \bigg(\frac{c}{c_s}\bigg)^2 (|\psi_p|^2-1)\,.
\label{eq:scalar_pot_ndim}
\end{eqnarray}

The dimensionless current densities and effective vector and scalar potentials are, respectively:

\begin{eqnarray}
{\bf J}_n &=& \frac{1}{2i} (\psi_n^*\nabla \psi_n-\psi_n\nabla\psi_n^*)\,; \nonumber\\
{\bf J}_p &=& \frac{1}{2i} (\psi_p^*\nabla \psi_p-\psi_p\nabla\psi_p^* ) - \frac{L_{ref}^2}{\xi_p^2\kappa} {\bf A}_{\rm eff}|\psi_p|^2\,;\nonumber\\
{\bf A}_{\rm eff}&=&{\bf A}+\frac{\xi_p^2}{L_{\rm ref}^2} \frac{\kappa}{2\alpha} ({\bf \Omega }\times {\bf r})\,;\nonumber\\
{ \phi}_{\rm eff}&=&{\phi}-\frac{\xi_p^2}{L_{\rm ref}^2} \frac{\kappa}{4\alpha} ({ \Omega^2 }{r^2})\,.
\end{eqnarray}

For the convenience of the reader, we define all the parameters and dimensionless ratios in Table (<ref>).

in Eqs.(<ref>)-(<ref>)     Description

$\alpha = \frac{c_s\xi_n}{\sqrt{2}L_{ref}V_{ref}}$      Coefficient of the kinetic term
in  Eqs.(<ref>) and (<ref>)

$\beta = \frac{c_sL_{ref}}{\sqrt{2}\xi_nV_{ref}}$       Coefficient of the nonlinear term
in  Eqs.(<ref>) and (<ref>)

$\thead{\mathfrak{G}=\frac{L_{ref}^3 2\sqrt{2}\pi Gm_n n_n}{V_{ref}c_s\xi_n}}$     Gravitational strength

$\thead{\mathfrak{g}=\frac{L_{ref}g_{np}}{\sqrt{2} \xi_n V_{ref}c_sm_nm_p}}$      Density-Density coupling strength

$\thead{\frac{n_p}{n_n}}$     Number-density ratio
(protons to neutrons)

$\thead{\gamma}$      Dimensional interaction strength

$\thead{\gamma_n=\gamma m_nn_n}$     Interaction coefficient for protons

$\thead{\gamma_p=\gamma m_pn_p}$     Interaction coefficient for neutrons

$\thead{\xi_n}=\frac{\hbar}{\sqrt{2m_n gn_n}}$     Coherence length for
neutron Cooper pairs

$\thead{\xi_p}=\frac{\hbar}{\sqrt{2m_p\alpha_s n_p}}$       Coherence length for
proton Cooper pairs

$\thead{\lambda_p}$     London penetration depth

$\thead{\frac{c}{c_s}}$     Ratio of speed of light to
speed of sound

$\thead{\kappa=\frac{\lambda_p}{\xi_p}}$     London ratio

$\thead{ \Omega}$     Dimensionless rotational speed

Definitions of all the dimensionless parameters and ratios appearing in Eqs. (<ref>)-(<ref>).

We use pseudospectral direct numerical simulations (DNSs) to solve Eq. (<ref>) and Eqs. (<ref>)-(<ref>) in a cubic domain, with side $L=2\pi$ and $N^3$ collocation points, and periodic boundary conditions in all three directions. We employ the Fourier expansion for the function $\Psi\equiv (\psi_n,\psi_p,A_x,A_y,A_z,\Phi,\phi)$ as follows

\begin{eqnarray}
\Psi ({\bf x}) = \sum_{ {\bf k} } \hat{\Psi}_{{\bf k}} \exp (i {\bf k} \cdot {\bf x}) \,,
\label{eq:Fourier}
\end{eqnarray}

and the $2/3$-rule for dealiasing, i.e., we truncate the Fourier modes by setting $\hat{\Psi}_{\bf k} \equiv 0$ for $ |{\bf k}| > k_{max}$ <cit.>, with $k_{max} = [N/3]$. Given current computational resources, it is well-nigh impossible
to use such a DNS with astrophysically realistic values (say for a pulsar) for the parameters and ratios in Table (<ref>).
Nevertheless, as we show below, it is possible to obtain a large body of results that are qualitatively relevant for (a) interactions between proton-superconductor flux tubes and neutron-superfluid vortices and (b) a minimal model for pulsars and their glitches <cit.>.

§ RESULTS

Our results are presented in the following Subsections:

* Section <ref>: we solve the imaginary-time versions of Eqs. (<ref>)-(<ref>) without any interactions ( $\gamma=0$), no crust potential ($V_{\theta}=0$), and $\Theta=0$, where $\Theta$ is the angle between the rotation axis and external magnetic field. Note that, even if there is no direct interaction between the neutron-superfluid and the proton-superconductor (i.e., $\gamma=0$), they interact indirectly through the gravitational Poisson equation (<ref>).

* Section <ref>:
we solve the imaginary-time  versions of Eqs. (<ref>)-(<ref>), but with $\gamma \neq0$,  $V_{\theta}=0$, and $\Theta=0$.

* Section <ref>:
we solve the imaginary-time  versions of Eqs. (<ref>)-(<ref>), but with $V_{\theta}=0$ and $\Theta=30^{\circ}$ and
(i) $\gamma=0$ [Sec. <ref>] and (ii) $\gamma\neq 0$ [Sec. <ref>].

* Section <ref>: we solve the real-time Eqs. (<ref>)-(<ref>) with non-zero crust potential ($V_{\theta}\neq0$), no direct interactions ($\gamma=0$), and $\Theta=0$. Note that the imaginary-time evolution in Secs. <ref>, <ref>, and <ref> has no dynamical significance; this evolution just provides us with a convenient way of obtaining the equilibrium configuration at very large imaginary time.

§.§ Imaginary-time study : $\gamma=0$, $V_{\theta}=0$, $\Theta=0$

We solve imaginary-time ($t\to -i t$) versions of Eqs. (<ref>)-(<ref>) together with Eqs. (<ref>)-(<ref>) with $\gamma=0$, $V_{\theta}=0$, and $\Theta=0$. The imaginary-time version of the Maxwell equation (<ref>) is given in Eq. (<ref>). For all the imaginary-time studies, we start with the initial conditions ${\bf ICI1}$, ${\bf ICI2}$, and ${\bf ICI3}$ given in Appendix <ref>. The neutron star in our model rotates with an angular velocity ${\bf \Omega} =\Omega \hat{z}$; for specificity, we choose $\Omega=2.5$. Both the neutron-superfluid and proton-superconductor subsystems also rotate with this frequency; if we include an external magnetic field ${\bf B}_{\rm ext }$, the proton-superconductor responds directly to it.
However, to isolate the effects of the rotation and the external magnetic field, it is useful to study the following three cases: (i) the neutron-superfluid rotates, but not the proton-superconductor, which is in an external magnetic field; (ii) both the neutron-superfluid and the proton-superconductor rotate, but there is no external magnetic field; and (iii) both the neutron-superfluid and the proton-superconductor rotate, and there is an external magnetic field. Clearly, only case (iii) is directly relevant to neutron stars.

Case(i): The neutron condensate rotates [${\bf \Omega} = \Omega \hat{z}$] and the non-rotating proton condensate is in an external magnetic field ${\bf B}_{\rm ext}=B \hat{z}$. In Figs. <ref>(a)-(b) and Figs. <ref>(c)-(d) we show via contour plots of $|\psi_n|^2$  and  $|\psi_p|^2$, respectively, that the neutron condensate is threaded by vortices and the proton-superconductor displays an Abrikosov flux lattice. Each vortex in this lattice carries the quantum of magnetic flux $\Phi_B = \oint {\bf A} \cdot d{\bf l}$. Initially, the magnetic field is confined near the boundary [Fig. <ref>(e)]. Eventually these vortices penetrate the condensate [Fig. <ref>(f)] and each one of them contributes the unit magnetic flux $\Phi_B = \oint {\bf A} \cdot d{\bf l}$ to the overall magnetic field, which is confined within quantized flux tubes, solely inside the proton superconductor.
Initially we insert flux tubes inside the proton condensate [Fig. <ref>(c) and ICI2 in Appendix <ref>] because, in the absence of rotation, the system cannot generate flux tubes at $T=0$.

(-90,90) (a)
(-90,110) Initial time: $t_i$
(-138,50) $|\psi_n|^2 $
(-116,50)  $\bf \rightarrow$

(-90,90) (b)
(-90,110) final time: $t_f$

(-90,90) (c)
(-138,50) $|\psi_p|^2 $
(-116,50)  $\huge \to$

(-90,90) (d)

(-90,90) (e)
(-138,50) $|{\bf B}|^2 $
(-116,50)  $\huge \to$

(-90,90) (f)
One-level contour plots of (${\bf a}$)-(${\bf b}$) the neutron Cooper-pair density $|\psi_n|^2$ and (${\bf c}$)-(${\bf d}$) proton Cooper-pair density $|\psi_p|^2$, obtained by using the imaginary time ($t\to -it$) versions of the GPPE (<ref>) and RTGLPE  (<ref>). Pseudocolor plots (${\bf e}$)-(${\bf f}$) of the magnetic field ${\bf B}=\nabla \times {\bf A}$ at the mid-plane $z=L/2$.
Here, the neutron condensate rotates with and angular velocity ${\bf \Omega} = \Omega \hat{z}$, with $\Omega=2.5$; the non-rotating proton condensate is placed in an external magnetic field ${\bf B}_{\rm ext}=B \hat{z}$, with $B=0.8$; neutron and proton Cooper pairs interact only through the gravitational potential (i.e., $\gamma=0$).

Case(ii): both the neutron-superfluid and the
proton-superconductor rotate [with ${\bf \Omega} = \Omega \hat{z}$], but there is no external
magnetic field [${\bf B}_{\rm ext}=0$]. The neutron condensate is threaded by vortices [Figs. <ref>(a)-(b)], beyond a critical angular velocity, as in Case (i). The behaviour of the proton superconductor presents a compelling contrast. For slow rotation, less than a critical $\Omega_c^p$ (we use $\Omega=2.5$ here), the proton superconductor assumes a spherical shape devoid of vortices, as illustrated in Fig. <ref>(a). Furthermore, within the superconductor, a uniform London field emerges (Fig. <ref>(c)), first from the superconductor's boundary at the characteristic length scale $\lambda_p$, which is the London penetration depth  [see Section <ref>]. In our simulation with Eq. (<ref>), the magnitude of the dimensionless London field ${\bf B}_L = \nabla \times {\bf A}_L$ is

\begin{equation}
{\rm B}_L=\frac{\xi_p^2}{L_{\rm ref}^2}\frac{\kappa}{\alpha}\Omega \simeq 1.12\ \Omega \,.
\end{equation}

The London magnetic field, a fundamental property of rotating superconductors, has been studied in a hydrodynamical model of a neutron-star interior in Ref. <cit.>. The generation of this uniform magnetic field is
facilitated by the macroscopic London current ${\bf J}_{\rm L}$, which is concentrated near the superconductor's boundary, as we show via 2D vector plots of $J_x$ and $J_y$  in Fig. <ref>(a).

(-95,88) (a)
(-70,100) $\Omega = 2.5 $
(-138,50) $|\psi_p|^2 $
(-116,50)  $\bf \rightarrow$

(-95,88) (b)
(-70,100) $\Omega>\Omega_c=4.5$

(-95,88) (c)
(-138,50) $|{\bf B}|^2 $
(-116,50) $\bf \rightarrow$

(-95,88) (d)

Pseudocolor plots of (${\bf a}$)-(${\bf b}$) the proton Cooper-pair density $|\psi_p|^2$ and (${\bf c}$)-(${\bf d}$) the magnetic field ${\bf B}=\nabla\times {\bf A}$  at the mid plane $z=L/2$. In columns 1 and 2, $\Omega=2.5$ and $\Omega>\Omega_c=4.5$, respectively. In these plots, we have ${\bf B}_{\rm ext}=0$. Neutron and proton Cooper pairs interact only indirectly through the gravitational potential (i.e., $\gamma=0$).

(-220,210) (a)
(-110,210) (b)
(-210,100) (c)
Two-dimensional (2D) vector plots of the proton Cooper-pair current densities $J_x$ and $J_y$ for (${\bf a}$) $\Omega=2.5$ and (${\bf b}$) $\Omega>\Omega_c=4.5$. (c) The magnitude of the radial component of the magnetic field ${\bf B}=\nabla\times {\bf A}$, normalized by the London magnetic field $B_L$, plotted as a function of the distance $r/R$ from the centre of the superconductor, where $R$ is the radius of the spherical proton condensate. The solid curve is from the imaginary-time DNS of Eqs. (<ref>)-(<ref>) and the dashed curve is the analytical relation (<ref>) for two values of $\chi = \frac{R}{\lambda_p}$, with $\lambda_p$ the London penetration depth of the superconductor.

If ${\bf B}_{\rm ext}=0$, the rotation is so slow that there are no flux tubes [$\Omega < {\rm \Omega}_c^p$], and $\gamma=0$, then we can write Eq. (<ref>) in the steady state as

\begin{equation}
\begin{aligned}
\nabla \times {\bf B}= \frac{q}{m_pc^2\epsilon_0}{\bf J}_s  -
\frac{qn_p}{c^2\epsilon_0} ({\bf \Omega \times {\bf r}})\,,
\end{aligned}
\end{equation}

where $n_p=|\psi_p|^2$. We now use the London equation

\begin{equation}
\nabla \times {\bf J}_s = -\frac{n_pq}{m_p}{\bf B}\,,
\label{eq:London}
\end{equation}

with ${\bf J}_s=m_pn_p {\bf v}$, to obtain

\begin{equation}
\begin{aligned}
\nabla \times \nabla \times {\bf v} = -\frac{1}{\lambda_p^2}({\bf v}- {\bf \Omega}\times {\bf r})\,,
\label{eq:internal_vecpot}
\end{aligned}
\end{equation}

where $\lambda_p = \sqrt{\frac{m_pc^2 \epsilon_0}{qn_p^2}}$ is the London penetration depth. [A similar relation has been used in Ref. <cit.> but for a multicomponent superconductor.] If we assume that the density distribution is spherically symmetric, then the proton superconductor has only the azimuthal component ${\bf v}=v{\bf e}_{\phi}$, so, by solving Eq. (<ref>), we get

\begin{equation}
\begin{aligned}
{\bf v} = \bigg[ { \Omega }r+ \frac{C}{r^2} (\sinh\bigg(\frac{r}{\lambda_p}\bigg)-\frac{r}{\lambda_p} \cosh\bigg(\frac{r}{\lambda_p}\bigg))\bigg]{\bf e}_{\phi}\,,
\label{eq:internal_vel}
\end{aligned}
\end{equation}

whence we obtain the radial component ${\rm B}_r$ of the magnetic field by using Eqs. (<ref>) and (<ref>):

\begin{equation}
\begin{aligned}
{\rm B}_r = \frac{m_p}{q}\bigg[ 2{ \Omega }+ \frac{2C}{r^3} \bigg(\sinh(\frac{r}{\lambda_p})-\frac{r}{\lambda_p} \cosh(\frac{r}{\lambda_p})\bigg)\bigg]\,.
\label{eq:internal_mag}
\end{aligned}
\end{equation}

We determine the constant $C$ by demanding ${\rm B}_{r={\rm R}} =0$, where ${\rm R}$ is the radius of the spherical condensate, because ${\bf B}_{\rm ext}=0$. Finally, we get

\begin{eqnarray}
\frac{{\rm B}_r}{{\rm B}_L}=1-\frac{1}{(r/\rm R)^3}\frac{\sinh(\frac{r}{\rm R}\chi)-\frac{r}{\rm R}\chi \cosh(\frac{r}{\rm R}\chi)}{\sinh(\chi)-\chi \cosh(\chi)}\,,
\label{eq:internal_field_analy}
\end{eqnarray}

where $\chi = \frac{\rm R}{\lambda_p}$, which we show via dashed lines in  Fig. <ref>(c) for
two representative values of $\chi$. The solid lines in
Fig. <ref>(c) give the results of our
DNS, which agree well with the results of our analytical approximation [Eq. (<ref>)]. We observe that, for large value of $\chi$ (small $\lambda_p$), the internal magnetic field is comparable to the London field ${\rm B}_L$; and it decreases as we decrease $\chi$.

Upon increasing the angular velocity ${ \Omega}$ beyond $\Omega_c^p \simeq 4.5$ [see Eq. (<ref>)], vortices begin to penetrate the proton superconductor, as we show in  Fig. <ref>(b). Each of these vortices carries a quantum of magnetic flux $\Phi_{\Omega} = \oint ({\bf \Omega }\times {\bf r}) \cdot d{\bf l} $, which is the flux because of the uniform London field ${\rm B}_L$ within the superconductor  (Fig. <ref>(d)).
Our DNS provides valuable insights into the dynamic interplay between rotation, vortices, and the London field in the proton superconductor. As vortices penetrate the condensate <cit.>, they reduce the uniform magnetic field through the system [This is in contrast to the behavior of the Abrikosov lattice, for a non-rotating superconductor in an external magnetic field <cit.>.]. The distribution of currents, crucial for generating and sustaining this magnetic field, is revealed by the 2D  vector plots of $J_x$ and $J_y$ in Fig. <ref>(b). Note that these currents are concentrated near the superconductor's boundary and around the vortices.

Case(iii): Both the neutron and proton condensates rotate with ${\bf \Omega}=\Omega\hat{z}$; and the proton condensate is subjected to an external magnetic field  ${\bf B}_{\rm ext}=B\hat{z}$. The equilibrium state of the neutron condensate resembles that of Case (ii), with vortices penetrating the system beyond the neutron critical angular velocity $\Omega_c^n$. The proton condensate manifests a London field ${\rm B}_L$ as in Case (ii). In Figs. <ref>(a)-(b) we present blue-scale plots of the magnetic field ${\bf B} = \nabla\times {\bf A}$, in the plane $z=L/2$ and at initial $t_i$ and final $t_f$ representative (imaginary) times; these plots show neutron vortices and proton flux tubes in red and black contours respectively.

At the initial (imaginary) time $t_i$ we include proton flux tubes (black) in the initial condition [Fig. <ref>(a)]; the magnetic field is concentrated near the boundary of the simulation box. In the equilibrium state, neutron vortices (in red) enter the condensate and organize themselves in the manner depicted in Fig. <ref>(b). Even though $\gamma = 0$, i.e., there is no direct interaction between the two components, the neutron-superfluid vortices and proton-superconductor flux tubes come close together, as we show in Fig. <ref>(b). This effective attraction follows from the coupling induced by the Poisson equation for gravitational potential [Eq.( <ref>)].

(-55,110)  $ t_i$
(-90,88) (a)

(-55,110)  $ t_f$
(-90,88) (b)
Pseudocolor plots, at the plane $z=L/2$, illustrating the magnetic field ${\bf B} = \nabla\times {\bf A}$ in blue at (a) initial $t_i$ and (b) final $t_f$ (imaginary) times, with neutron vortices and proton flux tubes indicated by red and black contours, respectively. Both neutron and proton
subsystems rotate with $\Omega=2.5$; and the external magnetic field is $B_{\rm ext}=0.8$.

(-110,90) (a)

(-110,90) (b)

(-110,90) (c)

Pseudocolor plots, at plane $z=L/2$, illustrating the magnetic field ${\bf B} = \nabla\times {\bf A}$ in blue [at equilibrium, i.e., the final imaginary time in our DNS], with neutron vortices and proton flux tubes shown via red and black contours, respectively, for (a) gravitational and attractive ($g'<0$) density-density interactions, (b) gravitational and repulsive ($g'>0$) density-density interactions, (c) gravitational, attractive density-density, and also current-current interactions. (d) Imaginary-time series plots of the magnetic-field energy $\rm E_{B}=|{\bf B}_{\rm ext}-\nabla\times {\bf A}|^2$ for zero current-current interaction ($\gamma=0$) in cyan and nonzero current-current interaction ($\gamma\neq0$) in blue. Both neutron and proton subsystems rotate with $\Omega=2.5$; and the external magnetic field is $B_{\rm ext}=0.8$.

§.§ Imaginary-time study: $\gamma\neq 0$, $V_{\theta}=0$, $\Theta=0$

We examine the interacting case of $\gamma\neq 0$
[in addition to the gravitational-potential-induced coupling in Case(iii)].
Consider first neutron-proton density-density interactions; their effective strength follows from the first term in the interaction Lagrangian (<ref>):

\begin{equation}
g'=g_{np} |\psi_n|^2 |\psi_p|^2\,.
\label{eq:gprime}
\end{equation}

For attractive (repulsive) couplings $g' < 0$ ($g' > 0$), the density minima of neutron-superfluid vortices and proton-superconductor flux tubes align (the maximum of the neutron-superfluid density aligns with the proton flux tubes).
[A similar analysis has been conducted in Ref. <cit.> without a gravitational interaction.] Given that the gravitational interaction is inherently attractive, the attractive case ($g' < 0$) only promotes this alignment:
If we start with the same initial condition as in Fig. <ref>(a), we obtain the final equilibrium state, shown in Fig. <ref>(a), with overlapping density minima of the neutron-superfluid vortices in red and of proton-superconductor flux tubes in black (as in Fig. <ref>(b), with only gravitation-induced interactions). In contrast, if $g' > 0$, the interplay between the repulsive density-density interaction and the attractive gravitational interaction
is such that cores do not overlap, as evident in Fig. <ref>(b).

We now incorporate the current-current interaction in the Lagrangian (<ref>). This leads to the entrainment of the proton-superconductor current because of the term $\frac{\gamma q}{c^2\epsilon_0} {\bf J}_n|\psi_p|^2$ in the vector potential (<ref>). This introduces a combination of gravitational, negative density-density ($g'<0$), and positive current-current interactions, whose effects we investigate by starting with the same initial condition as in Fig. <ref>(a). As imaginary time progresses, neutron-superfluid vortices (in red) and proton-superconductor flux tubes (in black) tend to merge and the minima of neutron and proton densities overlap because of the gravitational and negative density-density interaction (see Fig. <ref>(c)).
This entrainment causes neutron-superfluid vortices to drag proton-superconductor flux tubes and induces a magnetic field inside these vortices. This is clearly visible as intense blue spots in Fig. <ref>(c), inside vortices that do not overlap with proton-superconductor flux tubes. To illustrate the generation of this entrainment-induced magnetic field, we plot the magnetic-field energy $\rm E_B(t) = |{\bf B}_{\rm ext}-\nabla \times {\bf A}|^2$ versus imaginary time $t$ in Fig. <ref>(d) both with (blue curve) and without (cyan curve) current-current interaction. In the latter case $\rm E_B(t)$ is higher than in the former at large $t$.

§.§ Angle between ${\bf B}_{\rm ext}$ and ${\bf \Omega}$

So far we have examined cases with aligned  ${\bf B}_{\rm ext}$ and ${\bf \Omega}$, i.e., $\Theta=0$.
We turn now to $\Theta > 0$, which is the case in most pulsars.

§.§.§ No interactions: $\gamma=0$, $V_{\theta}=0$, and $\Theta=30^o$

As in Case(i) Subsection <ref> we
study a neutron condensate that rotates [${\bf \Omega} = \Omega \hat{z}$] and a non-rotating proton condensate in an external magnetic field ${\bf B}_{\rm ext}=B \hat{z}$.
We examine the case with no direct interactions, i.e., $\gamma=0$, no crust potential, i.e., $V_{\theta}=0$, and the representative values $\Theta=30^o$ and  $\Omega=2.5 > \Omega_c^n$, so vortices enter the neutron condensate, and the proton condensate is stabilized with an Abrikosov lattice.

In Fig. <ref>(a), we show
a one-level red contour plot of the neutron-superfluid vortices at the final imaginary time; these are aligned along the $z$-axis. By contrast, the proton-superconductor flux tubes, illustrated by cyan contour plots in Fig. <ref>(b), form an Abrikosov lattice, have their axes tilted at a fixed angle of $\Theta=30^\circ$ relative to the $z$-axis; the magnetic field manifests itself solely within these flux tubes [Fig. <ref>(c)].

(-70,60) (a)

(-70,60) (b)

(-70,60) (c)

One-level contour plots of $(\nabla \times (\rho v))^2$ for (${\bf a}$) neutron-superfluid vortices, and (${\bf b}$)  proton-superconductors flux tubes at the final imaginary time. (${\bf c}$) Volume plot (final imaginary time) of the magnetic field ${\bf B} = \nabla\times {\bf A}$. In these plots, the neutron condensate rotates with angular velocity ${\bf \Omega} = \Omega \hat{z}$, with $\Omega=4.0$, and a non-rotating proton condensate; there is an external magnetic field $B_{\rm ext}=0.8$ that makes an angle $\Theta=30^\circ$ with the $z$-axis. Both species interact only through the gravitational potential.

We now consider the counterpart of Case(ii) Subsection <ref>: both neutron and proton condensates rotate with an angular velocity ${\bf \Omega}=\Omega\hat{\bf z}$ and there is no external magnetic field [${\bf B}_{\rm ext}=0$].  Vortices, oriented along the $z$-axis, thread the neutron condensate. At small values of the imaginary time $t$, the proton-superconductor flux tubes are oriented at an angle $\Theta=30^o$ [Fig. <ref>(a)]. As $t$ increases, these flux tubes try to align themselves along the rotation axis $z$ [ Fig. <ref>(b)], and they do so ultimately [Fig. <ref>(c)]. This alignment is facilitated by the absence of an external magnetic field, which tends to counteract this reorientation. We also observe the generation of a London magnetic field ${\rm B}_L$ inside the proton condensate [Fig. <ref>(d)-(f)].

(-70,60) (a)

(-70,60) (b)

(-70,60) (c)

(-70,60) (d)

(-70,60) (e)

(-70,60) (f)

(-70,60) (g)

(-70,60) (h)

(-70,60) (i)

(a)-(c) One-level contour plots of $(\nabla \times (\rho v))^2$ for proton flux tubes and (d)-(f) pseudocolor plots of the magnetic field B at the midplane ($z=L/2$) at three representative imaginary times. Both neutron and proton subsystems rotate with the angular velocity ${\bf \Omega} = 4.0 \hat{z}$; and ${\bf B}_{\rm ext}=0$. (g)-(i) One-level contour plots of $(\nabla \times (\rho v))^2$ for proton flux tubes with ${\bf \Omega} = 4.0 \hat{z}$; and ${\bf B}_{\rm ext}=0.8$.  At the initial imaginary time in (a), the proton flux tubes make an angle $\Theta=30^{\circ}$ with the $z$-axis. Both species interact only through the gravitational potential.

Finally we consider the counterpart of Case(iii) Subsection <ref>: Both neutron and proton condensates rotate [${\bf \Omega}=\Omega\hat{\bf z}$] and there is an external uniform magnetic field ${\bf B}_{\rm ext}$, as in
a pulsar, but with interactions solely through the gravitational potential. At small values of the imaginary time $t$ [Fig. <ref>(g)], the proton-supercondcutor flux tubes are aligned at an angle $\Theta=30^o$ with the $z$-axis, and the magnetic field is concentrated primarily outside the condensate. As $t$ increases, the proton flux tubes try to orient themselves along the rotation axis [Fig. <ref>(h)], but
${\bf B}_{\rm ext}$ resists this alignment. Ultimately, these proton-superconductor flux tubes exhibit frustration [Fig. <ref>(i)] as they try both to align globally with the rotation axis and to adhere to ${\bf B}_{\rm ext}$, which makes an angle $\Theta=30^\circ$ with the z-axis. This frustration is akin to the glassy behavior of flux tubes, studied in Ref. <cit.>, without gravity but with a quadratic confining potential.
In the next Subsection we go beyond the study of Ref. <cit.> by incorporating the full Maxwell equations that lead to entrainment.

§.§.§ Non-zero interactions: $\gamma \neq 0$, $V_{\theta}= 0$, and $\Theta=30^o$

Consider now the counterpart of Subsection <ref>, i.e.,
direct interaction ($\gamma\neq0$) between neutron and proton condensates. We first examine attractive density-density interactions [Eq. (<ref>)], so $g' < 0$.
We show the evolution (in imaginary-time $t$) of the proton-superconductor flux tubes and the magnetic field in Figs. <ref>(a)- (c). At the initial time  $t=0$, the proton flux tubes are oriented at an angle $\Theta=30^o$ with respect to the $z$-axis [Fig. <ref>(a)]. As $t$ increases, these flux tubes exhibit a rapid realignment with the rotation axis because of the combined gravitational and attractive density-density interactions [Figs. <ref>(b)-(c)].

(-70,65) (a)

(-70,65) (b)

(-70,65) (c)

One-level contour plots of $(\nabla \times (\rho v))^2$ for proton flux tubes at three representative imaginary times in (a), (b), and (c). Both neutron and proton subsystems rotate with the angular velocity ${\bf \Omega} = \Omega \hat{z}$, where $\Omega=4.0$; and ${B}_{\rm ext}=0.8$, which makes an angle $\Theta=30^{\circ}$ with the $z$-axis. Furthermore, we have gravitational and attractive density-density ($g'<0$) interactions (first term in Eq. (<ref>)) between neutron and proton Cooper pairs.

To obtain entrainment, we must introduce the current-current interactions [the second term in Eq. (<ref>)]. The proton flux tubes evolve in imaginary time $t$ as in the previous case with attractive density-density interactions. However, the current-current interaction, $\frac{\gamma q}{c^2\epsilon_0} {\bf J}_n|\psi_p|^2$ in the vector potential (<ref>), leads to an entrained proton-superconductor current. This entrainment results in an induced magnetic field inside the neutron-superfluid vortices [gray and red isosurfaces, respectively, in Fig.<ref>].

One-level contour plots of $(\nabla \times (\rho v))^2$ for neutron-superfluid vortices, at the final imaginary time with a superimposed volume plot of the magnetic field ${\bf B}$ at the final time. Both neutron and proton subsystems rotate with the angular velocity ${\bf \Omega} = \Omega \hat{z}$, where $\Omega=4.0$; and ${B}_{\rm ext}=0.8$, which makes an angle $\Theta=30^{\circ}$ with the $z$-axis.  Gravitational and other interaction terms are included  [Eq. (<ref>))].

§.§ Real-time evolution

We now delve into the real-time equations (<ref>)-(<ref>).
We first follow the dynamics of the alignment of proton-superconductor flux tubes with the rotation axis [Fig. <ref>]. We begin with the configuration of vortices (cyan) and flux tubes (red)
shown in Fig.<ref>(a), which we obtain as the equilibrium state of the imaginary-time versions of Eqs.(<ref>)-(<ref>). The axes of rotation and the magnetic moment ${\bf m}$ make an angle $\chi$ [see the
schematic diagram in Fig.<ref>(b)]; we define them as follows:

\begin{eqnarray}
{\bf m}&=&\frac{1}{2}\int {\bf r} \times {\bf J}_p dV \,;\\
\cos(\chi) &=& \frac{{\bf \Omega}\cdot {\bf m}}{|{\bf \Omega}||{\bf m}|}\,.
\label{eq:mag_mom}
\end{eqnarray}

The angle $\chi$ depends on time $t$; for considerable lengths of time it shows minor fluctuations, but, occasionally, it changes dramatically, as we show in Fig.<ref>(c) via a plot of $\cos(\chi)$ versus $t$. The sudden change in $\cos(\chi)$, from positive to negative values, indicates that the magnetic moment undergoes reversals, which are familiar
in many other hydrodynamical <cit.>, dynamo <cit.>, geomagnetic <cit.>, and astrophysical <cit.> settings.
Measurements for several pulsars <cit.> indicate that $0^o\leq \chi \leq 90^o$; e.g., the two pulsars (PSR B1055-52 and PSR 1822-09) are aligned rotators, whereas another two (PSR B1702-19) are orthogonal rotators; the time scales of these observations are such that, for any given pulsar, $\chi$ has a steady, time-independent value.

In our minimal model the crust potential $V_{\theta}$ is a function of a single polar angle $\theta$ whose dynamics is given by Eq. (<ref>). To study the interplay of crust, neutron-superfluid, and proton-superconductors, we use Eqs. (<ref>-(<ref>) along with Eq. (<ref>) and initial conditions
that we obtain from the equilibrium states of imaginary-time studies of the previous Sections, with the same angular velocity $\Omega=\frac{d\theta}{dt}$ for the crust, the neutron condensate, and the proton condensate. For purposes of illustration, we consider $\Theta=0$, i.e., the angle between the rotation axis and the magnetic field is zero. The neutrons and protons interact only through the gravitational potential via the Poisson equation ($\gamma = 0$).

Real-time evolution: (a) One-level contour plot of $(\nabla \times (\rho v))^2$ for neutron vortices (in red) and proton flux tubes (in cyan) at the initial time. (b) Schematic diagram showing the angle $\chi$ between the rotation axis and the magnetic moment [Eq. (<ref>)]. (c) The evolution of the angle $\chi$ with time. Both neutron and proton subsystems rotate with angular velocity ${\bf \Omega} = \Omega \hat{z}$, where $\Omega=4.0$; and ${B}_{\rm ext}=0.8$, which makes an angle $\Theta=30^{\circ}$ with the $z$-axis. Insets (d) and (e) show illustrative proton flux-tube configurations before and after the reversal.

In Figs. <ref>(a)-(c), we show isosurface plots of the crust potential  in blue, isosurfaces of neutron-superfluid vortices in red, and proton-superconductor flux tubes in cyan at three representative times. At the initial time step (Fig.<ref>(a)), the system features $12$ neutron vortices and $6$ proton flux tubes. Given the parameters of our simulation, the proton-superconductor flux tubes are more effectively pinned by the crust potential than
neutron-superfluid vortices, in part because of the anchoring of the flux tubes to the strong external magnetic field ${\bf B}_{\rm ext}$. Over the duration of the simulation presented in Figs. <ref>(a)-(c), the number of neutron-superfluid vortices reduces by a factor of $2$ but the number of proton flux tubes remains unchanged. Both superfluid vortices and proton-superconductor flux tubes undergo differential rotation because of the friction coefficient $\alpha$ in equation (<ref>)) for the polar angle.

The angular momentum of the crust is given as ${\rm J}_c = I_c d\theta/dt$, with $I_c$ is the moment of inertia of the crust [see the Appendix <ref> for details]. The temporal evolution of ${\rm J}_c$ is complicated because of the subtle interplay between the friction, which slows down the crust, and the angular momentum in the neutron-superfluid vortices. When such a vortex is ejected from the collapsed pulsar, its angular momentum is transferred to the superfluid and thence to the crust.
Some neutron vortices linger close to proton-superconductor flux tubes because of the  Poisson-equation-induced gravitational attraction between them. This also affects the time-dependence of ${\rm J}_c$. Finally, we have an effective stick-slip dynamics for ${\rm J}_c$ that displays glitches whose statistics has properties that are akin to those seen in several pulsars and which yield pulsar glitches <cit.>.

One-level contour plots of the crust potential together with the neutron vortices (in red) and proton flux tubes (in cyan) at three different times in (a), (b), and (c) obtained by using the real-time GPPE (Eq. (<ref>)) and  RTGLE (Eq. (<ref>)). Both neutron and proton subsystems rotate with an angular velocity ${\bf \Omega} = \Omega \hat{z}$, where $\Omega=4.0$; and ${B}_{\rm ext}=4.0$, which is along the $z$-axis.

We now examine the analogues of pulsar glitches in our model, by following the methods developed in Ref. <cit.>. In particular, we present the time series of the angular momentum $ ({\rm J}_c-{\rm J}_{c_0})/{\rm J}_{c_0}$ of the crust in Fig.<ref>(a). This time series of $ ({\rm J}_c-{\rm J}_{c_0})/{\rm J}_{c_0}$ exhibits characteristic features that are associated with Self-Organized Criticality (SOC) <cit.>, which we have explored, in the context of gravitationally collapsed boson stars, in our earlier work <cit.>. Figures <ref>(b)-(d) present expanded views of specific segments [indicated by black boxes] of the time series in Fig.<ref>(a). From the time dependence of ${\rm J}_c$, we observe that the crust can either lose angular momentum to the superfluid or can gain angular momentum from it, because of the stick-slip dynamics mentioned above.

To characterize SOC, we quantify the statistics ${\rm J}_c$ as follows: We measure (a) the event size $\Delta  {\rm J}_c$, which is the difference between successive minima and maxima in ${\rm J}_c$, (b) the event-duration time $t_{ed}$, which is the time difference between successive minima and maxima of $ {\rm J}_c(t)$, and (c) the waiting time $t_{w}$, which is the time between successive maxima in $ {\rm J}_c(t)$. We then obtain cumulative probability distribution functions (CPDFs) of $\Delta  {\rm J}_c$, $t_{ed}$,
and $t_{w}$; as in Ref. <cit.>, the former two CPDFs exhibit power-law tails, whereas the last has an exponential tail. In Fig. <ref>(e), we plot the CPDF $Q(\Delta {\rm J}_c/{\rm J}_{c_0})$; it scales as $Q(\Delta {\rm J}_c/{\rm J}_{c_0})\sim (\Delta  {\rm J}_c/{\rm J}_{c_0})^{\beta}$, within the gray-shaded region. Therefore, the corresponding probability distribution function (PDF) scales as $P(\Delta {\rm J}_c/{\rm J}_{c_0})\sim (\Delta {\rm J}_c/{\rm J}_{c_0})^{\beta-1}$; for our representative run, we obtain the scaling exponent $\beta = 0.86\pm 0.15$, by using local-slope analysis. The CPDF of $t_{ed}$ shows the power law $Q(t_{ed}\Omega)\sim (t_{ed}\Omega)^{\gamma_t}$ in the gray-shaded region of Fig. <ref>(f), with an exponent $\gamma_t=2.5\pm 0.2$ for our run. The CPDF of $t_w$ shows the exponential form $Q(t_{w}\Omega)\sim \exp(-6.5 t_{w}\Omega)$ [Fig. <ref>(f)]. The qualitative forms of these CPDFs is similar to those seen in experiments, as has been noted in Ref. <cit.>, which uses the GPPE system without the proton superconductor and the Maxwell equations that we include. The values of the exponents $\beta$ and $\gamma$ lie close to those that have been observed for certain pulsars [e.g., PSR J 1825-0935 has glitch-size-PDF exponent $\simeq 0.36 $ ] <cit.>. In Ref. <cit.>, it has been noted that range of glitching sizes depends on $\Omega$. In addition, we find that these sizes also depend on $B_{\rm ext}$.

We have noted above that Poisson-equation-induced interaction between the neutron superfluid and the proton superconductor makes superfluid vortices approach superconductor flux tubes. As the crust decelerates, the neutron vortices leave the condensate abruptly. The associated jumps in $ {\rm J}_c(t)$ are somewhat sharper in time but smaller in amplitude than those in the GPPE model of Ref. <cit.>. Consequently, our values of the exponents $\beta$ and $\gamma$ are about $10\%$ larger than those in Ref. <cit.>, but still comfortably in the experimental range <cit.> $-0.13\lesssim -(\beta-1)\lesssim 2.4$. The inclusion of current-current and density-density interactions [Eqs. (<ref>)] may reduce the sizes of glitches, by slowing down the ejection of vortices from the condensate. Furthermore, we expect that the current-entrainment term in Eq. (<ref>), which induces a magnetic field inside neutron vortices, could reduce glitch sizes.

(-450,88) (a)
(-340,88) (b)
(-220,88) (c)
(-100,88) (d)

(-410,110) (e)
(-250,110) (f)
(-90,110) (g)
(a) Time series of the crust angular momentum $\rm (J_c-J_{c_0})/J_{c_0}$. (b), (c), and (d) are the zoomed versions of the rectangular regions shown in the preceding plots. Log-Log plots of ${\bf (e)}$ the CPDF $\rm Q(\Delta J_c/J_{c_0})$ of the event size and ${\bf (f)}$ the CPDF $\rm Q(\rm t_{ ed}\Omega)$ of the event duration. ${\bf (g)}$ semilog plot of the CPDF $\rm Q(\rm t_{ w}\Omega)$ of the waiting time. $\rm J_{c_0}$ and $\rm \Omega$ are the initial angular momentum and initial angular velocity of the crust, respectively.

§ CONCLUSIONS

We have developed a theoretical framework for studying the coupled motion of neutron-superfluid vortices and proton-flux tubes in a gravitationally collapsed condensate. In this framework we have employed (a) a 3D Gross-Pitaevskii-Poisson-Equation (GPPE) for neutron Cooper pairs, (b) the Real-Time-Ginzburg-Landau equation (RTGLE) for proton Cooper pairs, (c) the Maxwell equation for the vector potential ${\bf A}$, and (d) Newtonian gravity and interactions, both direct and induced by the Poisson equation, between the neutron and proton subsystems. For a pulsar we have included, in addition, a crust potential as in Ref. <cit.>. The recent studies in Refs. <cit.> use the Gross-Pitaevskii-Equation and Ginzburg-Landau equation (imaginary time) together with a static Ansatz for ${\bf A}$ in a harmonic trap. We have gone well beyond these earlier studies by including Newtonian gravity in the GPPE and RTGLE together with  the complete Maxwell equations for ${\bf A}$. To the best of our knowledge, this has not been attempted hitherto in the context of pulsars.

Our imaginary-time studies of the GPPE (<ref>) and the RTGLE (<ref>) reveals that, even in the absence of any direct interaction ($\gamma=0$), the neutron-superfluid vortices and proton-superconductor flux tubes interact gravitationally through the Poisson equation (<ref>).
By including the Maxwell equations, we demonstrate, for the first time, that neutron-superfluid vortices display an induced magnetization whose magnitude is proportional to $\frac{\gamma q}{c^2\epsilon_0} {\bf J}_n|\psi_p|^2$. This magnetization plays a crucial role in the expulsion of vortices from the pulsar.
The angle $\Theta$ is an important control parameter in our model.
For example, if $\Theta = 30^{\circ}$, proton flux tubes gradually endeavor to align themselves with the rotation axis over time [Figs.<ref>(a)-(c)], but they also exhibit a tendency to adhere to the external magnetic field;
this competition leads to frustration in the proton-superconductor flux tubes, which are no longer straight but become distorted [Figs. <ref>(b) and (h)].

The real-time dynamics of the GPPE (<ref>), RTGLE (<ref>), and the Maxwell equations (<ref>) can be applied qualitatively to pulsars. We must, of course, incorporate a pulsar-crust potential $V_{\theta}$, described, at the simplest level, by the polar angle $\theta$ in Eq. (<ref>). This provides a minimal model for studying pulsar glitches. Our investigation reveals that the proton-superconductor flux tubes remain anchored to the crust by the external magnetic field ${\bf B}_{\rm ext}$, while neutron-superfluid vortices leave the condensate and give rise to the glitching phenomenon, the complicated time evolution of the crust angular momentum ${\rm J}_c(t)$, which displays signatures of self-organised criticality (SOC). Although pulsar glitches have been obtained recently in the GPPE model <cit.>, they have not been studied in the presence of proton-superconductor flux tubes, whose dynamics affects pulsar glitches and the time series of ${\rm J}_c(t)$ significantly [as we can see by comparing Fig. <ref> with Fig. 4(b) in Ref. <cit.>].

The SOC that we obtain in our model for pulsars, which generalises the earlier work from our group <cit.>, is akin to what has been obtained in some pulsars; e.g., in the pulsar PSR J 1825-0935, the glitch-size exponent $\beta \simeq 0.36$. Given the simplicity of our model, this is indeed gratifying. In our simulations, some ratios match closely the values found in typical pulsars: In the outer core of a neutron star, the proton to neutron number density ratio is $\frac{n_p}{n_n}\simeq0.05$; we consider the range $0.5 \leq \frac{n_p}{n_n} \leq 1$. Furthermore, we choose the ratio of the neutron and proton coherence lengths $\frac{\xi_n}{\xi_p}=2$, which agrees with the value found in neutron stars <cit.>. The London parameter $\kappa$ for type II superconducting proton Copper pairs is chosen to be greater than $\frac{1}{\sqrt{2}}$, so that we have an Abrikosov phase. However, it is important to note that numerical studies cannot achieve spatial and temporal scales and resolutions that are in the ranges of direct relevance to pulsars. In particular, the radius of a typical neutron star is $\simeq 10$ km; by contrast, the core sizes of neutron-superfluid vortices are $\simeq 10^{-15}$m; the ratio of the speed of light to that of sound $\frac{c}{c_s}\simeq 10^6$, which is a challenge for any simulation. The number of neutron-superfluid vortices that thread a pulsar is estimated to be $N_v = 10^{16}$; this is far in excess of what can be simulated on even the world's biggest computers. The number of vortices in our model is given as $N_v = \frac{L}{2}\sqrt{\frac{\beta}{\alpha}}$, where $L=2\pi$ is the length of the simulation box and $\alpha$ and $\beta$ are given in Table <ref>. For the values we use in our simulations, namely, $\alpha=0.2$ and $\beta = 5$, we have $N_v\sim 15-20$ vortices. We can increase $N_v$ by increasing either the size of our simulation or the ratio $\beta/\alpha$ (or both); but these are limited severely by computational facilities.

§ ACKNOWLEDGMENTS

We thank the Indo-French Centre for Applied Mathematics (IFCAM), the Science and Engineering Research Board (SERB), and the National Supercomputing Mission (NSM), India for support, and the Supercomputer Education and Research Centre (IISc) for computational resources.

The imaginary-time version of the Maxwell equation (<ref>) is the following first-order partial differential equation [similar to the
vector-potential equation used in the formulation of the time dependent Ginzburg-Landau model of superconductivity <cit.>]:

\begin{eqnarray}
\frac{V_{\rm ref}^2}{c^2}\frac{\partial {\bf A}}{\partial t} - \nabla^2 {\bf A} - \nabla \times {\bf B}_{\rm ext} = \mathbb{P}\bigg[\frac{1}{\kappa}{\bf J}_p-\frac{\gamma_n}{\kappa}{\bf J}_n |\psi_p|^2 \bigg]\,.\nonumber\\
\label{eq:imaginary_time_vec_pot}
\end{eqnarray}

To solve imaginary-time ($t \to -i t$) versions of the GPPE (<ref>), the RTGLE (<ref>), and the Maxwell equation (<ref>), we use the following initial conditions:

* ICI1: The imaginary-time ($t \to -i t$) version of GPPE with ${\bf \Omega}=0$ is first evolved by using a uniform density distribution and with small superimposed perturbations. This gives a spherically collapsed condensate. We now use this collapsed state as an initial condition in the same equation but with a small value of ${\bf \Omega}$. We keep increasing ${\bf \Omega}$ in small steps until we get a collapsed object threaded by vortices.

* ICI2: For the RTGLE, we follow the procedure used in ICI1, but we insert vortices initially by choosing

\begin{eqnarray}
\psi_{pi} = \psi_{\rm uni}\times [\cos(kx)+i \cos(ky)]^n\,,
\end{eqnarray}

where $\psi_{\rm uni}$ is a uniform density distribution with small superimposed perturbations. Here, the integer $n$ denotes the multiplicity of a vortex; and $k$ is the number of vortices in the interval [$-\pi/n,\pi/n$].

* ICI3: For the imaginary-time version of the Maxwell equations (<ref>), we use the following initial condition:

\begin{eqnarray}
A_x &=& -\frac{1}{2}\times y_p\times B_{\rm ext}\,,\nonumber\\
A_y &=& \frac{1}{2}\times x_p\times B_{\rm ext}\,,
\end{eqnarray}

where $x_p$ and $y_p$ are the periodic versions of the coordinates [Eqs. (<ref>)] and $B_{\rm ext}$ is the uniform external magnetic field in the $z$-direction.

The total angular momentum of the system is ${\rm J} = {\rm J}_z+{\rm J}_c$, with ${\rm J}_c$, the crust angular momentum, and ${\rm J}_z$, the angular momentum of the system without the crust, which are, respectively,

\begin{eqnarray}
{\rm J}_c &=& I_c \frac{d\theta}{dt}\qquad {\rm{and}}\nonumber\\
{\rm J}_z &=& \int d^3  {\bf x} \ \psi_n^* (\hat{\bf e}_z\times {\bf r})\cdot (-i\hbar \nabla)\psi_n \nonumber\\
&+& \int d^3 {\bf x}\ \psi_p^* (\hat{\bf e}_z\times {\bf
r})\cdot (-i\hbar \nabla)\psi_p\,,
\end{eqnarray}

where $I_c$ is the moment of inertia of the crust. In the absence of friction [$\alpha=0$ in Eq. ((<ref>))], the total angular momentum is conserved in an infinite system. In the spatially periodic cubical domain that we consider, this conservation is only approximate because this domain does not have strict rotational invariance [see Ref. <cit.> for a detailed discussion].


# New MC determination of the critical coupling in $\phi^{4}_{2}$ theory.

Simone Bronzin1<EMAIL_ADDRESS>Barbara De Palma2
<EMAIL_ADDRESS>Marco Guagnelli1,2
<EMAIL_ADDRESS>1Via Antonini 20, Milano, Italy
2 Università degli Studi di Pavia, Via A. Bassi 6, 27100, Pavia, Italy
3 INFN, Sezione di Pavia, Via A. Bassi 6, 27100, Pavia, Italy

###### Abstract

We investigate the non–perturbative features of $\phi^{4}$ theory in two
dimensions, using Monte Carlo lattice methods. In particular we determine the
ratio $f_{0}\equiv g/\mu^{2}$, where $g$ is the unrenormalised coupling, in
the infinite volume and continuum limit. Our final result is
$f_{0}=11.055(14)$.

Suggested keywords

###### pacs:

12.38.Gc, 11.15Ha

††preprint: APS/123-QED

## Introduction

$\phi^{4}$ theory plays an important role in Quantum Field Theory as it
represents for example an extremely simplified model for the Higgs sector of
the Standard Model.

In $D=2$ dimensions the theory is super–renormalizable: the coupling constant
$g$ has positive mass dimensions $[g]=[\mu_{0}^{2}]$, where $\mu_{0}$ is the
(bare) mass parameter of the theory; this means that the ratio $f\equiv
g/\mu^{2}$, where $\mu^{2}$ is a renormalised squared mass in some given
renormalisation scheme, is the only physically relevant dimensionless
parameter we have to consider. Thanks to the super–renormalisability of the
theory we can use the unrenormalised coupling constant in the definition of
$f$, since in any case the renormalisation of $g$ amounts to a finite
constant.

In this paper we determine the value of $f\equiv g/\mu^{2}$ at the critical
point, that is the value of $f$ computed in the limit in which both $g$ and
$\mu^{2}$ go to zero. We follow the renormalisation scheme used in Loinaz and
Willey (1998); Schaich and Loinaz (2009), adopting the simulation technique
introduced in Korzec  _et al._ (2011), namely the _worm algorithm_ and we
compute the ratio $g/\mu^{2}$ using the same strategy implemented in Bosetti
_et al._ (2015); we present an improvement in the determination of the
critical value $f_{0}$, obtained thanks to the _gradient flow_ Monahan and
Orginos (2015), a technique that allows us to reach smaller values of the
coupling $g$ with respect to our previous work Bosetti  _et al._ (2015).

In the following, after briefly describing the model and the renormalisation
scheme chosen in order to extract $\mu^{2}$ at fixed $g$ in the infinite
volume limit, we will recall the main steps of the simulations, focusing on
the application of the gradient flow; we will then proceed to the continuum
limit extrapolation. In the end we will compare our results with our previous
determination of the same quantity and we will draw some conclusions.

##  I Lattice formulation

Let’s introduce the $\phi^{4}$ Lagrangian in the Euclidean space:

$\mathcal{L}_{E}=\dfrac{1}{2}\left(\partial_{\nu}\phi\right)^{2}+\dfrac{1}{2}\mu_{0}^{2}\phi^{2}+\dfrac{g}{4}\phi^{4}.$
(1)

In $D=2$ the Euclidean action is

$\mathcal{S}_{E}=\int d^{2}x\,\mathcal{L}_{E}.$ (2)

In order to obtain a dimensionless discretized action we put the system on a
2-dimensional lattice with spacing $a$ and linear size $L\equiv Na$. By
introducing the following parametrization

$\hat{\mu}^{2}_{0}=a^{2}\mu^{2}_{0},\qquad\hat{g}=a^{2}g.$ (3)

we have

$\mathcal{S}_{E}=\sum_{x}\left\\{-\sum_{\nu}\phi_{x}\phi_{x+\hat{\nu}}+\dfrac{1}{2}\left(\hat{\mu}_{0}^{2}+4\right)\phi_{x}^{2}+\dfrac{\hat{g}}{4}\phi_{x}^{4}\right\\},$
(4)

where $\phi_{x\pm\hat{\nu}}$ are fields at neighbor sites in the $\pm\nu$
directions.

In the following we will omit the “hat” on top of lattice parameters: all
quantities will be expressed in lattice units, i.e. they become dimensionful
when multiplied by appropriate powers of the lattice spacing $a$.

Numerical simulations are necessarily performed at finite values of $a$, $L$
and of the bare parameters of the theory: to obtain $f_{0}$ we have to
extrapolate our results to $L\to\infty$ (infinite volume or thermodynamic
limit) and to $a\to 0$ (continuum limit).

In order to approach the thermodynamic limit we fix $g$ to a given value in
lattice units, and this amounts to keep fixed the lattice spacing $a$; we then
simulate the system at several values of $N=L/a$. For each value of $N$ we
perform several simulations searching for a value of $\mu^{2}_{0}$ such that a
certain condition is satisfied. This condition, which we will describe in
details later on, is conceived in such a way that by extrapolating
$\mu^{2}_{0}(g,L/a)$ to the infinite volume limit we get a second order phase
transition point in the plane $(g,\mu^{2}_{0})$.

As we discuss in Bosetti  _et al._ (2015), in order to safely go to the
continuum limit, we have to work out a renormalisation of the mass parameter,
since $\mu_{0}^{2}$ in this limit diverges like $\log(a)$. Adhering to the
same renormalization procedure adopted in Loinaz and Willey (1998); Schaich
and Loinaz (2009), we determine the renormalised squared mass $\mu^{2}$
putting it equal to the solution, in the infinite volume limit, of the
equation

$\mu^{2}=\mu^{2}_{0}+3gA(\mu^{2}).$ (5)

where $A(\mu^{2})$ is the only 1–Particle–Irreducible divergent diagram in
$D=2$ on a $N\times N$ lattice:

$A(\mu^{2}_{0})=\dfrac{1}{(2\pi)^{2}}\int\dfrac{\textrm{d}^{2}p}{p^{2}+\mu^{2}_{0}}$
(6)

The condition (5) is equivalent to the introduction of a proper divergent
mass–squared counterterm in the action. We may finally extrapolate the
quantity $f\equiv g/\mu^{2}$ to $g\to 0$ in order to obtain $f_{0}$, the
critical value in the continuum limit.

Another useful parametrization of the action is the following:

$\begin{split}\mathcal{S}_{E}&=-\beta\sum_{x}\sum_{\nu}\varphi_{x}\varphi_{x+\hat{\nu}}+\sum_{x}\left[\varphi^{2}_{x}+\lambda(\varphi^{2}_{x}-1)^{2}\right]\\\
&=\mathcal{S}_{I}+\mathcal{S}_{Site},\end{split}$ (7)

where the relations between $(\mu_{0}^{2},\;g)$ and $(\beta,\;\lambda)$ are:

$\phi_{x}=\sqrt{\beta}\varphi,\qquad\mu_{0}^{2}=2\dfrac{1-2\lambda}{\beta}-4,\qquad
g=\dfrac{4\lambda}{\beta^{2}}.$ (8)

In eq.(7) there is an interaction term between neighbor sites,
$\mathcal{S}_{I}$, with a coupling constant of strength $\beta$ and a term
related to a single site, $\mathcal{S}_{Site}$.

###  I.1 Simulations

Now we outline the general computational strategy, focusing in particular on
the improvements with respect our previous work Bosetti  _et al._ (2015).

In our simulations we used the _worm algorithm_ Korzec  _et al._ (2011) and
used the lattice action given by (7). Operatively we fix a value of $\lambda$
and $L/a$ and search for a value of $\beta$ such that the physical condition

$mL=L/\xi=\textrm{const}=z.$ (9)

is matched for a given and fixed value of $z$. Condition (9) implies that the
correlation length $\xi$ of the system grows linearly with $L$: thus, when
$a/L\to 0$, we arrive at the critical point, where the correlation length
$\xi$ diverges if measured in units of the lattice spacing. We then perform
several simulations using different values of $N\equiv L/a$; for each couple
$(\lambda,N)$ we obtain a particular value of $\beta_{c}(\lambda,N)$ such that
the condition (9) is satisfied. $\beta_{c}(\lambda)$ is then obtained by
extrapolating our results to $a/L\to 0$. As explained and numerically
demonstrated in Bosetti  _et al._ (2015), Renormalisation Group arguments
ensure us that for small enough values of $a/L$ we can extrapolate
$\beta_{c}(\lambda,a/L)$ linearly in $a/L$.

Using the relations in (8) we compute $g(\lambda,\beta_{c})$ and
$\mu_{0}^{2}(\lambda,\beta_{c})$. Then, using the renormalization condition
(5), we get the value of $\mu^{2}(g)$ and hence the ratio $f\equiv g/\mu^{2}$.
This procedure is repeated for several values of $\lambda$ (and hence of $g$)
and finally, in order to obtain $f_{0}$, we extrapolate our results to $g\to
0$.

We now focus on the condition (9). In this work we introduce a slight
modification in the procedure for the computation of the mass parameter $m$:
it is implicitly defined by the condition

$R_{\rho}\equiv\dfrac{G_{\rho}(\tau,p^{*})}{G_{\rho}(\tau,0)}=\dfrac{m^{2}}{p^{*2}+m^{2}},$
(10)

where $p^{*}$ is the smallest momentum on the lattice. We decide to compute
the propagator $G_{\rho}(\tau,p)$ using the gradient flow technique in the
contest of scalar field theory (see for example Monahan and Orginos (2015)).
In particular, considering the action (7), we introduce a new scalar field
$\rho(x,\tau)$ depending on the space–time index $x$ and on the so–called
flow–time $\tau$. The flow–time evolution equation of $\rho(x,\tau)$ is

$\dfrac{\partial}{\partial\tau}\rho(x,\tau)=\partial^{2}\rho(x,\tau)$ (11)

where $\partial^{2}$ is the Laplace operator acting in the configurations
space. If we now impose the Dirichlet boundary conditions, that is
$\rho(x,0)=\varphi_{x}$, it is easy to write the exact solution, in the
momentum space, for the propagator of the field $\rho$ at flow–time $\tau$:

$G_{\rho}(\tau,\rho)=e^{-2\tau p^{2}}G(p).$ (12)

In this way we obtain a smearing effect of the original fields, since the
flow–time exponentially suppresses the ultraviolet modes. To the total
flow–time $\tau$ we can associate a smearing radius $r_{sm}=\sqrt{2d\tau}$,
where $d$ is the dimensionality of the space–time. For a certain value of
$\tau$, the ultraviolet suppression effect of the flow–time helps us to obtain
values which are closer to the continuum limit. In this way, at fixed
$\lambda$, we expect to safely extrapolate
$\beta_{c}(\lambda,a/L)\to\beta_{c}(\lambda,0)$ using not too large values of
$L/a$ with a consequently reduction of both the computational time and the
statistical errors.

Some preliminary simulations convinced us to assume the condition
$z_{\rho}=mL=1$ and to fix the value of $\tau$ such that at different $a/L$
values the smearing radius is equal to $L/4$. In this way we take advantage of
the smearing effect of the gradient flow and take under control finite volume
effects. In Fig. 1 we show the extrapolation of $\beta_{c}(L)$ at
$\lambda=0.25$: the blue line is the extrapolation obtained by fixing
$z_{\rho}=1$, while the red line is the one obtained fixing $z=4$ and that we
used in our previous work Bosetti  _et al._ (2015). As one can see, the new
results (triangular points) show a linear behaviour even for small lattice
sizes.

At fixed value of $\lambda$ we simulate the system for several values of
$L/a$, namely: $L/a=$ 32, 40, 48, 56, 72, 80, 96, 112, 128, 144, 192, 256,
320, 384, 448, 512. We perform our simulations for the coupling values
$\lambda=$ 0.005, 0.004, 0.003, 0.002, 0.001, 0.00075, 0.0005. The number of
thermalisation sweeps for all our simulations is several hundreds times the
autocorrelation time of (9). We keep under control the autocorrelation time of
our observables using a Python program described in Palma  _et al._ (2018),
based on Wolff (2004). We perform 1000 worm–sweeps between two consecutive
measurements and the numbers of measurements varies from $10^{4}\div 10^{5}$,
according to $L$ and $\lambda$.

Figure 1: Comparison between our previous results(lower set of points) and the
new ones at $\lambda=0.25$.

##  II Results

Table 1:  Infinite volume results of $\beta_{c}$ and $f(\lambda)$ with different linear lattice sizes. $\lambda$ | $\beta_{c}$ | fit type | $(L/a)_{\textrm{min}}$ | $(L/a)_{\textrm{max}}$ | $f(\lambda)$
---|---|---|---|---|---
$0.0005$ | $0.5019535(5)$ | linear | $192$ | $512$ | $10.9920(88)$
$0.00075$ | $0.5027800(10)$ | linear | $144$ | $384$ | $10.9754(79)$
$0.001$ | $0.5035613(10)$ | linear | $112$ | $384$ | $10.9258(79)$
$0.0015$ | $0.5050340(2)$ | linear | $112$ | $384$ | $10.8801(20)$
$0.002$ | $0.5064156(9)$ | linear | $112$ | $384$ | $10.8304(45)$
$0.003$ | $0.5089871(7)$ | linear | $128$ | $384$ | $10.7504(15)$
$0.004$ | $0.5113712(4)$ | linear | $128$ | $384$ | $10.6884(15)$
$0.005$ | $0.5136155(5)$ | linear | $112$ | $512$ | $10.6435(10)$

In Tab. 1 we report the infinite volume results at the several $\lambda$ we
simulate. These results are obtained performing a linear extrapolation in
which the smallest values of $L$ are excluded. All the extrapolations give a
final $\chi^{2}/d.o.f\sim 1$. We finally extrapolate $f(\lambda)\to f_{0}$ as
$\lambda\to 0$ by using both a linear fit functions (using only the last four
points) and a quadratic fit function using all points. Our final results are

$\displaystyle f_{0}^{quad}$ $\displaystyle=11.058(4)$
$\displaystyle\chi^{2}=0.71\qquad\text{with }5\text{ d.o.f},$ (13)
$\displaystyle f_{0}^{lin\,}$ $\displaystyle=11.053(13)$
$\displaystyle\chi^{2}=1.79\qquad\text{with }2\text{ d.o.f}.$ (14)

We take the almost perfect agreement of these two results as a numerical
evidence of the fact that a simple polynomial law is enough to describe the
behaviour of $f(\lambda)$ near $\lambda\to 0$.

We decide to take the mean value of the two results as our final value and the
difference between (13) and (14) as an estimate of the systematic error
involved in the extrapolation, adding it in quadrature to the statistical
error. We finally quote

$f_{0}=11.055(14)$ (15)

to be compared to our previous result:

$f_{0}=11.15(6)(3)$ (16)

Considering both the statistical and the systematic error in (16), the two
values are well compatible within 2–sigma level but the new result (13) has a
sensibly reduced error.

Figure 2: Final extrapolation of $f$ versus $g$.

##  III Conclusions

In table 2 we summarize some of the latest results of $f_{0}$ derived with
different approaches: the works Harindranath and Vary (1988); Lee  _et al._
(2001); Sugihara (2004); Milsted  _et al._ (2013); Elias-Miro  _et al._
(2017a, b) are based on Hamiltonian truncation (variational) methods, Borel
summability is applied in Serone  _et al._ (2018); Pelissetto and Vicari
(2015), where in Pelissetto and Vicari (2015) lattice results are used.
Finally in Wozar and Wipf (2012) lattice theory is simulated by using
non–local SLAC derivative. Since we are in a good agreement with our previous
result, the same considerations are still valid: our result is compatible with
the last six determinations (excluding Burkardt  _et al._ (2016); Bosetti  _et
al._ (2015)) at the 2$\sigma$–level, which come from different methods.

The gradient flow technique allows us to reach lower values of $\lambda$ with
respect to our previous work and to obtain a more precise estimation of
$f_{0}$.

Table 2: Sample of the results for the continuum critical parameter $f_{0}$ from the literature. DLCQ stands for _Discretized Light Cone Quantization_ , QSE diagonalization for _Quasi–Sparse Eigenvector_ diagonalization, DMRG for _Density Matrix Renormalization Group_ and for DLCH-FS _Diagonalized light-front Hamiltonian in Fock-Space representation_. Method | $f_{0}$ | year, Ref.
---|---|---
DLCQ | $5.52$ | 1988, Harindranath and Vary (1988)
QSE diagonalization | $10$ | 2000, Lee  _et al._ (2001)
DMRG | $9.9816(16)$ | 2004, Sugihara (2004)
Monte Carlo cluster | $10.8^{0.1}_{0.05}$ | 2009, Schaich and Loinaz (2009)
Monte Carlo SLAC derivative | $10.92(13)$ | 2012, Wozar and Wipf (2012)
Uniform Matrix product states | $11.064(20)$ | 2013, Milsted  _et al._ (2013)
Monte Carlo worm | $11.15(6)(3)$ | 2015, Bosetti  _et al._ (2015)
Borel summability–Lattice Results | $11.00(4)$ | 2015, Pelissetto and Vicari (2015)
DLCH-FS | $4.40(12)$ | 2016, Burkardt  _et al._ (2016)
Renormalised Hamiltonian | $11.04(12)$ | 2017, Elias-Miro  _et al._ (2017a, b)
Borel summability | $11.23(14)$ | 2018, Serone  _et al._ (2018)
This work | $11.055(14)$ | 2018

## References

  * Loinaz and Willey (1998) W. Loinaz and R. S. Willey, Phys. Rev. D 58, 076003 (1998).
  * Schaich and Loinaz (2009) D. Schaich and W. Loinaz, Phys.Rev. D79, 056008 (2009), arXiv:0902.0045 [hep-lat] .
  * Korzec  _et al._ (2011) T. Korzec, I. Vierhaus, and U. Wolff, Computer Physics Communications 182, 1477 (2011).
  * Bosetti  _et al._ (2015) P. Bosetti, B. De Palma, and M. Guagnelli, Phys. Rev. D 92, 034509 (2015).
  * Monahan and Orginos (2015) C. Monahan and K. Orginos, Phys. Rev. D 91, 074513 (2015).
  * Palma  _et al._ (2018) B. D. Palma, M. Erba, L. Mantovani, and N. Mosco, Computer Physics Communications (2018), https://doi.org/10.1016/j.cpc.2018.07.004.
  * Wolff (2004) U. Wolff (ALPHA), Comput. Phys. Commun. 156, 143 (2004), arXiv:hep-lat/0306017 [hep-lat] .
  * Harindranath and Vary (1988) A. Harindranath and J. P. Vary, Phys. Rev. D 37, 1076 (1988).
  * Lee  _et al._ (2001) D. Lee, N. Salwen, and D. Lee, Physics Letters B 503, 223 (2001).
  * Sugihara (2004) T. Sugihara, JHEP 0405, 007 (2004), arXiv:hep-lat/0403008 [hep-lat] .
  * Milsted  _et al._ (2013) A. Milsted, J. Haegeman, and T. J. Osborne, Phys. Rev. D 88, 085030 (2013).
  * Elias-Miro  _et al._ (2017a) J. Elias-Miro, S. Rychkov, and L. G. Vitale, Phys. Rev. D96, 065024 (2017a), arXiv:1706.09929 [hep-th] .
  * Elias-Miro  _et al._ (2017b) J. Elias-Miro, S. Rychkov, and L. G. Vitale, JHEP 10, 213 (2017b), arXiv:1706.06121 [hep-th] .
  * Serone  _et al._ (2018) M. Serone, G. Spada, and G. Villadoro, (2018), arXiv:1805.05882 [hep-th] .
  * Pelissetto and Vicari (2015) A. Pelissetto and E. Vicari, Physics Letters B 751, 532 (2015).
  * Wozar and Wipf (2012) C. Wozar and A. Wipf, Annals Phys. 327, 774 (2012), arXiv:1107.3324 [hep-lat] .
  * Burkardt  _et al._ (2016) M. Burkardt, S. S. Chabysheva, and J. R. Hiller, Phys. Rev. D 94, 065006 (2016).