FULL SIMULATION OF THE OPEN FIELD LINES ABOVE A PULSAR'S POLAR CAP. I. ACCELERATION

After 40 years of pulsar research, there are open questions, and still no consistent quantitative model. Models of pulsar activity are associated with an electron–positron pair production cascade, which requires particles' acceleration to a high Lorentz factor. Most polar cap cascade models assume a steady state—a stationary particle outflow. Due to the pair cascade, the magnetosphere is filled with electron–positron plasma, the accelerating electric field is screened, and the magnetosphere is assumed to be a force-free magnetosphere. Numerical models of the force-free pulsar magnetosphere presented by Contopoulos et al. (1999), Gruzinov (2005), Spitkovsky (2006), McKinney (2006), Komissarov (2006), Timokhin (2006), Kalapotharakos & Contopoulos (2009), and Bai & Spitkovsky (2010) indicate that the required particle accelerations needed to produce the observed radiation cannot occur when the force-free and/or ideal MHD condition E  ·  B = 0 is valid in the entire magnetosphere, and that a very small spatial region of activity is required. A good candidate for those physical conditions is the polar cap. The current density in these global force-free magnetosphere models is not a free parameter, but rather dictated by the global structure of the magnetosphere. The current density of the force-free magnetosphere models is incongruent with the one-directional Goldreich–Julian current density of the steady-state models (the current density adjustment problem). Here, we take into account time variability and, instead of making assumptions about and approximations for pair injection, we will calculate the pair production directly from the accelerating particles and let the simulation show how the current and charge densities in the open magnetic field lines are adjusted to maintain the global structure of the magnetosphere.

We consider one of the basic models for pulsar activity: an electron's beam extracted freely from the star surface above the polar cap region (electron polar cap Ω  ·  B > 0), using time-dependent, particle-in-cell (PIC) simulation, as presented by Birdsall & Langdon (1991). A complementary time-dependent PIC simulation, the Ruderman–Sutherland model, in which particles cannot be extracted from the star surface, is presented in Timokhin (2006, 2010).

In this paper, the electrons are extracted freely from the star surface and accelerate along the open field lines, due to the electric field induced by the rotation of the magnetized pulsar, as space charge limited flow (SCLF), as described by Michel (1974), Fawley et al. (1977), Cheng & Ruderman (1977), and Arons & Scharlemann (1979) (hereafter M74, FAS77, CR77, and AS79, respectively). Frame dragging, a general relativity (GR) effect, enables the acceleration of the electrons to much higher energies, as found by Muslimov & Tsygan (1992) and Muslimov & Harding (1997) (hereafter MT92 and MH97, respectively). Those primary relativistic electrons initiate the cascade of pair production. The initial calculations only consider pair production via curvature radiation (CR; AS79). Now, inverse Compton scattering (ICS) and synchrotron radiation (SR) are also considered important mechanisms for both particle energy loss and pair production (Dermer 1990; Sturner 1995; Zhang & Qiao 1996, 1997; Harding & Muslimov 1998; Hibschman & Arons 2001). Pair multiplicity (the number of pairs produced by a primary particle) from the ICS cascade is lower than that expected from the CR cascade by two orders of magnitude (Hibschman & Arons 2001).

We report, for the first time, the SCLF acceleration along the open magnetic field lines, using PIC numerical code, first without the GR frame-dragging correction and then with it. The pair cascade results, including CR, ICS, and SR, using the cross-sections from the literature, the Monte Carlo code, and the photon attenuation length (the distance a photon travels before pair production occurs), will be reported in a different paper. The cascade consequence, such as the positron return current (thus, we can calculate polar cap heating and X-ray emission), the photons' and electrons' observed spectra (photons and electrons that escape the magnetosphere after the cascade), the pulsar death line (pulsars without enough pair production), and the PFF height (pair formation front location), will be reported in a future paper.

The model and the equations are introduced in Section 2, numerical results are presented in Section 3, and conclusions appear in Section 4.

We consider a cylindrical conductive tube with constant radius R (which is the polar cap radius), with negative background Goldreich–Julian charge density ρ_GJ. Here, a brief discussion of the assumption of a constant tube radius is needed: the tube can be assumed to have a constant radius to a height of the order of the polar cap radius R_p. The simulation follows the particles (and photons) to a height of the order of 100 R_p. The assumption of a constant tube radius is still valid for particle acceleration, since the acceleration depends only on the difference ρ − ρ_GJ, and both are influenced in the same way from field line flaring. However, for pair production, the field line curvature cannot be ignored. To calculate the photon attenuation length (the distance the photon travels before pair production occurs), the angle between the photon direction and the magnetic field should be taken into account. A photon (emitted either by CR, ICS, or SR radiation) can produce a pair in a strong magnetic field via γB pair production, provided its energy exceeds the threshold energy ε_ph ⩾ 2mc²/sin ψ, where ε_ph is the photon energy and ψ is the angle between the photon direction and the magnetic field line. Since the photons are emitted roughly parallel to the magnetic field lines from accelerating particles, they are well below the threshold. The photons have to travel until their angle, ψ, with an adjacent curved field line is large enough to allow pair production. In the simulation, the radius of curvature at the photon location is calculated using the distance from the photon to the star. The angle ψ is then calculated from the radius of curvature and the initial angle ψ_i where the photon was emitted.

In every time step, electrons are continuously inserted into the base of the tube, as long as ρ_GJ at the base of the tube is not compensated. These are the extracted electrons from the star surface; positive ions are assumed to be bound. The tube is represented in the simulation as discrete cells on a grid. The fields are discrete at each cell, while the particle positions and velocity are continuous values, and weighting (interpolation) between the discrete and continuous values is done. The fields are calculated from Maxwell's equation by knowing the position and velocities of all particles. The forces on the particles are found using the electric field, and the positions and velocities of the particles are found using the equations of motion. This is repeated at each time step. Figure 1 shows one time step in the simulation (not including pair production).

Zoom In Zoom Out Reset image size

The equation of motion is:

$$\begin{equation} F = qE_z = \frac{{dp}}{{dt}} = \frac{{d(\gamma mv)}}{{dt}},\quad v = \frac{dz}{dt}, \end{equation} \tag{ 1 }$$

where γ is the Lorentz factor, z is the height above the surface, and q is the electron charge.

Maxwell's equations (FAS77) are:

$$\begin{equation} \begin{array}{l} \nabla \cdot E = 4\pi (\rho - \hat \rho) \\\\[-6pt] \nabla \cdot B = 0 \\\\[-6pt] \nabla \times B = \dsty\frac{{4\pi }}{c}(j - \hat j) + \frac{1}{c}\dsty\frac{{\partial E}}{{\partial t}} \\\\[-6pt] \nabla \times E = - \dsty\frac{1}{c}\frac{{\partial B}}{{\partial t}}. \\ \end{array} \end{equation} \tag{ 2 }$$

Well inside the light cylinder, $\hat \rho = \rho _{{\rm GJ}},\hat j = 0$.

Note that j ≠ 0 and evolves according to the particles entering or exiting the cell.

We rewrite those equations in polar coordinates, with cylindrical symmetry (∂/∂φ = 0), conductive walls (E_z(r = R) = 0), and transverse derivatives ∂/∂r → 1/R:

$$\begin{equation} \begin{array}{l} \dsty\frac{1}{c}\dsty\frac{{\partial E_z }}{{\partial t}} = \dsty\frac{{B_\varphi }}{R} - \dsty\frac{{4\pi }}{c}j \\\\[-6pt] \dsty\frac{1}{c}\dsty\frac{{\partial E_r }}{{\partial t}} = - \dsty\frac{{\partial B_\varphi }}{{\partial z}} \\\\[-6pt] \dsty\frac{1}{c}\dsty\frac{{\partial B_\varphi }}{{\partial t}} = - \dsty\frac{{\partial E_r }}{{\partial z}} - \dsty\frac{{E_z }}{R}. \\ \end{array} \end{equation} \tag{ 3a }$$

The two boundary conditions are B_φ(z = L) = 0 and E_r (z = 0) = 0, where L is the length of the tube.

We write those equations via the leap-frog method (Lyubarsky 2009):

$$\begin{eqnarray} &&\frac{1}{c}\left[ {B_\varphi \left(z,t + \frac{{\Delta t}}{2}\right) - B_\varphi \left(z,t - \frac{{\Delta t}}{2}\right)} \right]\nonumber\\ &&\quad = - \frac{{\Delta t}}{{\Delta z}}\left[ {E_r \left(z + \frac{{\Delta z}}{2},t\right) - E_r \left(z - \frac{{\Delta z}}{2},t\right)} \right] - \frac{{\Delta t}}{R}E_z (z,t)\nonumber \\ &&\frac{1}{c}\left[ {E_z (z,t + \Delta t) - E_z (z,t)} \right]\nonumber\\ &&\quad = \frac{{\Delta t}}{R}B_\varphi \left(z,t + \frac{{\Delta t}}{2}\right) - \frac{{4\pi }}{c}\Delta t_{_{} } J\left(z,t + \frac{{\Delta t}}{2}\right)\nonumber \\ &&\frac{1}{c}\left[ {E_r \left(z + \frac{{\Delta z}}{2},t + \Delta t\right) - E_r \left(z + \frac{{\Delta z}}{2},t\right)} \right]\nonumber\\ &&\quad = - \frac{{\Delta t}}{{\Delta z}}\left[ {B_\varphi \left(z + \Delta z,t + \frac{{\Delta t}}{2}\right) - B_\varphi \left(z,t + \frac{{\Delta t}}{2}\right)} \right]. \end{eqnarray} \tag{ 3b }$$

We initiate the electric fields by solving the Poisson equation (FAS77) with an empty tube (ρ = 0):

$$\begin{equation} \nabla ^2 \phi = \nabla _{||}^2 + \nabla _ \bot ^2 = \frac{{d^2 \phi }}{{dz^2 }} - \frac{\phi }{{R^2 }} = - 4\pi ({\rho - \rho _{{\rm GJ}} }), \end{equation} \tag{ 4 }$$

where Φ is the electric potential. Since the tube is initially empty, the magnetic field B_φ (t = 0) = 0.

The boundary conditions for the Poisson equation are Φ = 0 at the bottom of the tube (the star surface) and Φ' = 0 at the top of the tube. The additional boundary condition in the tube Φ = 0 on the tube walls is satisfied automatically by the way we rewrite Maxwell's equations: we assume conductive walls. This potential is the initial potential of a vacuum tube. After initiating the fields, the simulation evolves through the feedback between the fields and the particles' positions and velocities.

We get the SCLF condition Φ' = 0 at the bottom of the grid (the star surface) as the simulation runs. We also get a cloud of electrons at the bottom of the grid, as expected in the SCLF model.

The background Goldreich–Julian charge density ρ_GJ should be initiated, too. For flat spacetime, ρ_GJ is presented in AS79:

$$\begin{equation} \rho _{{\rm GJ}} = - \frac{{\Omega _* B_* }}{{2\pi _{} c}}\left({\frac{{R_* }}{r}} \right)^3 \left[ {\cos i + \frac{3}{2}\left({\frac{r}{{R_* }}} \right)^{1/2} \theta _* \sin i\cos \varphi } \right], \end{equation} \tag{ 5 }$$

where i is the inclination angle (the angle between magnetic and rotation axes), r = R_*+z, R_* is the star radius, and φ and θ_* are the magnetic colatitudes and azimuth angle of a field line where it intercepts the stellar surface (φ = 0, θ_* = ξ θ_p = 0.67 θ_p, chosen to be the field line where the plasma flow is centered; AS79). θ_p is the polar cap angle, Ω_* is the star's angular velocity, B_* is the magnetic field at the star surface, and c is the speed of light. The second term of Equation (5) becomes important far from the star, and is vanishing for an aligned rotator. We dropped the factor (R_*/r)³, since in our model the tube has constant radius R.

In the simulation, the above ρ_GJ is inserted into the grid, and then the Poisson equation is solved to initiate the fields. The simulation runs were described as before: at each time step, all the particles' locations and velocities are updated, as are the fields for the entire grid. The simulation's results (the blue line in Figures 2–4) are compared to the analytic expression given in M74 (the red line in Figures 2–4).

Zoom In Zoom Out Reset image size

MT92 and MH97 have shown that GR frame dragging is very important for particles acceleration. They derive ρ_GJ in curved spacetime (which is reduced to AS79's expression in flat spacetime):

$$\begin{eqnarray} \rho _{{\rm GJ}} &=& {-} \frac{{\Omega _* B_* }}{{2\pi _{} c\alpha }}\left({\frac{{R_* }}{r}} \right)^3 \frac{{f(\eta)}}{{f(1)}}\nonumber\\ &&{\times}\left[\!{\left({1\,{-}\,\frac{\omega }{{\Omega _* }}} \right)\cos i\,{+}\,\frac{3}{2}\frac{{H(\eta)}}{{f(\eta)^{1/2} }}\left(\!{\frac{r}{{R_* }}}\!\right)^{1/2} \theta _* \sin i_{} \cos \varphi } \right].\nonumber\\ \end{eqnarray} \tag{ 6 }$$

The GR information is stored in the two variables α and ω. α = (1 − r_g/r)^1/2 is the gravitational redshift, r_g = 2 GM_*/c² is the gravitational radius (r_g/R_* ∼ 0.4), and ω is the modification to the star's angular velocity due to the frame-dragging effect. Here, we modify ρ_GJ by replacing Ω_* with Ω_* −ω:

$$\begin{equation} \omega = \Omega _* \frac{2}{5}\left({\frac{{r_g }}{{R_* }}} \right)\left({\frac{{R_* }}{r}} \right)^3 = \Omega _* \frac{\kappa }{{\eta ^3 }}, \end{equation} \tag{ 7 }$$

where η = r/R_*. The functions f(η) and H(η) are given in MT92 and MH97, and are numerical factors for curved spacetime that are equal to 1 in flat spacetime. The sin i term is comparable in flat and curved spacetime, and much smaller than the cos i term for the relevant range of acceleration. Since the cosi term is significant at the relevant distance from the star, pulsars with small inclination angles accelerate particles to higher energies. Full GR treatment includes modifying Maxwell's equation and the equation of motion with the factor α as well, as described in MT92 and MH97, and additional effects described in Gonthier & Harding (1994; estimated to result in ∼10% corrections). As before, we drop the field line curvature effect in ρ_GJ.

In the simulation, the expression of Equation (6) is inserted into the grid instead of Equation (5), and then the Poisson equation is solved to initiate the fields and the simulation is run as described above. The simulation's result (the blue line in Figures 5 and 6) is compared to the following analytic expression (the green line in Figures 5 and 6):

$$\begin{equation} \gamma _{\max } = \frac{{e\;\Phi }}{{mc^2 }} + 1, \end{equation} \tag{ 8 }$$

where Φ is given in MT92 and MH97 (the field-line curvature effect in Φ should be removed).

The result in this paper is for a "standard" pulsar of 1 s period, magnetic field of B = 10¹² G, star radius R_* = 10⁶ cm, polar cap angle θ_p = $\sqrt(\Omega_{*}$ R_*/c) = 0.014, polar cap radius R_p = R_* θ_p ∼ 10⁴ cm, and star mass 1.4M_sun. The Goldreich–Julian charge density at the star surface is n_GJ = 6.6 × 10¹⁰ × cos i cm⁻³ (ρ_GJ = n_GJq = Ω_*Bcos i/2πc). The inclination angle i was taken to be i = 0 for the runs of flat spacetime SCLF, and i = 0.5 rad ∼28° for the GR frame-dragging correction. All those values are parameters of the simulation and can be changed for each run. The simulation values n_GJ (charge number density), q (electron charge), and m (electron mass) are normalized according to the plasma frequency $\omega _p = q\sqrt {4\pi _{} n_{{\rm GJ}} /m}$. Distances are normalized according to the polar cap/tube radius R. The speed of light in the simulation is taken to be c = 1.

The first step of the simulation is a run of flat spacetime SCLF. This step provides a relation between "real-world" values and simulation values. It is performed to "calibrate" the simulation values with the pulsar values, i.e., to indicate the values of (R, ω_p) that should be taken in the simulation, to represent the pulsar values.

The simulation's results for this step are presented in Figures 2–4; the axes are the Lorentz factor γ and the normalized distance from the star z/R. The simulation's result (the blue line in Figures 2–4) is compared to the analytic expression given in M74 (the red line in Figures 2–4). The simulation's result shows a "front." The initial vacuum potential found by solving the Poisson equation (Equation (4), which was used to initiate the electric field) is dropped significantly after the particles enter the tube. There is a "front" that ahead of it, where the particles have not yet reached and the potential (calculated from the electric field) and fields have the same values as the initial values. Behind the "front," where the particles are in the tube, the potential and fields are altered to the SCLF values, and the particles have the maximum Lorentz factor expected by SCLF, with an oscillatory pattern. In the "front," there are strong oscillations of the fields and the particle's Lorentz factor.

Table 1 and Figure 2 show results of different tube radii R and different plasma frequencies ω_p, such that the maximum Lorentz factor γ_max remains constant. The last two columns of Table 1 represent the "real world." The results of Figure 2 are all taken when the particles reach 30R. Moving from the left panel to the right panel (panels (a)–(c)), R increases (3000, 5000, and 10,000, respectively); thus, it takes more time steps "t" (90,000, 150,000, and 300,000, respectively) and more particles (78,161, 130,237, and 260,263, respectively), while the pattern looks the same. Here, the negative Lorentz factor indicates returning electrons, i.e., electrons that move backward.

Table 1. Analytic Results for Flat Spacetime SCLF with the Same γ_max

Parameter	Simulation	Simulation	Simulation	Simulation	Real World
				Fit for "Real World"
	R = 3000	R = 5000	R = 10,000	R = 3000	Rp ∼ 104 cm
$\gamma _{\max } = \sqrt 2 \omega _p R/c$	∼15	∼15	∼15	∼7000	∼7000
$\omega _p = q\sqrt {4\pi n_{GJ} /m}$	0.0035	0.0021	0.0011	1.66	1.5×1010 s−1

Download table as:  ASCIITypeset image

Table 2 and Figure 3 show results of different tube radii R, while ω_p is constant. The last two columns in Table 2 represent the "real world." The results of Figure 3 are taken when the particles reach 30R. Here again, R increases from the left panels to the right panels, which results in more time steps and more particles involved. As R increases while ω_p is constant (approaching "real-world" values), γ_max increases, too, and the shape becomes more oscillatory.

Zoom In Zoom Out Reset image size

Table 2. Analytic Results for Flat Spacetime SCLF with Varying R and Equal ω_p

Parameter	Simulation	Simulation	Simulation	Simulation	Real World
	R = 3000	R = 5000	R = 10,000	R ∼ 5.0× 106	Rp ∼ 104 cm
$\gamma _{\max } = \sqrt 2 _{} w_p R/c$	∼5	∼8	∼15	∼7000	∼7000
$\omega _p = q\sqrt {4\pi _{} n_{{\rm GJ}} /m}$	0.0011	0. 0011	0. 0011	0. 0011	1.5×1010 s−1

Download table as:  ASCIITypeset image

Figure 4 shows the results of the same (R, ω_p), with different distances from the bottom of the tube (the star surface). The results are shown for 30R, 50R, and 100R. Although it takes more time steps and more particles, the pattern looks the same. The bottom panels show the corresponding fields. The red line is the parallel accelerating electric field E_z, the blue line is E_r, and the green line is B_φ. Note the SCLF condition that appears at the bottom of the tube, i.e., the vanishing of the parallel electric field at the star surface.

Zoom In Zoom Out Reset image size

The second step in the simulation is run with the GR frame-dragging correction. The simulation results for this step are presented in Figures 5 and 6; the axes are the Lorentz factor γ(×10⁶) and the normalized distance from the star z/R. The simulation's result (the blue line in Figures 5 and 6) is compared to the analytic expression given in MT92 and MH97 (the green line in Figures 5 and 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the results of the "standard pulsar" parameters described above. Two different possibilities of the pair (R, ω_p) that fit the pulsar parameters are shown (6,253, 1.1) and (10,317, 0.66), respectively, using the SCLF calibration. The results are taken at two different distances from the bottom of the tube: ∼43R (upper panels) and ∼53R (lower panels). Since R increases from the left panels to the right panels, more time steps and particles are involved in the right panels, but the figures look the same (the two upper panels and the two lower panels). As can be seen, there is a "front." Behind the front, the particles have the Lorentz factor expected from the GR frame-dragging correction. These figures show an oscillatory pattern, with increasing oscillations as the distance from the star increases. The oscillations are coherent movements of many electrons that move fast, then more slowly, then faster, and so on. The blue line is the simulation's result, the green line is the analytic GR frame-dragging correction presented in MT92 and MH97, and the red line is the analytic flat spacetime SCLF presented in M74.

Figure 6 shows results of the same run as in Figure 5, but with smaller inclination angle i = 0.1 rad, resulting in acceleration to higher energies. The thick green line is the analytic expression with inclination angle i = 0.1 rad, the thin green line is the analytic expression with inclination angle i = 0.5 rad, and the blue line is the simulation's result.

In Figures 5 and 6, the particles are accelerating "forever," since there is no pair production and no energy loss in the simulation (for now). The first term in ρ_GJ in Equation (6), the cos i term, is dominant close to the star, when frame dragging is important. Here, the acceleration is strong. As the particle's distance from the star increases, acceleration approaches that of flat spacetime, when only the second term in Equation (6), the sin i term, is important (and comparable in curved and flat spacetime); thus, the particle's acceleration becomes more moderate. Pair production will occur as soon as it can occur, depending on the star parameters (e.g., magnetic field strength, field line curvature, star temperature, polar cap radius), which can be changed for each run of the simulation. Results of pair production will be reported in a future paper.

We have programmed a full simulation of particle acceleration above the open field lines of a pulsar's polar cap, using the equation of motion and Maxwell's equations with PIC numerical code.

The first part is a run with flat spacetime SCLF, as presented in M74, FAS77, CR77, and AS79. We got the SCLF condition at the bottom of the tube, i.e., the vanishing of the parallel electric field at the star surface. We also got a cloud of electrons at the bottom of the grid (the star surface), as expected for the SCLF model. This step enables us to calibrate the simulation distances and n_GJ, q, and m to the "real-world" pulsar values, through the tube/polar cap radius R and the plasma frequency ω_p. This step shows the increasing oscillatory pattern as the simulation parameter R ω_p/c increases, eventually resulting in the appearance of instability.

The second step is the run with the GR frame-dragging correction presented in MT92 and MH97. The electrons are accelerated according to the analytic estimates, reaching the Lorentz factor of the order of 10⁶. The particles' Lorentz factor shows an oscillatory pattern that increases with the distance from the star. The oscillations are coherent movements of many electrons that move fast, then more slowly, then faster, and so on.

The simulation shows a "front," which exists only when the particles propagate into the empty tube. The initial vacuum potential drops significantly after the particles enter the tube. Ahead of the "front," where the particles have not yet reached, the potential and fields have the initial values. Behind the "front." where the particles are in the tube, the potential, fields, and the particles' Lorentz factor have the theoretical values expected from the SCLF model (with or without the frame-dragging correction). In the "front," there are strong oscillations of the fields and particles' Lorentz factor.

The third step will be reported in a future paper, and will include pair production. Primary electrons accelerating along the open field lines emit CR photons, or upscatter thermal photons from the star via ICS. The CR or ICS photon can produce a pair via magnetic pair production γB → e⁻e⁺B. If one of the pair members is at a high Landau level, then it radiates its excess energy via SR. The SR photon can also produce a pair. Thus, the cascade of CR/ICS/SR photons and pairs is produced. Most polar cap cascade models assume a steady state—stationary particle outflow. Due to the pair cascade, the magnetosphere is filled with electron–positron plasma, the accelerating electric field is screened, and the magnetosphere is assumed to be a force-free magnetosphere. The current density in those global force-free magnetosphere models is not a free parameter, but rather dictated by the global structure of the magnetosphere. The current density of the force-free magnetosphere models is incongruent with the one-directional Goldreich–Julian current density of the steady-state models. Here, we will study the current adjustment along the tube, using time variability and boundary conditions (one boundary is the star surface and the other boundary is the plasma-free flow). A positron reverse current could be formed, to compensate for the Goldreich–Julian charge density. We will let the simulation show how the current and charge densities in the open magnetic field lines are adjusted to maintain the global structure of the magnetosphere. The multiplicity can be studied as well by keeping the generation number of each formed pair. After completing this stage, we can report the positron return current (thus, we can calculate polar cap heating and X-ray emission), the photons' and electrons' observed spectra (photons and electrons that escape from the magnetosphere after the cascade), the pulsar death line (pulsars without enough pair production), and the PFF height (pair formation front location).

I thank Eduardo Guendelman for useful discussions, valuable suggestions, and encouragement.

I also thank the anonymous referee for the careful review and valuable comments that helped improve this manuscript.