May 29, 2026
We calculate the entropy of a Dirac quantum field near a static spherically symmetric black hole in \(f(Q)\) gravity by combining two ingredients. The first ingredient is the residue-based Robson–Villari–Biancalana method, in which the Hawking temperature is written as a surface-gravity term plus a residue-induced correction. The second ingredient is the thin-film or brick-wall state-counting method for a Dirac field in a near-horizon background. Starting from the \(f(Q)\)-deformed metric \[\mathrm{d}s^{2}=-g(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{g(r)}+r^{2}\mathrm{d}\Omega^{2},\] we derive the Hamilton–Jacobi equation for the Dirac field, obtain the radial momentum, count the fermionic modes, and compute the free energy and entropy at the residue-corrected temperature. The final result shows that the Dirac-field entropy is still proportional to the horizon area after a proper cutoff is introduced, while the proportionality factor is corrected by the RVB residue through a cubic temperature factor. For the quadratic model \(f(Q)=Q+\alpha Q^{2}\), an explicit closed expression is obtained.
Keywords: \(f(Q)\) gravity, Dirac field, black hole entropy, RVB method, residue theorem, brick-wall model, thin-film model, nonmetricity.
Black hole entropy is one of the most direct bridges among gravity, quantum theory, and thermodynamics. The standard area law was proposed by Bekenstein and Hawking and was later related to Euclidean methods, Noether charge, and horizon thermodynamics [1]–[6]. Quantum matter fields near the horizon provide another statistical route to black hole entropy through the brick-wall or thin-film method [7]–[10]. In this picture, the entropy comes from the large density of quantum states close to the horizon, and a near-horizon cutoff is introduced to regulate the ultraviolet divergence.
The Dirac-field version of this calculation is especially useful because the leading WKB equation of a spinor field reduces to a Hamilton–Jacobi equation. Therefore, at leading order, the radial momentum of the Dirac field is controlled by the same near-horizon pole structure as scalar fields, while the fermionic statistics and spin degeneracy change the numerical coefficient [11]–[14]. This is the reason why the thin-film calculation of a Dirac field still leads to an entropy proportional to the horizon area after a suitable cutoff is chosen.
In the modified-gravity sector, \(f(Q)\) gravity is a nonmetricity-based extension of symmetric teleparallel gravity [15]–[18]. Its gravitational dynamics are governed by a function of the nonmetricity scalar \(Q\). Black hole thermodynamics in such nonmetricity-based theories is subtle, because entropy may be interpreted either through an effective first-law reconstruction or through a Noether-charge-like analysis [19]–[21].
The RVB method proposed by Robson, Villari, and Biancalana gives a complex-analytic way of extracting the Hawking temperature from the residue of a horizon pole [22]. In recent applications to \(f(Q)\) black holes, the RVB temperature can be written as the usual surface-gravity temperature plus a residue-induced shift [23], [24]. The purpose of this paper is to combine that RVB-corrected temperature with the Dirac thin-film entropy calculation. The result should be interpreted as a test-field, near-horizon, thermodynamic entropy branch rather than a universal Noether-charge theorem.
We use natural units \[G=\hbar=c=k_{\rm B}=1.\] The gravitational action of \(f(Q)\) gravity coupled to a Dirac field is written as \[I=\frac{1}{16\pi}\int \mathrm{d}^{4}x\sqrt{-g}\,f(Q)+I_{\rm D},\] where \(Q\) is the nonmetricity scalar and \(I_{\rm D}\) is the Dirac-field action. In the static spherically symmetric sector, we take \[\mathrm{d}s^{2}=-g(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{g(r)}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2}\right).\] The event horizon radius \(r_{+}\) is determined by \[g(r_{+})=0.\]
Following the residue-corrected \(f(Q)\) setup, we parameterize the metric function as \[g(r)=1-\frac{2M}{r}+\psi_{Q}(r),\] where \(\psi_{Q}(r)\) encodes the \(f(Q)\)-induced deformation. The horizon equation gives \[1-\frac{2M}{r_{+}}+\psi_{Q}(r_{+})=0.\] Solving this equation for \(M\), we obtain \[M(r_{+})=\frac{r_{+}}{2}\left[1+\psi_{Q}(r_{+})\right].\] Differentiating \(M(r_{+})\) with respect to \(r_{+}\), we obtain \[\frac{\mathrm{d}M}{\mathrm{d}r_{+}} = \frac{1}{2}\left[1+\psi_{Q}(r_{+})+r_{+}\psi_{Q}'(r_{+})\right].\] For later convenience, define \[\Xi(r_{+}) = 1+\psi_{Q}(r_{+})+r_{+}\psi_{Q}'(r_{+}).\] Then \[\frac{\mathrm{d}M}{\mathrm{d}r_{+}} = \frac{1}{2}\Xi(r_{+}).\] The radial derivative of the metric function at the horizon is \[g'(r_{+}) = \frac{1+\psi_{Q}(r_{+})}{r_{+}}+\psi_{Q}'(r_{+}).\] Using the definition of \(\Xi(r_{+})\), this can be rewritten as \[g'(r_{+}) = \frac{\Xi(r_{+})}{r_{+}}.\] Therefore the ordinary surface-gravity temperature is \[T_{0}(r_{+}) = \frac{g'(r_{+})}{4\pi} = \frac{\Xi(r_{+})}{4\pi r_{+}}.\]
In the RVB method, a simple horizon zero of \(g(r)\) becomes a simple pole of \(1/g(z)\) after complexification \(r\rightarrow z\). Hence \[\operatorname*{Res}_{z=r_{+}}\frac{1}{g(z)} = \frac{1}{g'(r_{+})}.\] The corresponding inverse Hawking temperature can be written as \[\beta_{0} = 4\pi \operatorname*{Res}_{z=r_{+}}\frac{1}{g(z)} = \frac{4\pi}{g'(r_{+})}.\] Thus \[T_{0} = \beta_{0}^{-1} = \frac{g'(r_{+})}{4\pi}.\]
For a more general residue correction, define a complex function \(F(z)\) associated with the metric function, the nonmetricity scalar, or another horizon-data function. The winding number is \[N_{\Gamma} = \frac{1}{2\pi \mathrm{i}}\oint_{\Gamma}\frac{F'(z)}{F(z)}\,\mathrm{d}z = \sum_{z_{k}\in \Gamma} \operatorname*{Res}_{z=z_{k}}\frac{F'(z)}{F(z)}.\] The residue-induced temperature shift is written as \[C_{\rm res} = \lambda_{\rm R} N_{\Gamma},\] where \(\lambda_{\rm R}\) is a normalization constant carrying the dimension of temperature. The RVB-corrected temperature is then \[T_{\rm RVB}(r_{+}) = \frac{g'(r_{+})}{4\pi}+C_{\rm res}.\] Using \(g'(r_{+})=\Xi(r_{+})/r_{+}\), this becomes \[T_{\rm RVB}(r_{+}) = \frac{\Xi(r_{+})}{4\pi r_{+}}+C_{\rm res}.\]
Before computing the Dirac-field entropy, it is useful to recall the geometric entropy branch generated by the same RVB temperature. The first law is assumed to be \[\mathrm{d}M=T_{\rm RVB}\,\mathrm{d}S_{\rm geo}.\] Therefore \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = \frac{1}{T_{\rm RVB}(r_{+})} \frac{\mathrm{d}M}{\mathrm{d}r_{+}}.\] Substituting \[\frac{\mathrm{d}M}{\mathrm{d}r_{+}} = \frac{1}{2}\Xi(r_{+})\] and \[T_{\rm RVB}(r_{+}) = \frac{\Xi(r_{+})}{4\pi r_{+}}+C_{\rm res},\] we obtain \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = \frac{\frac{1}{2}\Xi(r_{+})}{\frac{\Xi(r_{+})}{4\pi r_{+}}+C_{\rm res}}.\] Multiplying numerator and denominator by \(4\pi r_{+}\), one obtains \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = \frac{2\pi r_{+}\,\Xi(r_{+})}{\Xi(r_{+})+4\pi C_{\rm res}r_{+}}.\] Hence the residue-corrected geometric entropy is \[S_{\rm geo}(r_{+}) = \int^{r_{+}} \frac{2\pi u\,\Xi(u)}{\Xi(u)+4\pi C_{\rm res}u}\,\mathrm{d}u +S_{0}.\] When \(C_{\rm res}=0\), the integrand reduces to \(2\pi u\), so \[S_{\rm geo}(r_{+}) = \pi r_{+}^{2}+S_{0}.\] Choosing \(S_{0}=0\), we recover the Bekenstein–Hawking area law \[S_{\rm BH} = \frac{A_{+}}{4}, \qquad A_{+}=4\pi r_{+}^{2}.\]
For a small residue shift satisfying \[\left|4\pi C_{\rm res}r_{+}\right|\ll \left|\Xi(r_{+})\right|,\] we expand the denominator as \[\frac{1}{\Xi(r_{+})+4\pi C_{\rm res}r_{+}} = \frac{1}{\Xi(r_{+})} \left[ 1-\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right] +O(C_{\rm res}^{2}).\] Therefore \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = 2\pi r_{+} \left[ 1-\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right] +O(C_{\rm res}^{2}).\] After integration, the first-order entropy is \[S_{\rm geo}(r_{+}) = \frac{A_{+}}{4} - 8\pi^{2}C_{\rm res} \int^{r_{+}} \frac{u^{2}}{\Xi(u)}\,\mathrm{d}u +O(C_{\rm res}^{2}).\]
The Dirac action in curved spacetime is \[I_{\rm D} = \int \mathrm{d}^{4}x\sqrt{-g}\, \bar{\Psi} \left[ \mathrm{i}\gamma^{a}e_{a}^{\;\mu}D_{\mu}-\mu \right]\Psi,\] where \(\mu\) is the Dirac mass, \(e_{a}^{\;\mu}\) is the tetrad, and \(D_{\mu}\) is the spinor covariant derivative. A convenient tetrad for the metric is \[e_{a}^{\;\mu} = \operatorname{diag} \left( \frac{1}{\sqrt{g(r)}}, \sqrt{g(r)}, \frac{1}{r}, \frac{1}{r\sin\theta} \right).\] The Dirac equation is \[\left[ \mathrm{i}\gamma^{a}e_{a}^{\;\mu}D_{\mu}-\mu \right]\Psi=0.\]
To obtain the leading WKB equation, take the ansatz \[\Psi = a(x)\exp\left(\frac{\mathrm{i}}{\hbar}I(x)\right),\] where \(a(x)\) is a slowly varying spinor amplitude. Keeping only the leading order in \(\hbar\), the spin connection does not contribute to the principal Hamilton–Jacobi equation. Thus \[\left[ \gamma^{a}e_{a}^{\;\mu}\partial_{\mu}I-\mu \right]a(x)=0.\] Multiplying by \[\left[ \gamma^{b}e_{b}^{\;\nu}\partial_{\nu}I+\mu \right],\] and using the Clifford algebra, we obtain \[\left[ g^{\mu\nu}\partial_{\mu}I\partial_{\nu}I+\mu^{2} \right]a(x)=0.\] For a nontrivial spinor amplitude, the Hamilton–Jacobi equation is therefore \[g^{\mu\nu}\partial_{\mu}I\partial_{\nu}I+\mu^{2}=0.\]
Separate the action as \[I = -Et+W(r)+\Theta(\theta,\phi).\] The angular part is represented by the separation constant \(L^{2}\), namely \[\left(\partial_{\theta}\Theta\right)^{2} + \frac{1}{\sin^{2}\theta} \left(\partial_{\phi}\Theta\right)^{2} = L^{2}.\] Using the inverse metric components, the Hamilton–Jacobi equation becomes \[-\frac{E^{2}}{g(r)} + g(r)\left(\frac{\mathrm{d}W}{\mathrm{d}r}\right)^{2} + \frac{L^{2}}{r^{2}} + \mu^{2} = 0.\] Solving for the radial momentum gives \[k_{r}^{2}(r,E,L) = \left(\frac{\mathrm{d}W}{\mathrm{d}r}\right)^{2} = \frac{E^{2}}{g^{2}(r)} - \frac{1}{g(r)} \left( \mu^{2}+\frac{L^{2}}{r^{2}} \right).\] In the semiclassical angular-momentum approximation, \[L^{2}\simeq l(l+1).\] Thus \[k_{r}^{2}(r,E,l) = \frac{E^{2}}{g^{2}(r)} - \frac{1}{g(r)} \left[ \mu^{2}+\frac{l(l+1)}{r^{2}} \right].\]
The radial WKB quantization condition is \[n_{r}\pi = \int k_{r}(r,E,l)\,\mathrm{d}r.\] The total number of modes below energy \(E\) is therefore \[n(E) = \frac{g_{\rm D}}{\pi} \int_{r_{+}+\epsilon}^{r_{+}+\epsilon+\delta}\mathrm{d}r \int_{0}^{l_{\max}} (2l+1) k_{r}(r,E,l)\,\mathrm{d}l,\] where \(g_{\rm D}\) is the Dirac degeneracy. For a four-component Dirac spinor one may take \[g_{\rm D}=4.\] For a single two-spin-state fermion sector one may instead set \[g_{\rm D}=2.\]
The condition \(k_{r}^{2}\geq 0\) gives \[\frac{E^{2}}{g^{2}(r)} - \frac{1}{g(r)} \left[ \mu^{2}+\frac{l(l+1)}{r^{2}} \right] \geq 0.\] Multiplying by \(g^{2}(r)\), we get \[E^{2} - g(r)\mu^{2} - g(r)\frac{l(l+1)}{r^{2}} \geq 0.\] Therefore \[l_{\max}(l_{\max}+1) = r^{2} \left[ \frac{E^{2}}{g(r)}-\mu^{2} \right].\]
Now define \[y=l(l+1).\] Then \[\mathrm{d}y=(2l+1)\mathrm{d}l.\] The angular integral becomes \[\int_{0}^{l_{\max}} (2l+1) k_{r}(r,E,l)\,\mathrm{d}l = \int_{0}^{y_{\max}} \left[ \frac{E^{2}}{g^{2}(r)} - \frac{\mu^{2}}{g(r)} - \frac{y}{g(r)r^{2}} \right]^{1/2} \mathrm{d}y.\] Let \[A(r,E) = \frac{E^{2}}{g^{2}(r)} - \frac{\mu^{2}}{g(r)}.\] Let \[B(r) = \frac{1}{g(r)r^{2}}.\] Then \[y_{\max} = \frac{A(r,E)}{B(r)}.\] The integral is \[\int_{0}^{A/B} \left[A-By\right]^{1/2}\mathrm{d}y = \frac{2}{3B}A^{3/2}.\] Substituting \(A\) and \(B\), we find \[\int_{0}^{l_{\max}} (2l+1) k_{r}(r,E,l)\,\mathrm{d}l = \frac{2}{3} \frac{r^{2}}{g^{2}(r)} \left[ E^{2}-\mu^{2}g(r) \right]^{3/2}.\] Therefore the number of Dirac modes is \[n(E) = \frac{2g_{\rm D}}{3\pi} \int_{r_{+}+\epsilon}^{r_{+}+\epsilon+\delta} \frac{r^{2}}{g^{2}(r)} \left[ E^{2}-\mu^{2}g(r) \right]^{3/2} \mathrm{d}r.\]
Near the horizon, \[g(r)\rightarrow 0.\] Hence the mass term is subleading in the leading ultraviolet entropy. Thus \[\left[ E^{2}-\mu^{2}g(r) \right]^{3/2} = E^{3}+O(g).\] The leading mode number is \[n(E) = \frac{2g_{\rm D}E^{3}}{3\pi} I_{Q} +O(g),\] where \[I_{Q} = \int_{r_{+}+\epsilon}^{r_{+}+\epsilon+\delta} \frac{r^{2}}{g^{2}(r)}\,\mathrm{d}r.\]
For fermions, the free energy is \[F_{\rm D} = -\frac{1}{\beta} \sum_{s}\ln\left(1+e^{-\beta E_{s}}\right).\] In the continuum approximation, integration by parts gives \[F_{\rm D} = -\int_{0}^{\infty} \frac{n(E)}{e^{\beta E}+1}\,\mathrm{d}E.\] For the RVB-corrected thermal atmosphere, we take \[\beta=\beta_{\rm RVB}=\frac{1}{T_{\rm RVB}}.\] Substituting the leading mode number, \[F_{\rm D}^{\rm RVB} = - \frac{2g_{\rm D}I_{Q}}{3\pi} \int_{0}^{\infty} \frac{E^{3}}{e^{\beta_{\rm RVB}E}+1}\,\mathrm{d}E.\] The standard fermionic integral is \[\int_{0}^{\infty} \frac{E^{3}}{e^{\beta E}+1}\,\mathrm{d}E = \frac{7\pi^{4}}{120\beta^{4}}.\] Therefore \[F_{\rm D}^{\rm RVB} = - \frac{7g_{\rm D}\pi^{3}}{180} \frac{I_{Q}}{\beta_{\rm RVB}^{4}}.\] The entropy is \[S_{\rm D}^{\rm RVB} = \beta_{\rm RVB}^{2} \frac{\partial F_{\rm D}^{\rm RVB}}{\partial \beta_{\rm RVB}}.\] Using \[\frac{\partial}{\partial\beta_{\rm RVB}} \left( -\frac{1}{\beta_{\rm RVB}^{4}} \right) = \frac{4}{\beta_{\rm RVB}^{5}},\] we get \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}\pi^{3}}{45} \frac{I_{Q}}{\beta_{\rm RVB}^{3}}.\] Since \(\beta_{\rm RVB}^{-1}=T_{\rm RVB}\), the Dirac-field entropy becomes \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}\pi^{3}}{45} I_{Q} \left[ \frac{g'(r_{+})}{4\pi}+C_{\rm res} \right]^{3}.\]
This is the central Dirac-field result before rewriting the cutoff in proper-distance form. The residue correction enters through the cubic factor of the temperature because the free energy of a massless fermionic atmosphere scales as \(T^{4}\), while entropy scales as \(T^{3}\).
Assume that the horizon is nonextremal, so that \[g'(r_{+})\neq 0.\] Near the horizon, \[g(r) = g'(r_{+})(r-r_{+}) +O\left((r-r_{+})^{2}\right).\] Therefore \[g^{2}(r) = \left[g'(r_{+})\right]^{2}(r-r_{+})^{2} +O\left((r-r_{+})^{3}\right).\] In the thin-film region, \[r^{2} = r_{+}^{2}+O(r-r_{+}).\] Thus \[I_{Q} = \int_{r_{+}+\epsilon}^{r_{+}+\epsilon+\delta} \frac{r^{2}}{g^{2}(r)}\,\mathrm{d}r = \frac{r_{+}^{2}}{\left[g'(r_{+})\right]^{2}} \int_{\epsilon}^{\epsilon+\delta} \frac{\mathrm{d}x}{x^{2}} +\cdots,\] where \[x=r-r_{+}.\] The elementary integral is \[\int_{\epsilon}^{\epsilon+\delta} \frac{\mathrm{d}x}{x^{2}} = \frac{1}{\epsilon} - \frac{1}{\epsilon+\delta}.\] Therefore \[I_{Q} = \frac{r_{+}^{2}}{\left[g'(r_{+})\right]^{2}} \left[ \frac{1}{\epsilon} - \frac{1}{\epsilon+\delta} \right] +\cdots.\] Substituting this into the entropy formula gives the coordinate-cutoff expression \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}\pi^{3}}{45} \frac{r_{+}^{2}}{\left[g'(r_{+})\right]^{2}} \left[ \frac{1}{\epsilon} - \frac{1}{\epsilon+\delta} \right] \left[ \frac{g'(r_{+})}{4\pi}+C_{\rm res} \right]^{3} +\cdots.\]
Now introduce the proper near-horizon cutoffs \[h = \int_{r_{+}}^{r_{+}+\epsilon} \frac{\mathrm{d}r}{\sqrt{g(r)}} = 2\sqrt{\frac{\epsilon}{g'(r_{+})}}+\cdots,\] and \[H = \int_{r_{+}}^{r_{+}+\epsilon+\delta} \frac{\mathrm{d}r}{\sqrt{g(r)}} = 2\sqrt{\frac{\epsilon+\delta}{g'(r_{+})}}+\cdots.\] These relations imply \[\frac{1}{\epsilon} = \frac{4}{g'(r_{+})h^{2}},\] and \[\frac{1}{\epsilon+\delta} = \frac{4}{g'(r_{+})H^{2}}.\] Hence \[\frac{1}{\epsilon} - \frac{1}{\epsilon+\delta} = \frac{4}{g'(r_{+})} \left( \frac{1}{h^{2}}-\frac{1}{H^{2}} \right).\] Define \[\Delta_{h} = \frac{1}{h^{2}}-\frac{1}{H^{2}}.\] Then \[I_{Q} = \frac{4r_{+}^{2}}{\left[g'(r_{+})\right]^{3}} \Delta_{h} +\cdots.\] The proper-cutoff entropy becomes \[S_{\rm D}^{\rm RVB} = \frac{28g_{\rm D}\pi^{3}}{45} \frac{r_{+}^{2}}{\left[g'(r_{+})\right]^{3}} \Delta_{h} \left[ \frac{g'(r_{+})}{4\pi}+C_{\rm res} \right]^{3}.\]
Using \[A_{+}=4\pi r_{+}^{2},\] and factoring out the ordinary temperature \(g'(r_{+})/(4\pi)\), we get \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+\frac{4\pi C_{\rm res}}{g'(r_{+})} \right]^{3}.\] Since \[g'(r_{+})=\frac{\Xi(r_{+})}{r_{+}},\] the \(f(Q)\)-adapted expression is \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right]^{3}.\]
If the outer boundary of the thin film is much farther than the brick-wall cutoff, then \(H\gg h\), and \[\Delta_{h}\simeq \frac{1}{h^{2}}.\] The entropy reduces to \[S_{\rm D}^{\rm RVB} \simeq \frac{7g_{\rm D}}{2880\pi} \frac{A_{+}}{h^{2}} \left[ 1+\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right]^{3}.\] For a four-component Dirac field, \(g_{\rm D}=4\), so \[S_{\rm D}^{\rm RVB} \simeq \frac{7}{720\pi} \frac{A_{+}}{h^{2}} \left[ 1+\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right]^{3}.\]
For small \(C_{\rm res}\), the first-order expansion is \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{12\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right] +O(C_{\rm res}^{2}).\]
Now consider the quadratic deformation \[f(Q)=Q+\alpha Q^{2}.\] Following the simple \(f(Q)\)-deformed black hole ansatz, \[g(r)=1-\frac{2M}{r}+\alpha r^{2}.\] Thus \[\psi_{Q}(r)=\alpha r^{2}.\] The function \(\Xi(r)\) becomes \[\Xi(r) = 1+\psi_{Q}(r)+r\psi_{Q}'(r).\] Since \[\psi_{Q}'(r)=2\alpha r,\] we obtain \[\Xi(r)=1+3\alpha r^{2}.\] The horizon mass relation is \[M(r_{+}) = \frac{r_{+}}{2} \left( 1+\alpha r_{+}^{2} \right).\] The derivative of the metric at the horizon is \[g'(r_{+}) = \frac{1+3\alpha r_{+}^{2}}{r_{+}}.\] The RVB-corrected temperature is \[T_{\rm RVB}(r_{+}) = \frac{1+3\alpha r_{+}^{2}}{4\pi r_{+}} + C_{\rm res}.\]
The geometric entropy follows from \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = \frac{2\pi r_{+}(1+3\alpha r_{+}^{2})}{1+3\alpha r_{+}^{2}+4\pi C_{\rm res}r_{+}}.\] For small \(C_{\rm res}\), this becomes \[\frac{\mathrm{d}S_{\rm geo}}{\mathrm{d}r_{+}} = 2\pi r_{+} - 8\pi^{2}C_{\rm res} \frac{r_{+}^{2}}{1+3\alpha r_{+}^{2}} + O(C_{\rm res}^{2}).\] Integrating term by term gives \[S_{\rm geo}^{(\alpha)}(r_{+}) = \pi r_{+}^{2} - 8\pi^{2}C_{\rm res} \int^{r_{+}} \frac{u^{2}}{1+3\alpha u^{2}}\,\mathrm{d}u + O(C_{\rm res}^{2}) + S_{0}.\] The integral is evaluated by writing \[\frac{u^{2}}{1+3\alpha u^{2}} = \frac{1}{3\alpha} \left[ 1-\frac{1}{1+3\alpha u^{2}} \right].\] Therefore \[\int \frac{u^{2}}{1+3\alpha u^{2}}\,\mathrm{d}u = \frac{u}{3\alpha} - \frac{1}{3\alpha\sqrt{3\alpha}} \arctan\left(\sqrt{3\alpha}\,u\right).\] Choosing \(S_{0}=0\), we obtain \[S_{\rm geo}^{(\alpha)}(r_{+}) = \pi r_{+}^{2} - \frac{8\pi^{2}C_{\rm res}}{3\alpha}r_{+} + \frac{8\pi^{2}C_{\rm res}}{3\alpha\sqrt{3\alpha}} \arctan\left(\sqrt{3\alpha}\,r_{+}\right) + O(C_{\rm res}^{2}).\]
The Dirac-field entropy in the same quadratic model is obtained by substituting \[\Xi(r_{+})=1+3\alpha r_{+}^{2}\] into the general formula. Thus \[S_{\rm D,\alpha}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{4\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right]^{3}.\] For \(g_{\rm D}=4\), this becomes \[S_{\rm D,\alpha}^{\rm RVB} = \frac{7}{720\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{4\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right]^{3}.\] For \(H\gg h\), we have \[S_{\rm D,\alpha}^{\rm RVB} \simeq \frac{7}{720\pi} \frac{A_{+}}{h^{2}} \left[ 1+ \frac{4\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right]^{3}.\]
The first-order residue expansion is \[S_{\rm D,\alpha}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{12\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right] + O(C_{\rm res}^{2}).\] If \(\alpha r_{+}^{2}\ll 1\), then \[\frac{1}{1+3\alpha r_{+}^{2}} = 1-3\alpha r_{+}^{2} + O(\alpha^{2}).\] Thus \[S_{\rm D,\alpha}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ 12\pi C_{\rm res}r_{+} - 36\pi\alpha C_{\rm res}r_{+}^{3} \right] + O(C_{\rm res}^{2},\alpha^{2}).\]
In the test-field approximation, the total entropy can be written as the sum of the geometric entropy and the Dirac thermal-atmosphere entropy: \[S_{\rm total}^{\rm RVB} = S_{\rm geo}^{\rm RVB} + S_{\rm D}^{\rm RVB}.\] For the quadratic model, this gives \[S_{\rm total,\alpha}^{\rm RVB} = S_{\rm geo}^{(\alpha)} + S_{\rm D,\alpha}^{\rm RVB}.\] Substituting the explicit first-order expressions, we obtain \[\begin{align} S_{\rm total,\alpha}^{\rm RVB} = & \pi r_{+}^{2} - \frac{8\pi^{2}C_{\rm res}}{3\alpha}r_{+} + \frac{8\pi^{2}C_{\rm res}}{3\alpha\sqrt{3\alpha}} \arctan\left(\sqrt{3\alpha}\,r_{+}\right) \\ & + \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{12\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right] + O(C_{\rm res}^{2}). \end{align}\]
The first line is the residue-corrected geometric entropy reconstructed from the first law. The second line is the regulated Dirac-field entropy. The latter is cutoff dependent and should be interpreted as matter-field entanglement or thermal-atmosphere entropy. As usual in brick-wall calculations, this divergent contribution can be absorbed into a renormalization of the gravitational coupling or fixed by a physical cutoff.
First, when the residue correction vanishes, \[C_{\rm res}\rightarrow 0,\] the geometric entropy becomes \[S_{\rm geo}\rightarrow \frac{A_{+}}{4}.\] The Dirac-field entropy becomes \[S_{\rm D} \rightarrow \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h}.\] For \(g_{\rm D}=4\) and \(H\gg h\), this is \[S_{\rm D} \rightarrow \frac{7}{720\pi} \frac{A_{+}}{h^{2}}.\] This is the expected fermionic area-proportional brick-wall result.
Second, in the Schwarzschild limit, \[\alpha\rightarrow 0,\] we have \[\Xi(r_{+})\rightarrow 1.\] Hence \[S_{\rm D}^{\rm RVB} \rightarrow \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left( 1+4\pi C_{\rm res}r_{+} \right)^{3}.\] For small \(C_{\rm res}\), this becomes \[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left( 1+12\pi C_{\rm res}r_{+} \right) + O(C_{\rm res}^{2}).\]
Third, the residue correction modifies the Dirac entropy through the ratio \[\frac{T_{\rm RVB}}{T_{0}} = 1+\frac{4\pi C_{\rm res}}{g'(r_{+})} = 1+\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})}.\] Therefore \[S_{\rm D}^{\rm RVB} = S_{\rm D}^{(0)} \left( \frac{T_{\rm RVB}}{T_{0}} \right)^{3}.\] This compact formula explains why the RVB residue enters the Dirac entropy cubically.
We have derived the Dirac-field entropy of a static spherically symmetric black hole in \(f(Q)\) gravity using the RVB residue-corrected temperature and the thin-film state-counting method. The leading WKB equation of the Dirac field reduces to the Hamilton–Jacobi equation, so the near-horizon pole structure of the radial momentum controls the mode density. After introducing a proper brick-wall cutoff, the Dirac-field entropy is
\[S_{\rm D}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+\frac{4\pi C_{\rm res}r_{+}}{\Xi(r_{+})} \right]^{3}.\]
For the quadratic model \(f(Q)=Q+\alpha Q^{2}\), this becomes
\[S_{\rm D,\alpha}^{\rm RVB} = \frac{7g_{\rm D}}{2880\pi} A_{+}\,\Delta_{h} \left[ 1+ \frac{4\pi C_{\rm res}r_{+}}{1+3\alpha r_{+}^{2}} \right]^{3}.\]
Thus the Dirac-field entropy remains area proportional after regularization, while the RVB residue modifies the coefficient through a temperature-renormalization factor. This provides a direct bridge between the residue-based \(f(Q)\) black hole thermodynamics and the traditional Dirac thin-film entropy method.