Hartree-Fock-Bogoliubov theory for odd-mass nuclei with a time-odd constraint and application to deformed halo nuclei

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . We show that the lowest-energy solution of the Hartree-Fock-Bogoliubov (HFB) equation is the particle-number parity being even as long as the time-reversal symmetry is conserved in the HFB Hamiltonian. Based on this finding, we give a rigorous foundation of a method for solving the HFB equation to describe the ground state of odd-mass nuclei by employing a timereversal anti-symmetric constraint operator to the Hamiltonian, where one obtains directly the ground state as a self-consistent solution of the cranked-HFB-type equation. Numerical analysis is done for the neutron-rich Mg isotopes with a reasonable choice for the operator, and it is demonstrated that the anomalous increase in the matter radius of Mg is well described combining with the framework of the nuclear energy-density-functional method, revealing the deformed halo structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Subject Index D10, D11, D13


Introduction
Odd nuclei, composed of the odd number of either or both of protons and neutrons, unveil the unique features that one cannot observe in even-even nuclei. Even-even nuclei, for instance, have spin J π = 0 + in the ground state, where nucleons are paired off due to the correlation, while odd nuclei can have nonzero ground-state spin, where the last nucleon does not take part in the pair correlation and is responsible for the total spin. The spin gives us information of single-particle orbitals near the Fermi level. With increase in the structure information of unstable nuclei thanks to the recent advancement in radioactiveisotope beam technology [1], more and more exotic features in odd nuclei have been showing up: Highlights in the latest discoveries include the shape staggering in 181-185 Hg [2], and the deformed halo structures of 31 Ne [3,4] and 37 Mg [5,6]. A great theoretical challenge under these circumstances is to describe odd nuclei in a wide mass region of nuclear chart, where the pair correlation, the shape deformation, and the weak binding of nucleons are all considered in a unified manner.
Nuclear density-functional theory (DFT) or the self-consistent mean-field model has been extensively employed for describing the systematic features of not only the ground but excited states [7,8]. The nuclear landscape has been investigated in the framework of both nonrelativistic and relativistic energy-density functional (EDF) methods [9,10]. The Hartree-Fock-Bogoliubov (HFB) theory or the Kohn-Sham-Bogoliubov-de-Gennes scheme in DFT is capable of providing us a unified description of the ground-state properties for not only eveneven nuclei but odd nuclei taking the superfluidity and shape deformation into account [11].
In spite of the successful application of DFT, the calculations have been mostly restricted to even-even nuclei, and the odd nuclei remain largely unexplored [12]. This may be partly because the primal interest has been put on determining the drip lines [13]. With the recent advent of computational resources sufficient to perform global calculation in the framework of DFT, systematic odd-even alternations in atomic nuclei such as the odd-even staggering of the binding energies have attracted renewed interest [14].
There seem many obstacles to tackling the systematic investigation of odd nuclei in DFT, and some of them are the followings: (i) Coexistence of multiple levels in low energy; the high precision calculation, at the same time with high accuracy, is required for resolving the near degeneracy of several levels and identifying the ground state. This is not restricted to DFT but a challenge for any theoretical models. (ii) Non-vanishing time-odd densities; an EDF is a time-even scalar constructed from various densities, and includes the densities and currents that are odd with respect to time reversal to preserve the Galilean invariance and to describe properly e.g. the spin-dependent observables. While the time-odd densities automatically vanish for the ground state of even-even nuclei, they are nonzero in odd nuclei where the timereversal symmetry is intrinsically broken [7]. Allowing the breaking of time-reversal symmetry increases the computational cost [15]. Furthermore, EDFs commonly used in the practical calculations are phenomenologically constructed by using the properties of time-even states only. Thus, the coupling constants of time-odd densities are highly uncertain. (iii) Complexity of the blocking procedure; one cannot usually obtain the ground state of odd nuclei as the lowest energy solution and needs an additional procedure to excite one quasiparticle on top of the ground state of even-even nuclei [16]. Therefore, pragmatic techniques are needed in an actual calculation [17,18]. We are going to focus on (iii).
In this study, we show that the Bogoliubov transformation, the particle-number parity, and the time-reversal symmetry in HFB are closely related with one another. From these findings, we can give a rigorous foundation for a method initiated by Bertsch et al. [19] to describe the ground state of odd-particle system as the lowest energy state under an appropriate time-odd constraint in HFB theory. This method has a high affinity with DFT in the sense that either an odd-particle system or an even-particle system is described as the ground state uniformly. We apply this method to the neutron-rich Mg isotopes near the drip line and demonstrate that it produces the exotic behavior in radii observed experimentally.
The article is organized as follows: In Sec. 2, after recapitulating the basics of HFB theory, the relationships among the Bogoliubov transformation, the particle-number parity, and the time-reversal symmetry are presented. In Sec. 3, based on the relationships found in Sec. 2, we give foundation for a method describing the ground-state of odd-nuclei under a timeodd constraint. In Sec. 4, we give a numerical procedure for describing the weakly-bound neutron-rich nuclei by employing the Skyrme-type EDF with the inclusion of the time-odd fields. Then, results of the numerical analysis for Mg isotopes are presented. Finally, summary is given in Sec. 5.
Part of the preliminary results of this work are reported in Ref. [20].

Basics of HFB theory
We begin with recalling the basics of Hartree-Fock-Bogoliubov (HFB) theory. The notation used here follows Ref. [11]. The Bogoliubov quasiparticle (qp) creation and annihilation operatorsβ † k ,β k are defined as linear combinations of single-particle (sp) creation and annihilation operatorsĉ † k ,ĉ k : where indices k and l run over the whole configuration space (k = 1, . . . , M). Since we consider spin-1/2 particles, M is an even number. The Bogoliubov transformation between the qp and sp bases is represented by the 2M × 2M matrix as β In order to satisfy the fermion anticommutation relations for quasiparticles, W must be unitary: W † W = WW † = I 2M with I n representing the n × n identity matrix. The ground state wave function of the many-body system in HFB theory, or the HFB vacuum, |Φ is defined as the vacuum of the quasiparticles: The complete information about |Φ is contained by the density matrix ρ kl := Φ|ĉ † lĉ k |Φ and pairing tensor κ kl := Φ|ĉ lĉk |Φ , or by the generalized density matrix The unitarity of W guarantees that R is idempotent: R 2 = R. Under the idempotency of R, the variational principle with a particle-number constraint: δ Φ|Ĥ − λN |Φ = 0, whereĤ, N, and λ are the Hamiltonian of the system, the particle-number operator, and the chemical potential, respectively, leads to that R commutes with the HFB Hamiltonian where h and ∆ are the sp and pair Hamiltonian, respectively. It follows that R and H are simultaneously diagonalized, and thus the HFB equations are represented in a matrix form or where E := diag(E, −E) and E is a diagonal matrix of qp energies E k . Note that the HFB Hamiltonian inherently has the particle-hole symmetry: 4 where Σ x := 0 1 1 0 and O n represents the n × n zero matrix. It follows that when ϕ k is an eigenvector of H with eigenvalue E k , Σ x ϕ * k is also an eigenvector with eigenvalue −E k . Thus, the eigenvalues of H always come in pairs of opposite signs, but the theory says nothing about the individual signs of E k . Therefore, we have to choose for each k (k = 1, · · · , M) whether to take E k positive or negative. A naïve choice to describe the ground state of the system is to take all E k positive. Indeed, the state obtained by this choice is the lowest energy state in the sense that all the excitation energies are positive.
While the particle number is no longer conserved in HFB theory, the particle-number parity, or we call simply this as number parity, that is whether the number of constituent particles is even or odd, is still a good quantum number [11]. In what follows, we heuristically derive the explicit relations among the number parity π N , the Bogoliubov transformation matrix W, and the time-reversal symmetry of the HFB Hamiltonian H.

Particle-number parity and the Bogoliubov transformation
We first show a relationship between the number parity π N and the Bogoliubov transformation matrix W. The famous theorem of Bloch and Messiah says that a unitary matrix W of the form (2) can always be decomposed into three matrices of very special form [21]: Here C and D are unitary matrices andŪ,V are real matrices of the general form (1) 0 U (2) . . .
andŪ (p) ,V (p) are 2 × 2 matrices of the form where u p and v p satisfy the conditions: u p > 0, v p > 0, u 2 p + v 2 p = 1. One can explicitly construct the HFB vacuum |Φ in terms of the canonical basis defined asâ † k = l D lkĉ † l : Here |0 is the empty state defined asĉ k |0 = 0 for all k. The indexp represents an orbital paired with the orbital p, and N p indicates the maximum number of pairs in |Φ . N 1 represents the number of unpaired particles, and corresponds to the seniority number in the quasi-spin theory [22]. Depending on whether N 1 is even or odd, |Φ is a superposition of states with either even or odd particles, but never both. This means that the HFB vacuum |Φ is an eigenstate of the operatorP N = e iπN , whereN is the particle-number operator. The eigenvalue π N is a good quantum number, called the number parity [11].
where we used the fact that C and D are unitary matrices andW (p) :=Ū (p) +V (p) = (Ū (p) − V (p) ) T are orthogonal matrices. It follows that det W is nothing but the number parity: The fact that the determinant of W is either +1 or −1 implies that W can be unitarily transformed into an orthogonal matrix. In fact, for the unitary matrix W X := X WX † is a real orthogonal matrix. The HFB equations are transformed by the same X into a skew symmetric form of where H X := X HX † is a pure-imaginary skew-symmetric matrix and Since the both sides of (18) are 2M × 2M skew-symmetric matrices, we can take Pfaffian of the both sides. Then we obtain 1 For the lowest energy state obtained by taking all E k positive, one sees

Symmetries in HFB theory
Before investigating the time-reversal symmetry of the HFB Hamiltonian, we are going to give a general property of the symmetry in HFB theory. Since the HFB equations are nonlinear, the HFB Hamiltonian H does not necessarily hold the same symmetries as the Hamiltonian of the systemĤ. Nevertheless, certain symmetries are still conserved in HFB theory. Such symmetries are called self-consistent symmetries, and they often significantly reduce the dimension of the eigenvalue problem [23]. 1 Here we used the Pfaffian identities: For a 2n × 2n skew-symmetric matrix A and an arbitrary 2n × 2n matrix B, For an arbitrary n × n matrix C, Consider a symmetry transformation realized by a unitary or anti-unitary operatorÛ s which maps the sp space into itself by a M × M unitary matrix U s : or in the 2M-dimensional spacê where Here (· · · ) ( * ) denotes that the complex conjugate is taken ifÛ s is an anti-unitary operator.
Assuming thatÛ s is a symmetry operator of the system, that is [Ĥ,Û s ] = 0, we find the HFB Hamiltonian H transforms in the same way as the generalized density matrix R: Now suppose that the HFB vacuum |Φ is invariant up to a phase under the operationÛ s , i.e.Û s is a symmetry operator of the intrinsic system, it follows that Specifically, we are interested in suchÛ s that is generated by a hermitian particle-hole one-particle operatorŜ = kl S klĉ † kĉ l :Û s = e iθŜ or e iθŜK , where θ is a real parameter andK is the complex conjugation operator which leaves the sp basis |k = c † k |0 invariant:K |k = |k . In this case, U s = e iθS . Especially whenÛ s = e iθŜ , Eq. (26) leads to If this holds for an arbitrary θ, one sees Note that since the HFB vacuum is always an eigenstate of the number parity operator

Time-reversal symmetry and number parity
Using the results obtained above, we show that the time-reversal symmetry of the HFB Hamiltonian and the number parity of the HFB vacuum are directly related to each other.
Let us consider the case when the HFB Hamiltonian has time-reversal symmetry. Paying attention to the anti-unitarity of the time-reversal operatorT = exp(−iπŜ y )K, whereŜ y is the y-component of the total spin operator, we have the following equality relation from (26) where It follows that when ϕ k is an eigenvector of H with eigenvalue E k , the time-reversed state T ϕ * k is an independent eigenvector of H with the same eigenvalue E k . In other words, ϕ k and T ϕ * k are degenerated, i.e. the Kramers degeneracy shows up. Since the particle-hole symmetry (9) is always kept in HFB theory, Σ x ϕ * k and Σ x T ϕ k are also Kramers-degenerated eigenvectors with eigenvalue −E k . Therefore, the following unitary matrix can be taken as a solution of the HFB equations for H that satisfies (29): Then, the HFB equations read Multiplying U by T from the left and T T from the right, we get In the same way, one obtains T V T T = V * . Thus, it follows that T WT T = W * . Then, multiplying the both sides by W T from the left and T from the right and using the fact that W is a unitary matrix and T is a orthogonal matrix, we obtain This indicates that W is a symplectic matrix. It is known that the determinant of a symplectic matrix is +1. In fact, taking Pfaffian of the both sides of (37), it follows that det W pf T = pf T . And since pf T = (−1) M = 0, one sees Therefore, the Pfaffian of the X -transformed H satisfying (29) is given by The number parity of the lowest energy state, i.e. det W L.E. , obtained by the time-reversal symmetric H is then by substituting (39) into (21). Since both the numerator and the denominator of the right hand side of (40) are positive, we arrive at det W L.E. = +1.

Methodology for describing odd nuclei
We demonstrated above that the lowest energy state is always an even number-parity state, that is, an even particle system in HFB theory as long as the time-reversal symmetry is conserved for the intrinsic Hamiltonian. The procedure called the blocking method, which has been conventionally used to describe an odd number-parity state [12,16], can be viewed as follows: One solves the time-reversal symmetric HFB equations to generate a reference state with det W = +1, and then swaps one set of columns of W to obtain a state whose sign of determinant is reversed.
Alternatively, one can take a strategy to describe an odd number-parity state as the lowest energy state under the constraint which breaks time-reversal symmetry. The idea of obtaining odd nuclei by such a constraint was proposed by Bertsch et al., stimulated by the non-collective cranking method [19]. In what follows, we generalize this approach and encapsulate the essential point of this method.
Assume that the intrinsic system is invariant under a unitary transformation generated by a hermitian time-odd particle-hole type one-body operatorŜ = ij S ij c † i c j with S † = S, and T S * T T = −S. This means that the mean-field representation ofŜ commutes with the HFB Hamiltonian H: where Introducing a Lagrange multiplier λ s forŜ and constraining its expectation value, one minimizes the Hamiltonian under the constraint, or the Routhian; Then we can write the HFB Routhian as From the commutation relation (41), H can be block diagonalized for each eigenvalue of S.
Since S is proportional to the identity matrix in each block, the eigenvalues of H ′ are linearly shifted from those of H. Then, as shown below, the sign of pf H ′ X can change according to λ s and thus the number parity of the lowest energy state can vary from positive to negative.
Because of the time-odd character ofŜ, the eigenvalues of S always come in pairs of opposite signs. Let {x ±α n } n=1,2,··· be sets of orthonormalized eigenvectors of S with eigenvalue ±ω α (α > 0, ω α > 0), where n is a label that distinguishes states other than α. Then , 0 is a set of orthonormalized eigenvectors of S with eigenvalue ω α . An appropriate linear combination is a simultaneous eigenstate of S and H with an eigenvalue ω α and E α n , respectively. Thanks to the particle-hole symmetry of H, is also a simultaneous eigenstate of S and H with an eigenvalue ω α and −E −α n , respectively. Therefore, the HFB equations are block diagonalized for each eigenvalue of S as follows: where and Since χ α n and ϕ α n are simultaneous eigenstates of S and H, they are also eigenstates of H ′ with an eigenvalue E α n − λ s ω α and −E −α n − λ s ω α , respectively. In other words, where H ′α = H α − λ s ω α I α , and I α is the identity matrix for the block with eigenvalue ω α . In this way, one sees that the constraint on the intrinsic symmetryŜ does not change individual single qp states, but shifts only the qp energies according to the eigenvalues of S. Reflecting the time-odd character ofŜ, the qp energies of time-reversal pairs split in the opposite direction. Therefore, even if the original HFB Hamiltonian H has time-reversal symmetry, the Kramers degeneracy is resolved at λ s = 0. Then, we are going to show that the number parity of the lowest energy state can change according to λ s . The following unitary matrix diagonalizes H and H ′ simultaneously: Then the HFB equations for H and H ′ read where One sees for the lowest energy states Substituting this into (56), one sees det W = +1.
Substituting further this into (57), one sees Therefore from (59) we arrive at At λ s = 0, det W ′ L.E. = det W L.E. = +1, namely the lowest energy state has even number parity. For λ s = 0, the sign of det W ′ L.E. changes according to the magnitude relation between E α n and λ s ω α . The number parity of the lowest energy state thus changes according to the magnitude of λ s . 1 Schematic picture of the single-qp energies as functions of the parameter λ s for the cases employing the operator (a)Ĵ z , (b)Ĵ z /|J z |, and (c) |J z = 1/2 J z = 1/2| − |J z = −1/2 J z = −1/2| as the constraint operatorŜ. At λ s = 0, the single-qp levels with Ω are displayed as E Ω . In increasing λ s , the determinant of the lowest energy state keeps +1 for λ s < λ 1 , the sign changes at λ s = λ 1 , and the determinant is −1 for λ 1 < λ s < λ 2 . Figure 1 shows a schematic picture of the single-qp energies under the constraint. When λ s is set to an appropriate value λ 1 < λ s < λ 2 , here a pair of levels intersect the axis at λ s = λ 1 , an odd number-parity state is automatically obtained as the lowest energy state.

13
Determined by the operatorŜ are the orbital whose qp energy changes with an increase in λ s and the level which intersects the axis first. In this senseŜ is considered as a selector of the vacuum. Note that when two levels cut across the axis by increasing λ s , the number parity of the system becomes even corresponding to the two-qp excitation state.
Specifically, when the z-component of the total angular momentumĴ z is a symmetry of the intrinsic system, different vacuums are selected depending on the choice ofŜ. We here consider three cases. Case (i); whenĴ z itself is taken asŜ, each qp level rises or falls with the slope of the corresponding eigenvalue J z = Ω, see Fig. 1(a). As λ s increases, the orbital with high-Ω is preferably selected as the level to be excited. This is equivalent to the so-called noncollective cranking, where particles with high angular momentum about the symmetry axis are aligned. Case (ii); let us look at the case of usingĴ z /|J z | asŜ. SinceĴ z /|J z | is an operator that gives +1 for eigenstates with positive Ω and −1 for eigenstates with negative Ω, each qp level separates out with the slope of ±1 independent of the magnitude of Ω with an increase in λ s , see Fig. 1(b). Therefore, in this case, an orbital with a smaller qp energy is likely to be selected as the level to be excited. This corresponds to the two-Fermi level approach proposed in Ref. [19]. As mentioned in Ref. [19], this choice may not work well for the case in which the single qp level density is high near the Fermi level and the spherical systems where we have (2j + 1)-folded degeneracy. Introducing a kind of projection operator was conjectured in Ref. [19] to resolve the issue and we realize the practical method as Case (iii): Let us consider the case of using a time-odd projection operator |J z = Ω J z = Ω| − |J z = −Ω J z = −Ω| asŜ. This is a part ofĴ z /|J z |, and is an operator which works only for the state with a certain Ω. Therefore, in this case, only the level carrying the specified eigenvalue ±Ω splits for λ s = 0, and the levels having other eigenvalues ofĴ z do not change even if λ s increases, see Fig. 1(c). Thus, one can select the state of interest easily. The third choice is convenient for the practical use and we perform the calculations using this choice in the following investigation.

HFB equation for axially-symmetric nuclei with time-even and time-odd mean fields
The coordinate-space HFB equation obtained by employing the local EDF containing time-even and time-odd parts reads [24,25] where q stands for protons (p) and neutrons (n) in which the quasiparticles are assumed to be eigenstates of the third component of the isospin operator. The sp Hamiltonian h consists of the mean-field (Kohn-Sham) potentials Γ t composed of the time-even and time-odd isoscalar (t = 0) and isovector (t = 1) densities as Here, V Coul is the Coulomb potential, and the explicit expressions of Γ for the Skyrme-type EDF are shown in Appendix. Thanks to the time-reversal (anti-)symmetry of the potentials, one seesh We employ the pairing EDF that contains only the time-even densities as described below, so that we see We solve the HFB equation (64) by assuming the axial and reflection symmetries so that the quasiparticles are labeled by {Ω, π, q}, with Ω and π being the z-component of the total angular momentum and parity, respectively. In this case, the qp wave functions can be written in the form of ϕ q a,nΩπ (r σ) = ϕ q+ a,nΩπ (̺, z) e iΛ − φ χ 1/2 (σ) + ϕ q− a,nΩπ (̺, z) e iΛ + φ χ −1/2 (σ), (a = 1, 2), (70) where Λ ± = Ω ± 1/2 [26], and ̺, z, and φ are the cylindrical coordinates defining the threedimensional position vector as r = (̺ cos φ, ̺ sin φ, z), while z is the chosen symmetry axis. And, the wave functions satisfy the following symmetry The coordinate-space HFB equation has been solved under the assumption of the axial symmetry in many cases, however the time-reversal symmetry is often imposed [13,[27][28][29][30][31]. To keep the present paper self-contained, the mean-field potentials containing both the time-even and time-odd densities and currents in the cylindrical-coordinate representation are shown in Appendix. With the axial symmetry, the φ dependences of the sp and pair Hamiltonians are given by andh where l z = where h Ωπq ss ′ ,h Ωπq ss ′ , andh Ωπq ss ′ are defined by using Λ ↑ := Λ − and Λ ↓ := Λ + as h Ωπq ss ′ (̺, z) := h q ss ′ (̺, z; l z = Λ s ′ ),h Ωπq ss ′ (̺, z) :=h q ss ′ (̺, z; l z = Λ s ′ ), andh Ωπq ss ′ (̺, z) :=h q ss ′ (̺, z; l z = Λ s ′ ). To describe odd-A isotopes, we employ the time-odd projection operator to the states with a specific {Ω, π, q} quantum number, |Ω π q Ω π q| − |−Ω π q −Ω π q|, as the constraint operatorŜ. In other words, we introduce the Lagrange multiplier λ s for the specified {Ω, π, q} sector of the HFB equation (74) as We call below introducing λ s "blocking" since this procedure is equivalent to the traditional blocking method.

Numerical procedures
We solve Eq. (74) by diagonalizing the HFB Hamiltonian in the cylindrical-coordinate representation with the box boundary condition. We discretize the coordinates by ̺ i = (i − 1/2) × h (i = 1, 2, · · · , N ρ ) and z j = (j − 1) × h (j = 1, 2, · · · , N z ), with a lattice mesh size h = 0.8 fm, and use 30 points for N ρ and N z . The qp states are truncated according to the qp energy cutoff at 60 MeV, and the qp states up to Ω = 15/2 with positive and negative parities are included. The differential operators are represented by use of the 13-point formula of finite difference method. For diagonalization of the Hamiltonian or Routhian, we use the LAPACK dsyevx subroutine [32]. A modified Broyden's method [33] is utilized to calculate new densities during the self-consistent iteration. The Lagrange multiplier λ s is adjusted to and E q 2Ωπ are the lowest and the second lowest positive qp energies for the given {Ω, π, q} at each iteration so that only one pair of qp levels intersect the energy-zero axis. For the normal (particle-hole) part of a nuclear EDF, we employ the SLy4 functional [34]. The so-called naïve choice [12] is adopted for determining the coupling constants of the time-odd terms in the EDF except that the coupling constants of the terms of the form s · △s, where s is the spin density, are set to zero because the terms in some cases lead to divergences of the HFB iterative procedure [12]. For the particle-particle channel, we adopt the density functional (B1) which corresponds to the density-dependent contact interaction.
The parameters are set as V 1 = 1 and γ = 1 (surface pairing), and the pairing strength is taken as V 0 = −430 MeV fm 3 to reproduce approximately the experimental pairing gap of neutrons (1.28 MeV based on AME2016 [35]) of 35 Mg. The pairing gap is obtained by use of the three-point formula for the binding energy [36], and the calculated ∆ (3) n for 35 Mg is found to be 1.47 MeV.

Numerical results and discussion
To demonstrate the feasibility of our method, we perform the systematic calculation for the neutron-rich Mg isotopes with the mass number A = 34 -40. We exclude 32 Mg and 33 Mg in the present investigation, where the loss of spherical magic number of 20 has been under debate, because the shape fluctuation and the correlation beyond the mean-field approximation may be significant in 32 Mg [37][38][39], and many-particle many-hole states with different shape deformation may coexist in 33 Mg [40,41] as mentioned slightly below.
We tried blocking each of Ω π = 1/2 ± , 3/2 ± , 5/2 ± , and 7/2 ± orbitals for odd-mass isotopes, and the ground state was obtained by blocking the orbital with Ω π = 3/2 − , 5/2 − , and 1/2 − in 35 Mg, 37 Mg, and 39 Mg, respectively. The calculation may be in contradiction with the observation for 35 Mg, where J π = 3/2 − , a head of the K π = 1/2 − band, is suggested for the ground state [42]. We found that the binding energy obtained by blocking the Ω π = 1/2 − orbital is shallower by 1.0 MeV. For 37 Mg, the measurements suggest that the ℓ = 1 component is dominant in the ground state [5,6], while the Ω π = 5/2 − orbital contains angular momenta higher than ℓ = 3. Below, we are going to discuss 37 Mg on this point. It is noted here that the neutron superfluidity vanishes in 35,37 Mg. Figure 2(a) shows the calculated one-neutron separation energies S n compared with the experimental or evaluated data obtained from AME2016. A nice agreement within the error range can be seen. We reproduced the odd-even staggering and the instability of 39 Mg with respect to the neutron emission. It is noted that we need the binding energy of 33 Mg for the calculation of S n of 34 Mg. We obtained the near-degeneracy by blocking the 1/2 − and 7/2 − orbitals for 33 Mg. For 39 Mg, the calculated chemical potential is λ n = −0.98 MeV, and the qp energy of the blocked orbital is 0.83 MeV. Thus, the 39 Mg nucleus is bound, though quite loosely, in the present calculation. Furthermore, we could not find the bound-state solutions for 41,42 Mg.
That the deformation parameters for neutrons are lower than those for protons indicates the larger radii of neutrons than protons. The odd-even staggering in deformation cannot be recognized. Therefore, the odd-even staggering seen in the binding energy of these isotopes can be attributed to the pair correlation mainly.
Let us discuss the systematic feature in matter radii. Figure 3 shows the calculated matter radii r 2 m compared with the observation based on the reaction cross section measurement [50]. Except for 37 Mg, the present calculation reproduces the isotopic dependence observed experimentally. However, we see a systematic underestimation. This is mainly because the calculation gives a systematic over binding. The irregular dependence revealed in [34][35][36] Mg by the reaction cross section measurement [5] is well described; the matter radius of 35 Mg is smaller than the average of the radii of the neighboring isotopes of 34,36 Mg. The discrepancy between the calculation and the observation for 37 Mg is due to the suppression of spatial extension caused by the high centrifugal barrier of the blocked orbital of [312]5/2 stemming from the f 7/2 shell. As mentioned above, the experimental measurements suggest that the ground-state in 37 Mg is dominated by the p-wave [5,6]. When the deformation develops further, the [312]5/2 orbital crosses with the [321]1/2 orbital. The latter orbital is originating from the p 3/2 shell. Therefore, blocking the Ω π = 1/2 − orbital is worth Fig. 3 Calculated matter radii denoted by closed circles compared with the observation denoted by crosses with error bars taken from Ref. [50]. As in Fig. 2, the results obtained by blocking the Ω π = 1/2 − orbital are also shown for 37 Mg. The symbol for (1/2 − ) indicates the result obtained by ignoring the time-odd mean fields.
investigating. Indeed, the deformed halo structure in 37 Mg has been studied by assuming a high deformation of β 2 ∼ 0.5 in a deformed Woods-Saxon potential to put a neutron in the [321]1/2 orbital on the even-even 36 Mg nucleus [51]. We show in Fig. 3 the result obtained by blocking the Ω π = 1/2 − orbital for 37 Mg as indicated by a red circle. The matter radius increases by 0.08 fm by blocking the Ω π = 1/2 − orbital. Then, we can see a sudden enhancement from 36 Mg. As expected, the deformation develops that we can see in Fig. 2(b). The matter quadrupole deformation obtained is β 2 = 0.35, and we see an increase in β 2 by 0.03, which is far lower than the phenomenological value [51]. It is noted again that the enhancement in radius impedes the increase in deformation as far as the definition (76) is used. We mention here a way was proposed to disentangle the effects of the spatial extension and the quadrupole deformation [50]. Taking it into consideration, we are going to investigate in details the properties of neutron-rich Mg and Ne isotopes altogether in the forthcoming paper. The total binding energy calculated by blocking the Ω π = 1/2 − orbital is shallower by  the Ω π = 1/2 − and 5/2 − orbitals are shown. In the case of Ω π = 1/2 − , shown are also the result obtained by ignoring the time-odd mean fields.
0.53 MeV, which we see in Fig. 2(a). We found that the neutrons of 37 Mg obtained by blocking the Ω π = 1/2 − orbital are paired. The calculated chemical potential and the qp energy of the blocked orbital is −2.70 MeV and 2.20 MeV, respectively, while the sp energy of the last occupied neutron is −1.1 MeV obtained by blocking the Ω π = 5/2 − orbital. The asymptotic behavior of the last occupied orbital is given as rϕ 2,i (r) ∼ exp[− −2m(λ + E i )r/ ] for a paired system and ∼ exp[− √ −2mǫ i r/ ] for a unpaired system with ǫ i being the sp energy [52]. Thus, we have enhancement in the radius by blocking the Ω π = 1/2 − orbital though the chemical potential is not very shallow. This can be considered as unpaired-particle haloing [47].
To see the spatial structure of 37 Mg, we draw the calculated density distributions in terms of the equidensity lines on the ̺-z plane in Fig. 4(a). The contour lines are depicted in a logarithmic scale at 0.1 fm −3 down to 10 −7 fm −3 . In Fig. 4(b), shown are the density distributions of neutrons along the symmetry axis (at ̺ = 0.4 fm) also in a logarithmic scale. The density distributions obtained by blocking the Ω π = 5/2 − and 1/2 − orbitals are presented. In the case of Ω π = 5/2 − , the density distribution of neutrons is well localized in the center though the spatial extension is visible compared with that of protons, forming the neutron skin. Blocking the Ω π = 1/2 − orbital changes drastically the distribution of neutrons. A long tail emerges, interpreted as the neutron halo. The spatial distribution is extended toward the symmetry axis, forming a peanut shape. This is consistent with the previous calculation [47], and results from the p-wave dominance near the continuum threshold [53][54][55][56].

21
The time-reversal symmetry is intrinsically broken in the odd-mass isotopes possessing nonzero spin, so that the time-odd components in the mean field may be activated. Let us discuss finally the roles of the time-odd mean fields in 37 Mg by blocking the Ω π = 1/2 − orbital. The total binding energy is affected by 0.41 MeV, as shown in Fig. 2(a). Here the time-odd fields are set to zero, i.e. equivalent to the equal-filling approximation. This is only 0.16% to the total binding energy, and is negligibly small. A tiny effect on the nuclear mass has been brought out by the systematic calculation [12]. Accordingly, the deformation property is hardly influenced by the time-odd fields as shown in Fig. 2(b). The deformation parameter for neutrons is reduced by 0.01. The radius of neutrons thus calculated is lessened by 0.02 fm, while the protons are not affected. Then the matter radius is reduced by 0.02 fm as shown in Fig. 3. The chemical potential and the qp energy of the blocked orbital is −2.73 MeV and 1.92 MeV, respectively. Then, the qp energy of the last neutron is lowered by 0.28 MeV by ignoring the time-odd fields. It seems that this shift is also negligible. However, the asymptotic behavior of the halo structure is sensitive to the exponent of the qp wave function. Indeed, as shown in Fig. 4, the tail structure is influenced by the time-odd fields. Whether the time-odd mean fields enhance or reduce the halo structure depends on the EDF employed. The reaction observables sensitive to the outer surface of the halo nucleus can put a constraint on the time-odd coupling constants of the Skyrme EDF that are uncertain.

Summary
We have found relationships among the particle-number parity, the Bogoliubov transformation, and the time-reversal symmetry of the Hartree-Fock-Bogoliubov (HFB) Hamiltonian. Then we showed that the lowest-energy solution of the HFB equation is the particle-number parity being even as long as the time-reversal symmetry is conserved intrinsically. Based on this finding, we gave foundation of a method for solving the HFB equation to describe the ground state of odd-mass nuclei by employing an appropriate timereversal anti-symmetric constraint operator to the Hamiltonian. With this procedure, one can obtain directly the ground state of an odd-mass nucleus as a self-consistent solution of the cranked-HFB-type equation, while the ground state of an odd-mass nucleus is described as a one-quasiparticle excitation of a neighboring even-even nucleus in a usual procedure. This method is further applicable to the low-lying two-quasiparticle excitations in even-even nuclei. As a numerical example, we applied this method to the neutron-rich Mg isotopes close to the drip line, and showed that the anomalous increase in the matter radius of 37 Mg is well described when a neutron occupies the low-Ω orbital in the framework of the nuclear energy-density functional method. We found that the time-odd mean fields have little influence on the total binding energy, but an appreciable impact on the asymptotic behavior of the halo structure. and the particle-particle potential, or the pair Hamiltonian, is given bỹ h q σσ ′ (r) = δ σσ ′Ũ q (r), (B4) With the axial symmetry, the pairing density in the cylindrical coordinates r = (̺, φ, z) is given by employing the ansätze (70) as In the end, one obtains Γ even pair,t (r) = δ 0t U pair (̺, z) 0 0 U pair (̺, z) , Γ odd pair,t (r) = 0, andh q (r) = Ũ q (̺, z) 0 0Ũ q (̺, z) . (B8)