We perform a novel linear simulation of the feedback instability by applying a gyrokinetic model to the magnetosphere, where the geomagnetic field is modeled by a dipole magnetic field. In order to avoid huge numerical calculation costs, the effect of the fluctuating mirror force is neglected. It is found that kinetic effects, such as the finite Larmor radius effect and the electron Landau damping, strongly stabilize the feedback instability and form a peak in the perpendicular wave number spectrum of the growth rate. During the linear growth, it is observed that electron free energy has a much larger value than the electromagnetic field perturbation energy. Since such large electron free energy may lead to electron heating and acceleration, it is important to understand its generation mechanism. Analyses based on the energy conservation reveal that the coupling of the non-uniformity of the equilibrium fields along the geomagnetic field and the fluctuating electron distribution gives rise to thermodynamic power generating the electron free energy. A fine velocity space structure around an average Alfvén velocity of the fluctuating electron distribution mainly contributes to drive the thermodynamic power.

## I. INTRODUCTION

Auroral light emissions are caused by electrons falling into the ionosphere along the geomagnetic field. It is plausible that the strongly bright aurora called the breakup is due to the strong electric current, which is caused by a dynamo mechanism somewhere in the magnetosphere.^{1} On the other hand, various time-varying characteristics of the auroral fine structure suggest that the Alfvén wave traveling along the geomagnetic field is in a strong turbulent state, as indicated by spacecraft observations.^{2} There are various mechanisms for excitation of the Alfvén wave in the magnetosphere: the oscillation of the geomagnetic field triggered by the substorm,^{3} the solar wind pressure disturbance at the magnetopause,^{4} the Kelvin–Helmholtz instability at the magnetospheric boundary,^{5} and the feedback instability that occurs in the magnetosphere–ionosphere coupling system. The feedback instability is categorized into two types: one is due to the standing Alfvén wave called the field line resonance with a low frequency (1–10 mHz) and a long wavelength ( $\u223c105$ km),^{6} and the other is due to the Alfvén wave trapped in the vicinity of the ionosphere called the ionospheric cavity mode (or the ionospheric resonator) with a high frequency (0.1–1 Hz) and a short wavelength ( $\u223c104$ km).^{7}

In the following, we focus on the low-frequency and long-wavelength feedback instability associated with the field line resonance. In order to explain the behavior of the so-called quiet aurora, various studies have been done for the feedback instability.^{8} In the high temperature region of the magnetosphere, kinetic effects, such as the finite Larmor radius (FLR) effect and the electron Landau damping, on the Alfvén wave become important, which is the so-called kinetic Alfvén wave.^{9,10} In addition, the strong non-uniformity along the geomagnetic field could affect the characteristics of the Alfvén wave and the wave–particle interaction with electrons. In the research context considering the non-uniformity, the magnetosphere has often been modeled by the dipole magnetic field. The feedback instability with the Alfvén wave propagating in the dipole magnetic field has been investigated using fluid models.^{11–18} It has been discussed how the wave–particle interaction between the Alfvén wave and electrons occurs in the presence of the non-uniformity of the dipole magnetic field using kinetic models.^{19–21} There has also been discussion about how the non-uniformity of the dipole magnetic field affects the characteristics of the kinetic Alfvén wave, for example, the amplitude of the parallel electric field.^{22–27}

The importance of the electron Landau damping for the Alfvén wave in the magnetosphere has been pointed out.^{10,28,29} In our preceding studies, in order to investigate the kinetic effects on the Alfvén wave, we derive a gyrofluid model of the magnetosphere from the gyrokinetic model, where the geomagnetic field is modeled by a uniform magnetic field.^{30} In the derivation of the gyrofluid model, the dispersion relation of the kinetic Alfvén wave^{10,31} is used to close the electron parallel pressure perturbation. Using the gyrofluid model, it is clarified that the electron Landau damping strongly suppresses the linear growth rate of the feedback instability^{30} and plays an important role in the nonlinear saturation of the feedback instability.^{32}

However, fluid models cannot reveal how the energy is transferred to electrons through the electron Landau damping. In addition, in the presence of the non-uniformity of the geomagnetic field, it is not straightforward to derive the dispersion relation of the kinetic Alfvén wave in a simple algebraic form, which is necessary for constructing the gyrofluid model. These issues motivate the direct application of the gyrokinetic model to the magnetosphere as an initial value problem rather than through the dispersion relation. Then, a simulation method is developed, which solves the gyrokinetic model of the magnetosphere as an initial value problem and properly treats the magnetosphere–ionosphere boundary condition.^{33} For the uniform model of the geomagnetic field, it is confirmed that simulation results given by the above-mentioned method reasonably agree with those predicted by the dispersion relation.^{33} In this study, this methodology is applied to the dipole model of the geomagnetic field.

The remainder of this paper is organized as follows. Models of the magnetosphere and ionosphere are shown in Sec. II. A linear simulation of the feedback instability for the dipole model of the geomagnetic field is performed in Sec. III. It is observed that kinetic effects strongly stabilize the feedback instability and form a peak in the perpendicular wave number spectrum of the growth rate. During the linear growth of the instability, it is observed that the electron free energy is much larger than the energy of the electromagnetic field perturbation. Analyses based on the energy conservation are performed to clarify the generation mechanism of the large electron free energy. A summary is given in Sec. IV.

## II. MODEL EQUATIONS

### A. Nonlinear gyrokinetic model of magnetosphere

We consider a magnetic flux tube along the geomagnetic field. By considering a sufficiently thin flux tube, we can assume that perturbations are periodic in the plane perpendicular to the geomagnetic field. Then, the so-called $\delta f$-gyrokinetic equation is introduced. To focus on the dynamics of the feedback instability in the parallel direction, we neglect perpendicular gradients of the background magnetic field, the equilibrium density, the equilibrium temperature, and the equilibrium pressure for simplicity. Only the perpendicular gradient of the equilibrium electrostatic potential, which is a driving source of the feedback instability, is taken into account.

^{34}

^{35,36}

*e*is the elementary charge, $\u2207\u22a5=\u2207\u2212b0(b0\xb7\u2207)$, $\u2207\u2225=b\xb7\u2207$, $fe=fe0+\delta fe$, $\varphi =\varphi 0+\varphi \u0303$, where $\varphi 0$ is the equilibrium electrostatic potential, $\varphi \u0303$ is the electrostatic potential perturbation, $A\u0303\u2225$ is the perturbation of the parallel component of the vector potential, and

*c*is the velocity of light. Hereafter, the perturbation amplitude in wave number space is represented by the subscript

*k*.

^{33}

^{,}

^{30}For later use, we derive a gyrofluid model: considering the zeroth and first order moments of $v\u2225$ for Eq. (5), the electron continuity equation and the generalized Ohm's law (or the parallel electron momentum equation) are given by

^{30}

*a*and

*b*, $\delta e=c/\omega pe$ is the electron inertial length, and $\omega pe=(4\pi ne0e2/me)1/2$ is the electron plasma frequency. Different from Eq. (29) in Ref. 33, the fluctuating mirror force term characterized by the perpendicular electron density fluctuation is not involved in Eq. (14).

### B. Linearized gyrokinetic model of the magnetosphere

*x*–

*y*plane is the plane perpendicular to the geomagnetic field, and

*z*is the parallel position. In our coordinates, the non-uniformity along the geomagnetic field is considered, but the curvature and perpendicular gradient of the geomagnetic field are omitted. This is consistent with the assumption of $\u2207\u22a5B0=0$, since the curvature is written by $(1/B0)\u2207\u22a5B0$. The ionosphere and magnetic equatorial plane are located at $z=0$ and $z=\u2113$, respectively. However, as mentioned later, the fluid buffer is introduced in the range of $0\u2264z\u2264z0$. Thus, the gyrokinetic model is applied to the range of $z0\u2264z\u2264\u2113$. The conservation of the magnetic flux, i.e., $B0L\u22a52=const.$, is held along the magnetic field, where $L\u22a5(z)$ is the perpendicular scale length depending on the value of

*z*. In our coordinates, the perpendicular wave number is defined by $k\u22a52(z)=kx2(z)+ky2(z)$, where $kx(z)=2\pi p/L\u22a5(z)$ and $ky(z)=2\pi q/L\u22a5(z)$, and ${p,q}$ are the mode numbers in the

*x*and

*y*directions, respectively. Similarly in Ref. 33, but considering the

*z*-dependence of the equilibrium fields, a linearized version of Eq. (5) is given by

Here, we review the dispersion relation of the kinetic Alfvén wave. In the case of the uniform equilibrium along the geomagnetic field, the dispersion relation of the kinetic Alfvén wave [Eq. (46) in Ref. 30] can be derived from Eqs. (11), (12), and (15) by considering $\u2202/\u2202t\u2192\u2212i\omega $ and $\u2202/\u2202z\u2192ik\u2225$, where $\omega $ is the frequency, and $k\u2225$ is the wave number parallel to the geomagnetic field. The derived dispersion relation is basically identical to that in Ref. 10 and involves effects of the ion finite Larmor radius, the electron Landau damping, the electron inertia, and the electron pressure. However, as mentioned in the Introduction, it is not straightforward to derive the dispersion relation of the kinetic Alfvén wave in the presence of the non-uniformity of the equilibrium along the geomagnetic field. Therefore, it is necessary to solve Eqs. (11), (12), and (15) as an initial value problem.

### C. Law of conservation of energy in magnetosphere

^{37,38}where $Eg$ is the energy accompanied by the gyrocenter motion, and $E\varphi $ is the energy associated with the finite Larmor radius effect. In the total energy source $SGK$, $Sg$ is the source of the electron free energy, and $S\varphi $ is the electric power.

### D. Parity separation of fluctuating distribution

### E. Linearized two-fluid model of ionosphere

^{33},

*h*is the ionosphere depth, $\mu P=qiI/miI\nu in$ is the Pedersen mobility, $qiI$ is the ion charge, $miI$ is the ion mass, $\nu in$ is the ion-neutral particle collision frequency, $\mu H=c/B0I$ is the Hall mobility, and $\alpha $ is the recombination coefficient. For the details of the above-mentioned model, see the original literatures, Refs. 39 and 40, and our brief review of the derivation in Ref. 30.

In addition, we mention that recent studies on the feedback instability associated with the ionospheric cavity mode discuss an effect of the vertical shear of the horizontal ion velocity (due to change in the collision frequency of ions and neutral particles along the altitude of the ionosphere).^{41,42} Since the height-integrated ionosphere model is used in this study, such an effect coming from the degree of freedom in the height direction is out of scope.

## III. NUMERICAL SIMULATION

### A. Parameter and equilibrium profile

The geomagnetic field is approximated by the dipole magnetic field. For the dipole magnetic field, the magnitude of the magnetic field is given by $B0=(Bp/2)(RE/R)3(1+3\u2009cos\u20092\theta )1/2$, where *R* is the radial position, $\theta $ is the latitude angle, $Bp=0.6\u2009G$, and $RE$ is the earth radius. By solving a streamline equation for the dipole magnetic field, we consider the coordinate *z* along the magnetic field extending from $R=RE$ and $\theta =\pi /9\u2009rad$ and remind that $z=0$ is the ionosphere and $z=\u2113$ is the magnetic equator, as shown in Fig. 1. For the above-mentioned dipole magnetic field, we obtain $B0|z=0=B0I=0.573\u2009G$, $B0|z=0=B0M=4.82\xd710\u22124\u2009G$, and $\u2113=6.88\xd7104\u2009km$.^{16} Hereafter, the subscript M indicates the value at the magnetic equator.

Almost the same as in Refs. 30 and 33 (except the value of $B0I$), the parameters in the ionosphere are set as follows: $L\u22a5|z=0=LI=100\u2009km$, $\alpha =10\u221213\u2009m3s\u22121$, $h=30\u2009km$, $nI0=5\xd71010\u2009m\u22123$, $\mu H=1.75\xd7104\u2009m2V\u22121s\u22121$, $\mu P=0.25\mu H$, $EI0=\u2212\u2207\u22a5\varphi I0=EI0y\u0302$, and $EI0=20 mV m\u22121$. Similarly, using the same density and temperature as in Refs. 30 and 33, the parameters in the magnetosphere are set as follows: $vAM=1.05\xd7103 km s\u22121$, $TeM=TiM=600\u2009eV$, $vteM=1.03\xd7104 km s\u22121$, $n0M=1\u2009cm\u22123$, $L\u22a5|z=\u2113=LM=LIB0I/B0M$, $\u2207\u22a5\varphi 0=\u2207\u22a5\varphi I0$, $\rho iM=51.8\u2009km$, and $\delta eM=5.31\u2009km$. For convenience, we define $\tau A=\u2113/vAM=65.5\u2009s$.

Following Ref. 21, we use the observations in Ref. 43 for the density profile along the geomagnetic field in the magnetosphere. According to the observations in Ref. 43, the density is approximately proportional to the magnitude of the geomagnetic field. There is no reliable observation for the temperature, but it is roughly several hundred to a thousand eV at the magnetic equator and a few eV near the Earth region, which is roughly inversely proportional to the magnitude of the geomagnetic field. Therefore, we model $n0\u221dB0$ and $Te=Ti\u221dB0\u22121$, and these are consistent with the assumption of Eq. (2). Taking into account the relation of $B0L\u22a52=B0MLM2$, the *z* dependence of the main parameters in the magnetosphere is given as follows: $vA(z)/vAM=(B0(z)/B0M)1/2$, $vte(z)/vteM=(B0M/B0(z))1/2$, $\rho i(z)/L\u22a5(z)=(\rho iM/LM)(B0M/B0(z))$, and $\delta e(z)/L\u22a5(z)=\delta eM/LM$. Figure 2 shows the *z* profiles of the Alfvén velocity and the electron thermal velocity.

### B. Linear simulation

Numerical methods for the buffer layer near the ionosphere and boundary conditions are described in the Appendix. In the following, linear simulations are performed using Eqs. (38) and (39) for the magnetosphere model, Eqs. (43) and (44) for the ionosphere model, and Eqs. (16), (17), and (A1) for the buffer layer model. Note that, to facilitate numerical calculation, the equations in the rest frame of the $E\xd7B$ drift are used. As in Ref. 33, in order to calculate $\u2202A\u0303\u2225,k/\u2202t$ on the right-hand side of Eq. (39), Eq. (17) is used. The time integral is calculated by the fourth-order Runge–Kutta method. The time step is $\Delta t/\tau A=10\u22124$. The *z* derivative is calculated by the finite difference method. The integral in the $v\u2225$ space is calculated by the trapezoidal rule. A step size in the *z* direction is $\Delta z/\u2113=1/256$. Considering $vteM/vAM=10.3$, a numerical calculation range in the $v\u2225$ direction is chosen to be $\u221250\u2264v\u2225/vAM\u226450$. Since the required resolution in the $v\u2225$ direction changes along *z*, applying the smallest resolution to the entire magnetosphere increases the numerical cost. In order to avoid this issue, unequal interval step sizes in the $v\u2225$ direction $\Delta v\u2225i=(1+a)i\Delta v\u22250\u2009(i=0,1,2,\u2026,N\u22121)$ are used, where $a=0.012\u2009077\u20098$, $N=400$, and $\Delta v\u22250=5\xd710\u22123vAM$.

Figure 3 shows the time evolution of real and imaginary parts of $n\u0303I,k$ for $(kx,ky)=(0,40)2\pi /L\u22a5$. It is observed that the initial perturbation grows exponentially with the frequency, i.e., there exists the phase difference between the real and imaginary parts.

Figure 4 shows the time evolution of $Eg$, $E\varphi $, and $EA\u2225$ for $(kx,ky)=(0,40)2\pi /L\u22a5$. Similar to Fig. 3, each energy also grows exponentially. It is found that the electron free energy $Eg$ is much larger than the electromagnetic field perturbation energy $E\varphi +EA\u2225$. This is a prominent feature, since $Eg$ and $E\varphi +EA\u2225$ are comparable in the case of the uniform magnetic field.^{33}

Figure 5 shows the time evolution of $\u222bSgdt$ and $\u222bS\varphi dt$ for $(kx,ky)=(0,40)2\pi /L\u22a5$. We confirmed that $\u222bSgdt+\u222bS\varphi dt$ agrees well with $Eg+E\varphi +EA\u2225$, which indicates that the energy conservation law Eq. (19) is satisfied. In Fig. 5, it is observed that $\u222bSgdt$ is much larger than $\u222bS\varphi dt$. Therefore, it is found that the generation of the large $Eg$ in Fig. 4 is mainly caused by $Sg$.

Figure 6 shows the $ky$ dependence of the growth rate and frequency for $kx=0$. In Fig. 6, linear simulation results using the gyrokinetic model are represented by dots, and those using the MHD model [Eqs. (16), (17), and (A1) in the limit of $k\u22a5\rho i\u21920$] are represented by a solid line. Similarly in the case with the uniform magnetic field equilibrium,^{30} it is observed that the growth rate has a peak if kinetic effects are taken into account. However, a difference from the uniform equilibrium case is that the peak locates at the higher $ky$ region. This is probably due to the increase in the average Alfvén velocity and the shorter length of the magnetic field line, which reduces the propagation time of the Alfvén wave along the geomagnetic field. In other words, for the same equilibrium electric field, only high wavenumber perturbations can resonate with the Alfvén waves with the short timescales. On the other hand, it is observed that the frequency does not change significantly even in the presence of the kinetic effects.

Figure 7 shows the contours of real and imaginary parts of $\delta ge,k$ in the *z*- $v\u2225$ space for $(kx,ky)=(0,40)(2\pi /L\u22a5)$ at $t/\tau A=10$. Similarly, in the case with the uniform magnetic field equilibrium,^{30} a fine structure (the ballistic mode) in the $v\u2225$ direction is observed, which is a characteristic of the phase mixing causing the Landau damping. The *z* dependence differs largely depending on the value of $v\u2225$. On the other hand, although the thermal velocity greatly changes along the geomagnetic field, the region where $\delta ge,k$ exists in the $v\u2225$ space hardly changes along the geomagnetic field. In particular, there exists the strong peak in the region of $|v\u2225|/vAM\u22723$. These imply that the velocity space structure is mainly dominated by the phase mixing, rather than by the thermal velocity. In addition, an average value of the Alfvén velocity is given by $v\xafA=(1/\u2113)\u222b0\u2113vAdz=1.96vAM$.

Figure 8 shows the *z* profiles of real and imaginary parts of $n\u0303e,k$, $j\u0303\u2225k$, and $p\u0303\u2225e,k$ for $(kx,ky)=(0,40)(2\pi /L\u22a5)$ at $t/\tau A=10$. The region of $z/\u2113\u22640.2$ is the FLR-MHD buffer layer, and that of $z/\u2113\u22650.2$ is the gyrokinetic model layer. It is observed that the profiles are smoothly connected at $z/\u2113=0.2$. It is found that the *z* dependence of $n\u0303e,k$, $j\u0303\u2225k$, and $p\u0303\u2225e,k$ are fairly different, even taking into account that each perturbation is normalized. This is consistent with the characteristic that the *z* dependence differs largely depending on the value of $v\u2225$ in Fig. 7.

*n*can be easily obtained by direct calculations using Mathematica, for example. The probabilist's Hermite polynomials satisfy the orthogonality of $\u222b\u2212\u221e\u221eHenHem\u2009exp(\u2212v2/2)dv=(2\pi )1/2n!\delta m,n$, where $\delta m,n$ is the Kronecker delta. Using the orthogonality, $an$ in Eq. (45) can be solved as follows:

*n*increases, $an$ includes higher-order fluid perturbations. Finally, by substituting Eq. (45) into Eq. (22) and considering the orthogonality, $Eg$ is expanded as follows:

Figure 9 shows the Hermite spectrum of $Eg,n$ for $(kx,ky)=(0,40)2\pi /L\u22a5$ at $t/\tau A=1$, $t/\tau A=5$, and $t/\tau A=10$. In Fig. 9, at $t/\tau A=10$, where the linear growth rate is converged, the lower-order $n=0,1,2$ perturbations, i.e., the electron density perturbation, the electron parallel velocity perturbation, and the electron parallel pressure perturbation, are found to have almost the same energy as the electromagnetic field perturbation in Fig. 4. It is observed that the Hermite spectrum has a long tail in the higher-order $n\u22653$ region, although each energy is not large. Therefore, it is clarified that the large electron free energy corresponds to the state where the energy is distributed to the higher-order fluid perturbations.

*z*profile of the thermodynamic power density $Dg$ for $(kx,ky)=(0,40)2\pi /L\u22a5$ at $t/\tau A=8$, $t/\tau A=9$, and $t/\tau A=10$. It is observed that $Dg$ is positive in most

*z*regions, although $Dg$ is negative in the small

*z*region. In the integrand of Eq. (42), a combination of equilibrium quantities $(Tev\u2225/ge0)L\u22a52B02$ is always a positive real value. Therefore, whether the integrand contributes positive or negative depends on the sign of the following quantity:

Figure 11 shows the contour of $\sigma g$ in *z*- $v\u2225$ space for $(kx,ky)=(0,40)2\pi /L\u22a5$ at $t/\tau A=10$. Comparing Figs. 10 and 11, it is found that the sign relationship between $Dg$ and $\sigma g$ is almost the same. In Fig. 11, it is observed that there exists the prominent peak in the range of $v\u2225/vAM\u22723$. Therefore, it is clarified that the fine velocity space structure around an average Alfvén velocity ( $\u223c2vAM$) of the fluctuating electron distribution mainly obtain the electron free energy from the parallel gradients of the equilibrium fields. The generation of the thermodynamic power enhances the localized velocity space structure as shown in Fig. 7. Such localized velocity space structure is observed as the generation of the higher-order fluid perturbation energies, as shown in Fig. 9, which is similar to how the high-wave number Fourier components appear around the spatially localized structure.

Finally, Fig. 12 schematizes the physical mechanism of the large electron free energy generation. In the feedback instability, the kinetic Alfvén wave is excited by the electric power to the magnetosphere through the magnetosphere–ionosphere coupling. In the development of the kinetic Alfvén wave, the fine velocity space structure appears due to the phase mixing as shown in Fig. 7. Even in the uniform equilibrium limit, such a fine structure exists, and the thermodynamic power is driven at the magnetosphere–ionosphere boundary.^{33} Contrarily, when the equilibrium is non-uniform, the coupling of the fine structure and the equilibrium non-uniformity drives the positive thermodynamic power in most regions of the magnetosphere, as shown in Fig. 10. As a result, the electron free energy has a large value compared with the electromagnetic perturbation energy as shown in Fig. 4. Since the electron free energy is generated for electrons having velocities characterized by the wave–particle interaction, rather than thermal velocity as shown in Fig. 11, the free energy is stored as those of higher-order fluid perturbations as shown in Fig. 9.

## IV. SUMMARY

The method applying the gyrokinetic model to the magnetosphere as an initial value problem is revisited. The magnetosphere–ionosphere coupling method developed for the uniform magnetic field model of the geomagnetic field in our previous study is extended to the dipole magnetic field model of the geomagnetic field.

It is numerically difficult to completely resolve the velocity space for the entire magnetosphere, because the density and the temperature largely change along the geomagnetic field. For this reason, we introduce a gyrokinetic model pre-integrated for the magnetic moment space. Furthermore, the fine structure in the parallel velocity space due to the phase mixing causing the Landau damping is numerically incompatible with the parallel velocity differentiation involved in the fluctuating mirror force term. In order to resolve this issue, we introduce an approximate form of the fluctuating distribution that does not explicitly give rise to the fluctuating mirror force term.

Using the above-mentioned magnetosphere model, we consider the magnetosphere–ionosphere coupling and simulate the linear evolution of the feedback instability. Similar to the previous studies applying a uniform magnetic field model to the geomagnetic field, it is confirmed that the linear growth rate is strongly suppressed by the kinetic effects, which give rise to the peak of the growth rate in the perpendicular wave number spectrum. On the other hand, a major difference from the case using the uniform magnetic field is that the electron free energy is much larger than the electromagnetic field perturbation energy in the linear growth of the feedback instability. As a result of the analysis using the Hermitian expansion of the fluctuating distribution, it is found that the electron free energy is distributed not only to the lower-order fluid perturbations such as the density, parallel velocity, and parallel pressure, but also to higher-order fluid perturbations. It is newly revealed that the thermodynamic power generating the electron free energy is driven by combining the equilibrium non-uniformity along the geomagnetic field with the fine velocity space structure of the fluctuating distribution due to the phase mixing. It is also found that the distribution function perturbation has a localized structure in velocity space, which is due to the intensive free energy transfer to electrons with velocities about the average Alfvén velocity. Such a localized structure in velocity space is observed as the appearance of the high-order fluid perturbations. The above results suggest that the non-uniformity in the magnetosphere plays an important role in generating the electron free energy. The electron free energy could be converted into the electron kinetic energy, i.e., acceleration or heating of electrons, by some mechanism, such as energy dissipation mechanisms or the mirror force.

In this study, the fluctuating mirror force is omitted for the numerical calculation reasons. It is necessary to consider how to deal with the magnetic moment space, which requires extreme numerical costs. In addition, the buffer region is introduced near the ionosphere to mimic the collisional effect. The validity of this approximation needs to be checked. Furthermore, by extending this study, it is possible to perform nonlinear simulations to investigate how the non-uniformity in the magnetosphere affects the nonlinear evolution of the feedback instability, that is, the spontaneous structure formation of the aurora. In order to more accurately predict the frequency of the feedback instability, it is possible to introduce a more realistic stretched geomagnetic field model taking into account the effects of the solar wind.^{45,46} By considering the non-uniformity perpendicular to the geomagnetic field, i.e., perpendicular gradients of equilibrium fields, it is possible to examine the effects of magnetic and diamagnetic drifts on the feedback instability, although those effects are estimated to be minor. These are the issues to be addressed in future work.

## ACKNOWLEDGMENTS

The computation in this work has been done using the facilities of the Center for Cooperative Work on Computational Science, University of Hyogo. S.N. would like to acknowledge the Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu University. This work was supported by JSPS KAKENHI Grant No. 21K03502 [Grant-in-Aid for Scientific Research (C)].

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Seiya Nishimura:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (lead). **Ryusuke Numata:** Methodology (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: BUFFER LAYER AND BOUNDARY CONDITION

^{47}extends the entire magnetosphere, then, the $v\u2225$ profile of $\delta ge,k$ greatly deviates from that of $ge0$ in the vicinity of the ionosphere. Then, near the ionosphere, the thermodynamic power density $Dg$ takes a very large value, and the higher-order fluid perturbations become dominant. We consider that such phenomenon is unphysical. This is because, in a more realistic situation, the collision between particles becomes effective near the ionosphere, and the generation of the fine structure in the $v\u2225$ space due to the ballistic mode would be greatly suppressed. In fact, the ionosphere is well described by the two-fluid model without the higher-order fluid perturbations. As an example for evaluating the effect of the collision, for our parameters and equilibrium profiles, the magnetic Reynolds number $S=4\pi L\u22a52/c2\tau A\eta \u2225$ ( $\eta \u2225$ is the parallel resistivity) is $S\u223c3.2\xd71012$ at the magnetic equator ( $z=\u2113$), while, near the ionosphere, $S\u223c4.2\xd7108\u2009(z=0.2\u2113)$, $S\u223c1.7\xd7107\u2009(z=0.1\u2113)$, and $S\u223c6.5\xd7104\u2009(z=0)$. Therefore, the vicinity of the magnetic equator is in the collisionless state, while the region of $0\u2264z\u22640.2\u2113$ is considered to be in the weakly collisional state, where the collision operator in the gyrokinetic model works effectively. In this study, we introduce a buffer region in the range of $0\u2264z\u2264z0=0.2\u2113$ to take into account the effect of the collisional dissipation on the fine structure in the $v\u2225$ space, rather than introducing a collision operator in the gyrokinetic model. In the buffer region, we apply the gyrofluid model, where $p\u0303\u2225e,k$ is given by the following relation: