# Superconducting Phases in Neutron Star Cores
## 3 Superconducting Ground States
### The Numerical Model  
Whether we choose to work with the Helmholtz free energy, \(\mathcal{F}\), or with the Gibbs free energy, \(\mathcal{G}\), we obtain the same system of Euler-Lagrange equations:  
\[\kappa^{2}\mathbf{\
abla}\times(\mathbf{\
abla}\times\mathbf{A})=\mathrm{ Im}\bigg{\{}\psi_{\mathrm{P}}^{*}(\mathbf{\
abla}-\mathrm{i}\mathbf{A})\psi_{ \mathrm{P}}+\frac{h_{1}}{\epsilon}\psi_{\mathrm{h}}\psi_{\mathrm{P}}^{*}(\bm {\
abla}-\mathrm{i}\mathbf{A})(\psi_{\mathrm{n}}^{*}\psi_{\mathrm{P}})\bigg{\}}\,, \tag{16}\] \[\
abla^{2}\psi_{\mathrm{n}}=R^{2}(|\psi_{\mathrm{n}}|^{2}-1)\psi_ {\mathrm{n}}+\alpha(|\psi_{\mathrm{P}}|^{2}-1)\psi_{\mathrm{n}}\] \[-h_{1}\psi_{\mathrm{P}}(\mathbf{\
abla}+\mathrm{i}\mathbf{A})^{2}( \psi_{\mathrm{p}}^{*}\psi_{\mathrm{h}})-\psi_{\mathrm{n}}\
abla^{2}\bigg{(} \frac{h_{2}-h_{1}}{2}|\psi_{\mathrm{P}}|^{2}+\frac{h_{3}}{2\epsilon}|\psi_{ \mathrm{n}}|^{2}\bigg{)}\,,\] (17) \[(\mathbf{\
abla}-\mathrm{i}\mathbf{A})^{2}\psi_{\mathrm{P}}=(|\psi_{ \mathrm{P}}|^{2}-1)\psi_{\mathrm{P}}+\frac{\alpha}{\epsilon}(|\psi_{\mathrm{n} }|^{2}-1)\psi_{\mathrm{P}}\] \[-\frac{h_{1}}{\epsilon}\psi_{\mathrm{n}}(\mathbf{\
abla}-\mathrm{i} \mathbf{A})^{2}(\psi_{\mathrm{n}}^{*}\psi_{\mathrm{P}})-\psi_{\mathrm{P}}\
abla ^{2}\bigg{(}\frac{h_{2}-h_{1}}{2\epsilon}|\psi_{\mathrm{n}}|^{2}+\frac{h_{3}}{ 2}|\psi_{\mathrm{P}}|^{2}\bigg{)}\,. \tag{18}\]  
However, the appropriate boundary conditions for these two experiments are different, and also depend on the particular size and shape chosen for the domain. Without loss of generality, we will assume from here on that the magnetic field is oriented in the \(z\)-direction, and that all variables are independent of \(z\); so our domain will be some region within the \(xy\)-plane. We solve a discretized version of Equations (16)-(18), which are obtained by minimizing a discrete approximation to the free energy on a regular grid in \(x\) and \(y\). The gauge field is included via a Peierls substitution, in order to maintain gauge invariance. The equations are solved using a simple relaxation method, and the grid resolution is repeatedly refined until a sufficient level of accuracy has been obtained. Additional details on the numerical algorithm can be found in Appendix A. In the present study, we are not interested in the effect of physical boundaries on the phase diagram, and so we would ideally use an infinite domain, but for numerical calculations the domain must of course be finite. Moreover, we cannot use periodic boundary conditions, because in the presence of fluxtubes neither \(\mathbf{A}\) nor \(\psi_{\mathrm{p}}\) is spatially periodic. Instead, we must use quasi-periodic boundary conditions ([PERSON] and [PERSON], 1998), which involves specifying not only the size of the domain, \(L_{x}\times L_{y}\), say, but also the number, \(N\), of magnetic flux quanta within the domain. Working in the symmetric gauge, the quasi-periodic boundary conditions for our dimensionless variables take the form  
\[\mathbf{A}(\mathbf{x}+\mathbf{L}) =\mathbf{A}(\mathbf{x})+\frac{N\pi}{L_{x}L_{y}}\mathbf{e}_{z} \times\mathbf{L}\,, \tag{19}\] \[\psi_{\mathrm{P}}(\mathbf{x}+\mathbf{L}) =\psi_{\mathrm{P}}(\mathbf{x})\exp\bigg{(}\mathrm{i}\frac{N\pi}{L_ {x}L_{y}}\mathbf{e}_{z}\times\mathbf{L}\cdot\mathbf{x}\bigg{)}\,,\] (20) \[\psi_{\mathrm{n}}(\mathbf{x}+\mathbf{L}) =\psi_{\mathrm{n}}(\mathbf{x})\,, \tag{21}\]  
where \(\mathbf{L}\) represents either of the translational symmetries \((L_{x},0)\) or \((0,L_{y})\). The effect of these boundary conditions is to impose a total magnetic flux of \(2\pi N\) within our domain (in our dimensionless units, the quantum of magnetic flux is \(2\pi\)), and thus a mean magnetic flux of \(\overline{B}=2\pi N/(L_{x}L_{y})\). By solving Equations (16)-(18) subject to these boundary conditions, in domains of varying size, we can thus obtain the Helmholtz free energy, \(\mathcal{F}\), as a function of the mean magnetic flux. We note that the choice of domain aspect ratio affects the configuration of any fluxtube array that forms.