{'fraction_non_alphanumeric': 0.03517977615219583, 'fraction_numerical': 0.016543024788688708, 'mean_word_length': 4.5562257647565705, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three-dimensional Boltzmann-Hydro code for core-collapse in massive stars II. The Implementation of moving-mesh for neutron star kicks

We present a newly developed moving-mesh technique for the multi-dimensional Boltzmann-Hydro code for the simulation of core-collapse supernovae (CCSNe). What makes this technique different from others is the fact that it treats not only hydrodynamics but also neutrino transfer in the language of the 3+1 formalism of general relativity (GR), making use of the shift vector to specify the time evolution of the coordinate system. This means that the transport part of our code is essentially general relativistic although in this paper it is applied only to the moving curvilinear coordinates in the flat Minknowski spacetime, since the gravity part is still Newtonian. The numerical aspect of the implementation is also described in detail. Employing the axisymmetric two-dimensional version of the code, we conduct two test computations: oscillations and runaways of proto-neutron star (PNS). We show that our new method works fine, tracking the motions of PNS correctly. We believe that this is a major advancement toward the realistic simulation of CCSNe.


INTRODUCTION
There are a number of observational and theoretical indications that the inner engine of core collapse supernovae (CCSNe) is highly non-spherical. Recent observational developments include three-dimensional direct mapping of 44 Ti by NuStar, which reveals that the inner parts of the ejecta of SN 1987A and Cas A have experienced large-scale mixing and convection (Grefenstette et al. 2014;Boggs et al. 2015). This is consistent with the earlier evidence from polarimetric observations that the explosion is not spherical in general and becomes more so as ones sees deeper inside (see e.g., Wang & Wheeler (2008) and references therein). High spatial velocities of pulsars are suggested to result from the recoil in asymmetric explosions (Motch et al. 2009).
On the theoretical side, there have been various mechanisms proposed as the cause for these multi-dimensional features, which may be important also for the explosion mechanism itself. The stellar rotation may be the simplest. Unsteady accretions of turbulent matter due to strong convection in the last stage of stellar evolution could be the seed perturbations for the asymmetry of CCSNe, which would set off different kinds of hydrodynamical instabilities during the stalled-shock phase and enhance neutrino heating and increase turbulent energies in the post-shock flow . The nascent field of gravitational wave astronomy will be capable of directly investigating such asymmetric dynamics in the vicinity of proto-neutron star (PNS). The multidimensional modelling of CCSNe is hence indispensable to unveil the mechanism of CCSNe.
The neutrino radiation-hydrodynamic simulations of CCSNe have made a remarkable progress with ever increasing computational resources in the last few decades (for the current status, see e.g., Ott et al. (2008); Takiwaki et al. (2014); Kuroda et al. (2015); Dolence et al. (2015); Bruenn et al. (2014); Lentz et al. (2015); Müller et al. (2012); Just et al. (2015)). Although the present interest of supernova society is directed mainly to 3D hydrodynamical aspects, neutrino transfer is certainly one of the most important ingredients in realistic modelling of CCSNe. Since neutrinos are not in thermal equilibrium with matter in general, the time evolution of the neutrino distribution function at each spatial location should be determined in principle by solving the Boltzmann equation in six-dimensional phase space with both special and general relativistic effects taken into account properly. Both high numerical cost and technical difficulties have prevented us from conducting the ab-initio simulations, however. As a matter of fact, various approximations have been employed even in the most sophisticated multi-dimensional simulations, which include the multi-group flux-limited diffusion (MGFLD) (Mezzacappa et al. 1998;Zhang et al. 2013;Dolence et al. 2015), isotropic diffusion source approximation (IDSA) (Liebendörfer et al. 2009;Suwa et al. 2010;Takiwaki et al. 2014;Pan et al. 2016), two-moment method (Just et al. 2015;Skinner et al. 2015;O'Connor & Couch 2015), fastmulti-group transport (FMT) method (Müller & Janka 2015), variable Eddington tensor method (Müller et al. 2012) and 2D Boltzmann transport method without v/c corrections (Livne et al. 2004;Ott et al. 2008); some of them (Müller et al. 2012;Bruenn et al. 2014; Lentz et al. 2015) employ the ray-by-ray plus approximation further. It should be stressed that some of these approximations have been validated only under spherical symmetry and their performance in non-spherical situations, which no doubt prevail in the post-bounce supernova core, are still uncertain (see some also ; Skinner et al. (2015) for the validation of the ray-by-ray approximation). In fact, some recent simulations have yielded outcomes, that seem at odds with each other, the reason for which may be the different approximations they adopted for neutrino transfer (Müller et al. 2012;Takiwaki et al. 2014;Bruenn et al. 2014;Skinner et al. 2015). Multi-dimensional simulations with a Boltzmann solver, to which we refer in the following as Boltzmann-Hydro simulations, are hence indispensable to validate these approximations. They are obviously crucial to address the CCSNe mechanism.
Motivated by these facts, we have tackled the development of a Boltzmann-Hydro solver for the last few years. Sumiyoshi & Yamada (2012) constructed a three dimensional Newtonian Boltzmann solver, which was latter coupled to a hydrodynamics solver with self-gravity and, more importantly, was extended to fully accommodate special relativity . The latter paper demonstrated in 1D simulations of CCSNe the capability of our new code based on two different energy grids: Lagrangian-remapping grid and laboratory-fixed grid. Having been fine-tuned for massive parallel supercomputers, this code in the 2D version is currently being run on the K supercomputer for axisymmetric simulations of CCSNe. This paper is based on our solution to the problem we encountered in these computations.
The problem is the following: the nascent PNS starts to receive random kicks shortly after core bounce when the matter that has experienced prompt convection falls onto the PNS. It is further kicked around later by the matter that has undergone hydrodynamical instabilities such as the standing accretion shock instability (SASI) or neutrino-driven convection (see also Janka & Mueller (1994); Nordhaus et al. (2010); Scheck et al. (2006)). The PNS then moves at velocities of the order of 100km/s and is temporarily dislocated by several kilometers from the coordinate center (see also Nordhaus et al. (2012)). It is also important to remember that if the shock wave is successfully revived, asymmetric ejecta will certainly produce PNS kicks. The fixed polar coordinates are not appropriate to follow these proper motions of PNS, since dislocated spheres cannot be reproduced very well on this grid. Indeed, we found that the simulation without any special treatment either crashed or resulted in unphysical outcomes once the PNS moves a few km away from the mesh center.
The previous studies treated this problem rather pragmatically.
For example, Hanke et al. (2013);Takiwaki et al. (2014);Bruenn et al. (2014); Lentz et al. (2015); Müller (2015)) restricted the motion of PNS by artificially imposing spherical symmetry in the PNS. Although no numerical problem may have occured in these methods, it might have discarded some potentially important physics along. On the other hand, MPA group employed a moving-mesh technique, adding an extra velocity to the advection terms in the equations of motion that compensates the PNS motion (Scheck et al. 2006). In their experimental simulations Nordhaus et al. (2010) also used a similar moving-mesh method to track the PNS motion, remapping the grid. It is stressed, however, that the latter two groups completely ignored the effects of the moving-mesh on neutrino transfer.
In this paper, we propose an entirely new method to precisely treat this issue not only for hydrodynamics but also for neutrino transfer. The basic idea is like in the previous papers (Scheck et al. 2006;Nordhaus et al. 2010) to globally translate the coordinates so that the mass center of PNS should stay always very close to the mesh center. It should work, since the PNS remains almost spherical even when it oscillates or runs away 1 . The important thing is that the basic equations both on neutrino transfer and hydrodynamics should be modified on this moving grid, since it is not an inertial system. Hydrodynamics equations are easy to extend (see e.g., Scheck et al. (2006)) 2 , while the modification of Boltzmann equation is more complicated because the coordinate acceleration affects the neutrino transport in nontrivial ways in the six-dimensional phase space. Then the general relativistic description is a natural and probably the unique choice for handling various effects correctly. We deal with this within the 3+1 formulation of GR, using the conservation form of the general relativistic Boltzmann equation (Shibata et al. 2014). This means that the current upgrade of the neutrino transport module in our code is equivalent to a GR extension, which was actually planned as the next improvement to our Boltzmann solver. It should be mentioned, however, that in this paper the GR transport code is applied only to the flat Minknowski spacetime, since the treatment of gravity in our code is still Newtonian; the Newtonian version of the hydrodynamics code is employed in this work for the same reason (see also Section 4). The GR capability of the code will be demonstrated elsewhere (Nagakura et al. 2016).
This paper is organized as follows. In Section 2, we reformulate the Boltzmann equation on the moving-mesh with the language of the 3+1 formalism of GR. We also explain the numerical implementation of this GR extension to our previous special relativistic (SR) code in Section 3. Then, the feedbacks to hydrodynamics are described in Section 4. We validate our new method with two tests: PNS oscillations around and runaways from its original position. The results are presented in Section 5. Finally we conclude this paper with a summary in Section 6. Throughout this paper, Greek and Latin subscripts denote space-time and space components, respectively. We use the metric signature of − + ++. Unless otherwise stated, we work in units with c = G = 1, where c and G are the light speed and gravitational constant, respectively.

BASIC EQUATIONS
Following Shibata et al. (2014), we start with the conservation form of the Boltzmann equation in general rel- ativity: where g, x α are the determinant of the metric, coordinates of spacetime, respectively, and f is the neutrino distribution function; e α (µ) (µ = 0, 1, 2, 3) denote a set of the tetrad bases for a local orthonormal frame; ℓ i are directional cosines for the direction of neutrino propagation with respect to e α (i) (see Fig. 1 in Shibata et al. (2014)). The three components of ℓ i can be written as whereθ andφ stand for the polar and azimuthal angles (Lindquist 1966). We further define coordinates q i in momentum space-q 1 = ν, q 2 =θ and q 3 =φ with ν being the neutrino energy in this local orthonormal frame and also expressed as ν ≡ −p α e α (0) with the four momentum of neutrino, p α . In this paper neutrinos are assumed to be massless. ω (0) , ω (θ) , ω (φ) are given as As shown in Shibata et al. (2014), these ω's can be expressed with the Ricci rotation coefficients. S rad on the right hand side of Eq. (1) originates from the collision term for neutrino-matter interactions.
In the 3+1 formulation of GR, the line element is expressed as where α, β i and γ ij denote the lapse function, shift vector and spatial 3-metric, respectively. In our extended Boltzmann code, the time-like basis e α (0) is chosen so that it should coincide with the unit vector n α normal to the spatial hypersurface with t = const. This choice is a natural extension from our previous SR Boltzmann solver (see Section 3 for more details). Then three other spatial tetrad bases are taken so that they should be tangential to the spatial hypersurface. In this paper we assume that the spacetime is flat and is foliated with flat spatial hypersurfaces, on which we deploy the polar coordinates (x 1 = r, x 2 = θ, x 3 = φ). Then non-vanishing components of the 3 metric are γ rr = 1, γ θθ = r 2 and γ φφ = r 2 sinθ 2 . The spatial tetrad bases are chosen so that the e (1) be parallel to the radial coordinate, and e (2) be tangential to the surface spanned by ∂ t and ∂ θ , and e (3) be orthogonal to the other two: We refer to this orthonormal frame as the O-frame in the following. In accord with the above foliation of spacetime we set α = 1. We utilize the shift vector to deal with the motion of the spatial coordinates (see Figure 1). In fact, we set β i =V i , whereV i is approximately the velocity of PNS measured in the O-frame (see the next section for details). Note that the employment of the globally uniform shift vector in this paper should be compatible with the use of other gauge conditions for the shift vector in possible applications of the current formulation to (dynamical) curved spacetimes. This completes the description of Eq. (1). We now turn to its numerical implementations.

NUMERICAL IMPLEMENTATION
3.1. Shift vector Let us suppose that the basic equations are somehow finite-differenced and all hydrodynamics and space-time quantities are obtained up to the n-th time step. The average velocity of PNS at this time step (V i(n) ) is then given via the linear momentum (P) and mass of PNS (M) as where ρ, v i o and dV PNS denote the density, 3-velocity of matter (measured in the O-frame) and volume element in PNS, respectively. The PNS is defined to be the region, where the angle-averaged density (ρ) is larger than 10 13 g/cm 3 . The time derivative of the velocity, or the acceleration of PNS, at the same time step is given by the following relation: where ∆t (n) is the interval between the (n + 1)-th and n-th time steps.
Note that we do not use these V i and dV i /dt as they are for the following reasons. First, V i obtained in this way shows glitches from time to time when the PNS surface traverses an interface of the radial mesh points. Second, if the tracking of PNS motions were to be perfect, the acceleration of PNS should be determined iteratively, since the velocity of PNS at the next time step should be consistent with this acceleration but is obtained only after the advancement of the step. Such iterative process would be very time-consuming. Fortunately, however, it is unnecessary to exactly trace the motion of PNS and it turns out that the following approximate treatment suffices to deal with the proper motion of PNS.
We define the approximate PNS velocityV i as follows: with dV i(n− 1 2 ) /dt is given by Eq. (7) (but the backward difference); C (n) and D (n) are the terms that allow some deviations of the coordinate velocity and/or origin (denoted here by X i(n) m ) from those of PNS and thus avoid the glitch; T is the recovering time and is taken to be 0.1ms in this paper. C (n) and D (n) also prevent secular drifts of PNS. In fact C n works as a damper to prohibit large differences between two velocities, whereas D n serves as an attractor to ensure that the coordinate origin tends to the mass center of the PNS. As an additional measure to ensure smooth coordinate motions, we do not update the value of dV i(n) /dt when the PNS surface trespasses an interface of radial mesh points. As demonstrated later, we find that employingV as the shift vector in combination with the evaluation of dV i(n) /dt given above is indeed sufficient to solve the problems with the proper motion of PNS.
Although the shift vector field thus obtained is spatially uniform, its derivatives with respect to r, θ and φ are non-vanishing, since the coordinates are curvilinear. The explicit form of the Ricci rotation coefficients are rather involved (although calculations are straightforwardly) numerically in the code. This will be useful indeed, since we are required to evaluate Ricci rotation coefficients for numerically obtained metrices in truly GR simulations.

Modifications to SR code
Although the GR Boltzmann equation, Eq. (1), has a simple form, the consistent treatment of the advection and collision terms is complicated even for the flat spacetime. In Nagakura et al. (2014), we surmounted the difficulties by introducing two energy grids: Lagrangian remapped grid (LRG) and Laboratory fixed grid (LFG). We also devised for the SR code some other numerical techniques (e.g., a semi-implicit method for temporal sweep). It is therefore desirable in the GR extension to the current SR code that we should retain these features as much as possible. In the following, we describe how we achieved it.
We first consider the Boltzmann equation (1) in flat spacetime. The tetrad bases, Eqs. (5), are reduced in this case to where the subscripts "F" hereafter implies quantities in the flat spacetime. Then we can evaluate the ω variables in Eq.
(3) as: Substituting these results into Eq. (1) and using the determinant of the metric for the fixed polar coordinates in the flat spacetime ( −g (F ) = r 2 sinθ), we reproduce the SR Boltzmann equation we employed in Nagakura et al. (2014).
As we mentioned earlier, in the GR extension we want to retain the various features already implemented in our SR Boltzmann-Hydro code. This is most easily achieved by casting Eq. (1) into the following form: It should be apparent that these four variables in Eq. (13) can be regarded as the GR corrections to the SR equation. This allows us to directly utilize our SR Boltzmann code in the GR extension. Although we employ the moving spherical coordinates in the Minkowski spacetime in this paper, the GR-extended code can accommodate any metric and gauge conditions, evaluating various GR terms numerically. Note also that, unlike other advection terms, the energy-derivative term represents gravitational redshift, a purely GR effect, which we calculate on the LRG  3 . The treatment of the collision terms is also similar to that in the flat spacetime. Since the collision terms can be most easily calculated in the fluid-rest frame, we first evaluate them in this frame and then Lorentz-transform them to the laboratory frame, which is identical to the O-frame in the current formulation. This is done with the tetrads corresponding to these frames. We denote the tetrad bases of the fluid-rest frame asêμ, which is expressed with e (ν) aŝ where Λ stands for the Lorentz boost transformation. The components of Λ are given as where I 3 denotes the 3 × 3 identity matrix, v (i) and γ are defined with the tetrad bases of the laboratory frame and the fluid 4-velocity u as where η µν denotes the Minkowski metric. The 4momentum of neutrino is also projected onê (μ) . Then the energy shift and aberration are determined by the Doppler factor given above and our SR formulation can be directly passed over to the GR-extended code (see Section 4 in Nagakura et al. (2014)).

FEEDBACK TO MATTER
In this section, we describe the feedback from neutrino interactions to hydrodynamics in detail. We first present the fully GR formulation, and then take the Newtonian limit. It is noted that the hydrodynamics part of our code is also fully GR but its Newtonian version is employed in this paper, since the self-gravity part of the code is still Newtonian.
The basic equations for matter dynamics consist of the conservation laws of baryon number, energy-momentum and electron number, which are written as, respectively, where ρ 0 , T µν (hd) and N ν (e) denote the rest-mass density of baryons, energy-momentum tensor of matter and the 3 Other advection terms are calculated on LFG (see Nagakura et al. (2014) for more details). Note also that the energy derivative term disappears in the current study for the movingmesh in the flat spacetime, since ω (0) becomes trivially zero.
electron number 4-current, respectively. The right-hand sides of the latter two equations represent the feedback, which are related to the collision term of Boltzmann equation, Eqs (1) or (12), as follows: In these expressions, dV p (= νsinθdνdθdφ) denotes the invariant volume in the momentum space. The subscript "i" indicates the neutrino species, which we omit hereafter for simplicity.
In actual simulations, we first evaluate the tetrad components of G in the fluid-rest frame aŝ where the hat indicates variables in the fluid-rest frame, i.e.,Ĝ (μ) = G ·ê (μ) etc.. Then, the coordinate components of G can be expressed viaĜ (ν) andê (ν) as: whereê µ (ν) denotes the coordinate components ofê (ν) . In the 3+1 formalism the basic equations for matter dynamics can be expressed as where various variables are defined as follows: (see also Eq. A2 in Nagakura & Yamada (2008)). In the above equations, Y e , p, h, g µν γ, K ij are electron fraction, pressure and specific enthalpy of matter, 4-dimensional metric of spacetime, the determinant of 3-dimensional metric of space and extrinsic curvature, respectively. In this paper, instead of employing these fully GR equations, we adopt their Newtonian approximations, which can be derived by taking the weak gravitational field limit, ignoring the time derivative of gravitational potential and the space derivatives of 3-dimensional space metric (see Nagakura et al. (2011)). Then basic equations can be reduced the spherical coordinates to where each term is given as (see also Eqs.(12)-(16) in Nagakura et al. (2014)). In the above expressions, √ g(= r 2 sinθ), ψ andβ j denote the volume factor for the spherical coordinates, the Newtonian gravitational potential, and the time derivative of shift vector, respectively. Note that W a represents the acceleration of the coordinates.

Numerical setup and input physics
In this section, we validate our new formulation of the moving-mesh technique by performing 2D axisymmetric Boltzmann-Hydro simulations for the toy model of PNS oscillations and runaways from the original position. A thorough investigation of the code performance, including the GR capability, will be reported in the forthcoming paper (Nagakura et al. 2016). As the initial conditions of these tests, we employ a configuration of a supernova core at 100ms after core bounce, which is obtained by a 1D spherically symmetric simulation of core-collapse, bounce and stagnation of shock wave with the same code. We map the resultant 1D configuration onto the 2D grid.
We use a 11.2M ⊙ progenitor (Woosley et al. 2002). For the 1D simulation, a non-uniform grid of N r = 384 points in the radial direction covers the region of 0 ≤ r ≤ 5000km, while the momentum space is discritized with a grid of N ν = 20 points in the energy region of 0 ≤ ν ≤ 300 MeV and Nθ = 10 points covering the polar angle from 0 to π. For the 2D simulations, on the other hand, we deploy N r = 192 radial grid points and uniformly distributed N θ = 32 angular grid points in the entire meridian section; for the momentum space we use N ν = 20 energy grid points and Nθ(= 4) × Nφ(= 4) angular grid points. Although this is admittedly a rather coarse mesh both in space and momentum space, it is not a serious issue for the current purpose: a proof-of-principle kind of code tests. We take into account three species of neutrinos: ν e , νē, and ν x , which are electron-type neutrinos, electron-type antineutrinos and heavy-lepton neutrinos (µ and τ neutrinos and their anti-particles collectively), respectively.
As for the input physics, our current Boltzmann-Hydro code has updated some treatments of microphysics from Sumiyoshi & Yamada (2012); Nagakura et al. (2014). One of them is an incorporation of the multi-nuclear species equation of state (EOS) by Furusawa et al. (2011Furusawa et al. ( , 2013a. This tabulated EOS provides us with not only  Figure 2 but for the post-bounce phase. Note that the horizontal axis is not the mass coordinate but the radius in these panels. thermodynamics quantities but also abundances of nuclei with mass numbers up to A ∼ 1000 in nuclear statistical equilibrium or NSE, which are then employed to obtain the rate of electron captures by heavy nuclei (see below). Incidentally the EOS includes also the information on the abundances of light elements. We are currently studying possible roles of their interactions with neutrinos in the post-bounce phase of CCSNe.
Neutrino-matter interactions have been also improved from Sumiyoshi & Yamada (2012). One of the upgrades is the full account of non-isoenergetic scatterings between neutrinos and electrons and positrons. Unlike the spherically symmetric case, the interaction rate depends onφ, the fact which prevents a direct application of the method used in Mezzacappa & Bruenn (1993); Sumiyoshi et al. (2005). We hence obtain the interaction rate by direct numerical integrations combined with the Chebychev expansions of Polylogarithms (Kolbig et al. 1970). It is important to note that an implicit treatment of non-isoenergetic scatterings is highly expensive both in computational time and memory. We hence handle the neutrino scatterings on electrons and positrons explicitly. Since they are rather minor, having smaller reaction rates, compared with other reactions such as emissions, absorptions and scatterings on nucleons, this poses no problem.
Another important upgrade is the treatment of electron captures by nuclei as mentioned earlier. We tabulate the reaction rates based on the results by Juodagalvis et al. (2010) and the approximation formula of Langanke & Martínez-Pinedo (2000) and Langanke et al. (2003) with the mass fractions of heavy nuclei being taken from the Furusawa's EOS.
As shown in Lentz et al. (2012), these two updates are critically important for CCSNe. As a matter of fact, the deleptonization would be erroneously suppressed during the infall phase if they were neglected, which would then result in a larger mass of the inner core and a stronger shock wave (see below and also Mezzacappa & Bruenn (1993); Langanke et al. (2003); Hix et al. (2003); Lentz et al. (2012)).
Figures 2-6 display the result of the 1D simulation. Figure 2 plots the distributions of density, radial velocity, electron fraction, lepton fraction, entropy per baryon and temperature at different times in the pre-bounce phase. We find that the mass of the inner core is somewhat less than 0.6 M ⊙ due to the unsuppressed deleptonization, which is consistent with other 1D computations (see, e.g., Lentz  bounce counter parts are shown in Figure 3 for different times, where T b denotes the time after bounce in this figure. Figure 4 presents the energy fluxes measured in the laboratory frame (i.e., neutrino luminosities) at r = 422 km. Note that the well-known bounce feature, i.e., a slight decrease followed by a quick rise in the luminosity of electron-type neutrino reaches this radius at T b ∼ 4ms. The prominent neutronization burst of electron-type neutrinos can be clear seen in the upper panel, while the luminosities of other neutrinos start to rise somewhat later. Note that the production of electron-type anti-neutrinos is initially suppressed by high electron fractions around the neutrino sphere (Y e ∼ 0.3). In fact, the luminosity of heavy-lepton neu-  Sumiyoshi et al. (2005) in 1D simulations for the 15M ⊙ progenitor model by Woosley et al. (2002). The radial velocity, entropy, temperature and electron fraction at core bounce are shown clockwise from the upper left panel. The red line corresponds to the results of the current code, while the green line represent ones for the other code. trino goes up in ∼ 10ms, while that of electron-type anti-neutrinos increases gradually for ∼ 50ms. Figure 5 shows the time evolution of the mean energy for each neutrino species at r = 422 km. The mean energy is defined as We also show the trajectories of some important radii for the post-bounce phase in Figure 6. The shock expands initially but is stagnated around 170km at T b ∼ 80ms, while the gain radius starts to deviate from the shock radius at T b ∼ 20ms. The trajectories of the points that have the densities of ρ = 10 11 , 10 12 and 10 13 g/cm 3 are also shown in this figure. They will serve as a rough guide to the size of PNS as a function of time. All of these features are in qualitative agreement with those observed in previous studies (see e.g., Lentz et al. (2012); Müller & Janka (2014)). The data at T b = 100ms are used as the initial condition for the subsequent 2D simulations.
As a quick validation of the core part of our code, we show in Figure 7 some results of the comparison with another code, in which we ran 1D simulations twice for the 15M ⊙ progenitor model by Woosley et al. (2002), using the current code and another 1D but fully GR Boltzmann-Hydro code developed by (Sumiyoshi et al. 2005). All input physics are identical between the two computations but the latter is fully general relativistic and employs the Lagrangian formulation. In the figure, we show some key quantities at core bounce. In each panel, the red (green) line gives the result of our new (Sumiyoshi et al. (2005)) code. Considering the differences just mentioned, we think the two results agree with each other reasonably well.

PNS Oscillation
It is known that the center of PNS oscillates with velocities of the order of 100 km/s and periods of several tens of milliseconds (see e.g., Nordhaus et al. (2012)). Mimicking such a situation, we start a 2D simulation by adding a velocity in the z-direction as ∆v z = 100km/s in the region of r < 30km. This simulation is carried out for 100ms, which is long enough for the purpose of this study. Figure 8 shows the time evolution of PNS velocity (upper panel), trajectory of the origin of moving-mesh (middle panel), and the z-coordinate of the mass center of PNS on this moving-mesh (lower panel). One can see from the middle panel that the PNS moves in the positive z direction initially owing to the velocity added at the beginning. It is also clear from the upper panel that the PNS is decelerated and the direction of motion is reversed after a few ms. It is noted the PNS moves from the original position by ∼ 1km by the time of t ∼ 15ms as shown in the middle panel. This is not a small distance and the subsequent evolution can not be calculated without the moving-mesh technique. As a matter of fact, we conducted the same simulation on the ordinary nonmoving grid, and found that it ended up with a numerical crash with unphysical matter distributions around the coordinate origin.
At t ∼ 20ms the PNS again changes the direction of motion and returns to the origin. Although the period is variable in time, the PNS experiences two cycles of oscillations in this simulation. It is also important to note that the moving-mesh nicely traces the motion of PNS (see the lower panel). This leads to the successful Boltzmann-Hydro simulation on the spherical polar grid. We confirmed indeed that all the unphysical features observed in the simulation on the fixed mesh have gone with the moving-mesh technique.
5.3. Runaway motions of PNS In realistic simulations the shock wave expands asymmetrically after successful shock revival and expels the envelope anisotropically. The PNS then undergoes a recoil (Janka & Mueller 1994;Nordhaus et al. 2012). In this section, we mimic such runaways of PNS very crudely. The initial configurations of matter and neutrinos are the same as those used in the previous section. In this case, however, we do not add velocity perturbations. Instead, we continuously add by hand the external acceleration of 10 11 cm/s 2 in the positive z-direction within the region of r < 30km on the moving grid. The simulation was carried out for t = 10ms. Figure 9 is the counter part of Figure 8 for the present case. As expected, the PNS moves continuously in the positive z-direction. It is seen in the top panel that the PNS velocity reaches ∼ 3000km/s at the end of the simulation, which is much larger than the realistic kick velocity of a few hundred km/s. In spite of this rather extreme runaway of PNS, the moving-mesh tracks it very well, as shown in the lower panel of the figure. In fact, the distance between the mass center of PNS and the origin of the moving-mesh remains less than 10 −2 km, which is close enough to avoid numerical problems. Incidentally, oscillations that are evident after t ∼ 4ms in the lower panel are ascribed to the deformations of PNS.
A series of snapshots of the entropy distribution in the meridian plane are shown in Figure 10. The low entropy region with deep blue colors in this figure corresponds to the central un-shocked part of PNS. It is apparent that the PNS moves upward with time. Note that this figure is drawn on the fixed coordinates, which coincide with the moving coordinates initially. In each panel of the same figure, we put two concentric circles, which represent the moving-mesh. It is confirmed again that they trace the PNS closely. In Figure 11, on the other hand, the (number) density contour is drawn in color for electron-type neutrinos. One can see that the neutrino density is slightly non-spherical owing to the deformation of PNS. It is more important, however, that the neutrinos are comoving with PNS. This is of course a consequence of the neutrino trapping, which occurs in the optically thick region. In the Boltzmann-Hydro simulation, however, it is highly non-trivial and is in fact ensured by the combination of the following two conditions: (1) neutrinos are isotropically distributed in the fluid-rest frame; (2) the neutrino distribution in the O-frame is related accurately with that in the fluid-rest frame by a Lorentz transformation. The result we have just presented is yet another demonstration that our code is working properly.
Last but not least, we mention the conservation of linear momentum in our code. It is a well known fact that it is difficult in general for hydrodynamics codes like ours that adopt curvilinear coordinates to enforce the conservation of linear momentum. It is evident indeed in Eqs (31) or (42) that these equations can not be written in the conserved form even in the absence of gravity and neutrino interactions. Not to mention, the gravity term written not in the conservation form also attributes to the violation of momentum (and energy) conservation (see also Müller et al. (2010)). Note, however, that as formulated in Section 3.1, we do not use the conservation law to evaluate the PNS velocity and its acceleration. Regardless, we checked quantitatively the violation of linear momentum in our code by conducting another 2D purely hydrodynamical simulation (same initial conditions as previous tests but adding 1% random density perturbation) for 100ms. We found that the numerical error is equivalent to ∼ 10km/s of the PNS's kick velocity, which is not negligibly small but still much smaller than the typical velocity of ∼ 100km/s. Considering the rather coarse grid employed in this paper and the purpose of this paper, we may conclude that our code performed well.

SUMMARY
In this paper, we have presented a novel method to deal with motions of PNS on spherical polar coordinates. It is based on a moving-mesh technique, as far as the neutrino transport part is concerined, it is essentially equivalent to the general relativistic extension of the special relativistic Boltzmann solver we developed earlier. In fact, the Boltzmann equation is reformulated in the 3+1 formalism of GR although the GR code thus obtained is applied only to the flat spacetime coupled with the Newtonian hydrodynamics code and self-gravity module in this paper. The shift vector is utilized to specify the movement of spatial coordinates so that they could track the PNS motion approximately. As a matter of fact, without such a technique we encountered a numerical crash with unphysical features emerging at the coordinate origin, as the PNS is dislocated from the original position. Since the coordinate origin stays very close to the mass center of PNS with the moving-mesh technique, we expect that in more realistic simulations the extended Boltzmann-Hydro code will be able to treat the violent oscillations and ultimate runaway of PNS in the post-bounce phases of CCSNe on spherical polar coordinates.
We have also described in detail the numerical implementations of the GR extensions to our SR Boltzmann code, which was constructed on two energy grids so that it could deal with both the advection and collision terms correctly. It turns out that the extension is rather straightforward thanks to the use of appropriate tetrads. The two energy grids and the transformations between them, which are employed in the SR Boltzmann code, are nicely identified with these tetrads and their Lorentz transformations, respectively.
In Section 5, we have validated our method by applying it to two test problems: toy models of PNS oscillations and a runaway, which mimic very crudely more realistic post-bounce simulations. We have demonstrated that the code has nice tracking capabilities and can follow the evolutions without any problems such as those we encountered when the fixed spatial grid was deployed. Incidentally, the extended code is currently being applied to realistic 2D simulations of CCSNe and the results and a thorough code validation will be reported elsewhere soon. We are also planning to conduct truly GR Boltzmann simulations of neutrino transport in a black hole spacetime. It should be apparent that any (local) gauge conditions can be imposed in addition to the uniform shift we employed in this paper. The next step is to couple the GR neutrino transport code with a solver of the Einstein equations. Note that the hydrodynamics code is already GR (see Nagakura & Yamada (2008)) although the Newtonian version was used in this paper. Such an integrated code, once completed, will certainly broaden the scope of application much beyond CCSNe.
We are grateful to A. Juodagalvis for providing the data of electron capture rates on heavy nuclei. H.N. acknowledges to M. Shibata, Y. Sekiguchi and H. Okawa for valuable comments and discussions. H.N. also thanks Werner Marcus for proofreadings. The numerical computations were performed on the supercomputers at K, at AICS, FX10 at Information Technology Center of