# Orbital tomography of hybridized and dispersing molecular overlayers

Thomas Ules  Daniel Lüftner  Eva Maria Reinisch  Georg Koller  Peter Puschnig
<EMAIL_ADDRESS>Michael G. Ramsey<EMAIL_ADDRESS>University of Graz, Institute of Physics, NAWI Graz, Universitätsplatz 5, 8010
Graz, Austria

###### Abstract

With angle resolved photoemission experiments and _ab-initio_ electronic
structure calculations, the pentacene monolayers on Ag(110) and Cu(110) are
compared and contrasted allowing the molecular orientation and an unambiguous
assignment of emissions to specific orbitals to be made. On Ag(110), the
orbitals remain essentially isolated-molecule like, while strong substrate-
enhanced dispersion and orbital modification are observed upon adsorption on
Cu(110). We show how the photoemission intensity of extended systems can be
simulated and that it behaves essentially like that of the isolated molecule
modulated by the band dispersion due to intermolecular interactions.

###### pacs:

33.60.+q,31.15.ae,73.20.-r,68.43.-h

##  I Introduction

The effects arising upon adsorption of $\pi$ conjugated molecules on metal
substrates are of great interest due to their importance in organic
electronics. The ability to identify electronic states with specific molecular
orbitals and to determine their energy ordering is vital to the understanding
of overlayer/substrate systems. This can be particularly difficult for strong
chemisorptive interactions where significant broadening or level splitting
occurs or when intermolecular dispersing bands are formed. In the last decade
_ab-initio_ electronic structure calculations, particulary within density
functional theory (DFT), have become almost indispensable in the
interpretation of experimental results. However, it is becoming recognized
that DFT results can be misleading owing to approximations for exchange-
correlation effects which may severely affect predicted adsorption geometries
and/or the electronic structure of organic/metal interfaces Puschnig  _et al._
(2011). The pentacene/Cu interface is such a case. Despite numerous
experimental and theoretical studies, there is as yet no consensus on the
orbital assignment. Here we demonstrate how the angle resolved photoemission
(ARUPS) technique, that is becoming known as orbital tomography Ziroff  _et
al._ (2010); Puschnig  _et al._ (2011); Stadtmüller  _et al._ (2012); Wießner
_et al._ (2012); Reinisch  _et al._ (2014); Lüftner  _et al._ (2014a); Feyer
_et al._ (2014), can provide a definitive assignment of the emissions even for
strongly interacting extended two-dimensional (2D) systems and give insight
into the nature of dispersion and hybridization.

The pentacene (5A) monolayer on Cu(110) has been extensively studied. From
their ARUPS investigations Seki et al. Yamane  _et al._ (2007, 2008) concluded
that selection rule arguments cannot explain the photoemission behavior and
suggested strong hybridization with the substrate could be modifying the
orbital symmetry Yamane  _et al._ (2008). They also suggested the appearance
of dispersing, interface-induced states arising from substrate interactions
Yamane  _et al._ (2008). Ferretti et al. Baldacchini  _et al._ (2007);
Ferretti  _et al._ (2007) introduced the possibility that the mixing of 5A
molecular orbitals with the Cu substrate leads to ”interaction states
localized at the interface”, where their calculations suggest partial
occupation of the lowest unoccupied molecular orbital (LUMO) whose symmetry
might not be directly related to the original molecular state. On the same
Cu(119) vicinal surface, scanning tunneling microscopy (STM) images also
suggest partial occupation of the LUMO and dispersive electronic states
associated with a perturbed electron charge density distribution. Annese  _et
al._ (2008) Also a recent combined DFT and ARUPS study of the 5A/Cu(110)
system concludes that only partial LUMO occupation takes place Müller  _et
al._ (2012).

The simple relationship between the angular distribution of the photoemission
current and the Fourier transform of the emitting molecular orbital has been
shown to be reasonable for a number of molecular adsorbate systems on various
noble metal surfaces Puschnig  _et al._ (2009); Ziroff  _et al._ (2010);
Puschnig  _et al._ (2011); Lüftner  _et al._ (2014a); Feyer  _et al._ (2014).
This allows molecular orientations to be determined Puschnig  _et al._ (2009);
Feyer  _et al._ (2014), molecular orbital energy ordering to be deduced
Puschnig  _et al._ (2011); Willenbockel  _et al._ (2013) and even the
reconstruction of the molecular orbitals in real space Lüftner  _et al._
(2014a). For these systems, the angular-dependent emissions can essentially be
accounted for by the photoemission from isolated molecules, thus
intermolecular orbital overlap plays a minor role. However, in extended 2D
overlayer systems, a description in terms of isolated molecular orbitals is no
longer strictly appropriate. Here, with the comparison between pentacene
monolayers on Ag(110) and Cu(110), we show how orbital tomography and the
Fourier transform description can be applied to extended systems with strongly
dispersing emissions. In so doing, we provide a definitive description of the
5A/Cu ARUPS and show that the LUMO is in fact fully occupied and displays a
substrate-induced dispersion which is significantly larger than that reported
for similar organic overlayer systems Wießner  _et al._ (2013a, b). Moreover,
unlike all previous tomography studies we here show evidence for a
modification of orbital shape on adsorption.

##  II Experimental Details

Photoemission experiments were performed at BESSYII using a toroidal electron-
energy analyzer described previously Broekman  _et al._ (2005) which was
attached to the beamline U125/2-SGM of the synchrotron radiation facility
BESSY II, Helmholtz-Zentrum-Berlin. Photon energies of 30 and 35 eV and an
incidence angle of $\chi$ = 40∘ with respect to the surface normal were used.
The polarization direction is in the specular plane, which also contains the
photoelectron trajectory measured. Emitted photoelectrons are recorded
simultaneously with polar angles $\theta$ of $-80^{\circ}$ to $+80^{\circ}$
with respect to the surface normal in an energy window of 1 eV. The energy
window is divided into 40 individual slices which goes well below the
analyzer’s intrinsic energy resolution of 150 meV. Note, however, for the
presented momentum maps we take only the data from the side of the emission
direction that is pointing in the orientation of the electric field vector of
the incident photons, this is referred to as the positive side (see Fig. 1),
i.e., $\theta$ = 0∘ to +80∘. This maximizes the polarization factor appearing
in the photocurrent cross section (see Eq. 4). It is noted that the molecular
features are enhanced relative to the substrate emissions on the positive side
(compare Fig. 4c). To obtain the full ( kx, ky ) range for the presented
momentum maps at constant binding energy, azimuthal scans are made by rotating
the sample around the surface normal in 1∘ steps for an azimuthal angle range
$>180^{\circ}$ and then imposing the substrate’s twofold symmetry to obtain
the full 360∘. The angular emission data are then converted to parallel
momentum components kx and ky using the relation

$\displaystyle
k_{x}=\sqrt{\frac{2m_{e}}{\hbar^{2}}E_{\mathrm{kin}}}\sin\theta\cos\phi$ (1)

$\displaystyle
k_{y}=\sqrt{\frac{2m_{e}}{\hbar^{2}}E_{\mathrm{kin}}}\sin\theta\sin\phi,$ (2)

to create the momentum maps. The Cu(110) and Ag(110) substrates were prepared
in the conventional way by a sequence of sputter-annealing cycles. 5A
molecules were evaporated from an effusion cell onto the surfaces at room
temperature with the amount monitored by a microbalance. The resulting
monolayer LEED structures are $\left(\begin{array}[]{cc}3&-1\\\
1&4\end{array}\right)$ for pentacene on Ag(110) and
$\left(\begin{array}[]{cc}6.5&-1\\\ -0.5&2\end{array}\right)$ for pentacene on
Cu(110). In both cases mild annealing improved the order. On Cu(110) the
monolayer was annealed at 200∘C which is above the 5A sublimation temperature
whereas on Ag(110) the molecule substrate bond is weaker and the annealing
temperature must not exceed 140∘C so as not to desorb from the monolayer.

Figure 1:  Sketch of an ARUPS experiment declaring the geometry of the
experimental setup. The incoming photon with the energy $h\nu$, the incidence
angle $\chi$ and vector potential A excites an electron from the initial state
$\psi_{i}$ to the final state $\psi_{f}$. This final state is characterized by
the kinetic energy $E_{\mathrm{kin}}$ and the momentum vector k and the
outgoing photoelectron is detected as a function of $E_{\mathrm{kin}}$ and
emission direction, defined by the polar angle $\theta$ and the azimuthal
angle $\phi$. Where forward emissions are denoted for by a positive parallel
momentum and backward by a negative, respectively.

##  III Computational Details

###  III.1 Density functional calculations

All theoretical results presented here are obtained within the framework of
density functional theory (DFT) using the VASP code Kresse and Hafner (1993);
Kresse and Joubert (1999). We have performed three types of calculations:
firstly for the isolated pentacene molecule, secondly for a two-dimensional,
extended free-standing layer of pentacene, and thirdly for monolayers of
pentacene adsorbed on Ag(110) and Cu(110) surfaces.

The isolated molecule calculations were performed using a supercell with a
minimum of 15 Å vacuum between pentacene’s periodic replica. We use the
generalized gradient approximation (GGA) Perdew  _et al._ (1996) for exchange-
correlation effects, and the projector augmented wave (PAW) method Blöchl
(1994). The simulated momentum maps of the pentacene HOMO and LUMO orbitals
shown in Fig. 2 are obtained as Fourier transforms of the respective Kohn-Sham
orbitals as described previously Puschnig  _et al._ (2009).

Electronic structure calculations for the freestanding monolayer of 5A have
been carried out using the repeated slab approach with a vacuum layer of 20 Å
between adjacent layers. The generalized gradient approximation (GGA) Perdew
_et al._ (1996) is used for exchange-correlation effects, and the projector
augmented waves (PAW) Blöchl (1994) approach was used allowing for a
relatively low kinetic energy cut-off of about 400 eV. We use a Monkhorst-Pack
$6\times 12\times 1$ grid of $k$-points Monkhorst and Pack (1976), and a
first-order Methfessel-Paxton smearing of 0.1 eV Methfessel and Paxton (1989).
The simulation of momentum maps for extended systems will be described in the
subsequent section.

Finally, we have also performed DFT calculations of pentacene monolayers
adsorbed on Ag(110) and Cu(110). The substrate is taken into account within
the repeated slab approach by using five metallic layers with an additional
vacuum layer of 15 Å between slabs. To avoid spurious electrical fields, a
dipole layer is inserted in the vacuum region Neugebauer and Scheffler (1992).
In the case of 5A/Ag(110), we have taken into account the experimental LEED
structure $\left(\begin{array}[]{cc}3&-1\\\ 1&4\end{array}\right)$ mentioned
above and relaxed the atomic positions of the molecule and the first metallic
layer considering van-der-Waals interactions by employing the empirical
correction scheme according to Grimme Grimme (2006). Exchange correlation
effects were treated either within the GGA Perdew  _et al._ (1996) or within
the hybrid functional HSE Heyd  _et al._ (2006) with k-point meshes of
$9\times 6\times 1$ and $6\times 4\times 1$, respectively, and a first-order
Methfessel-Paxton smearing of 0.1 eV. For 5A/Cu(110), we have chosen a
commensurate structure $\left(\begin{array}[]{cc}6&0\\\ 0&2\end{array}\right)$
and an adsorption site similar to a previous study Müller  _et al._ (2012).
Also, here the electronic structure is calculated within the GGA and HSE using
k-point meshes of $9\times 12\times 1$ and $4\times 6\times 1$, respectively,
and using a first-order Methfessel-Paxton smearing of 0.1 eV.

###  III.2 Simulation of ARUPS maps

The Kohn-Sham energies $\varepsilon_{n\mathbf{q}}$ and orbitals
$\psi_{n\mathbf{q}}$ of the relaxed structures serve as input for the
subsequent simulation of ARPES intensity maps within the one-step model of
photoemission Feibelman and Eastman (1974). Here, the angle-resolved
photoemission intensity $I(\theta,\phi;E_{\mathrm{kin}},\omega)$ is a function
of the azimuthal and polar angles $\theta$ and $\phi$, respectively, the
kinetic energy of the emitted electron $E_{\mathrm{kin}}$, and the energy
$\omega$ and polarization $\mathbf{A}$ of the incoming photon:

$\displaystyle I(\theta,\phi,E_{\mathrm{kin}};\omega,\mathbf{A})$
$\displaystyle\approx$
$\displaystyle\sum_{n}^{\textrm{occ}}\sum_{\mathbf{q}}^{\textrm{BZ}}\left|\langle\psi_{f}(\theta,\phi;E_{\mathrm{kin}})|\mathbf{A}\cdot\mathbf{p}|\psi_{n\mathbf{q}}\rangle\right|^{2}$
(3) $\displaystyle\times$
$\displaystyle\delta(\varepsilon_{n\mathbf{q}}+\Phi+E_{\mathrm{kin}}-\omega),$

This formula can be viewed as a Fermi’s golden rule expression, in which the
photocurrent $I$ results from a summation over all occupied initial states
$\psi_{n\mathbf{q}}$, characterized by the band index $n$ and Bloch vector
$\mathbf{q}$, weighted by the transition probability between the initial state
and a final state. For the transition operator $\mathbf{A}\cdot\mathbf{p}$,
the dipole approximation is assumed, where $\bf{p}$ and $\bf{A}$ are the
momentum operator and the vector potential connected to the incoming photon.
The $\delta$ function ensures energy conservation where $\Phi$ denotes the
work function.

We further approximate the final state $\psi_{f}$ by a plane wave Gadzuk
(1974). As outlined in more detail in a previous paper Puschnig  _et al._
(2009), and also noted earlier Shirley  _et al._ (1995); Mugarza  _et al._
(2003), this approximation allows us to greatly simplify the evaluation of the
matrix element appearing in Eq. 3. If we denote the wave vector of the of the
final, free-electron state by $\mathbf{k}$, thus
$E_{\mathrm{kin}}=\frac{\hbar^{2}}{2m}k^{2}$, Eq. 3 simplifies to

$\displaystyle I(k_{x},k_{y},E_{\mathrm{kin}};\omega,\mathbf{A})$
$\displaystyle\approx$
$\displaystyle\sum_{n}^{\textrm{occ}}\sum_{\mathbf{q}}^{\textrm{BZ}}\left|\mathbf{A}\cdot\mathbf{k}\right|^{2}\left|\langle
e^{i\mathbf{kr}}|\psi_{n\mathbf{q}}\rangle\right|^{2}$ (4)
$\displaystyle\times$
$\displaystyle\delta(\varepsilon_{n\mathbf{q}}+\Phi+E_{\mathrm{kin}}-\omega).$

We obtain the simple result that the matrix element from a given initial state
$n\mathbf{q}$ is proportional to the square modulus of the Fourier transform
of the initial state wave function $\psi_{n\mathbf{q}}$ modulated by the
weakly angle-dependent factor $|\bf{A}\cdot\bf{k}|^{2}$ which depends on the
angle between the polarization vector $\bf{A}$ of the incoming photon and the
direction $\mathbf{k}$ of the emitted electron.

Note Eq. 4 can be applied to single molecules as well as to extended states
such as organic layer adsorbed on metallic surfaces. In the former case the
summation of the Brillouin zone (BZ) reduces to just one point, the $\Gamma$
point, while for the latter situation, the possible dispersing bands are taken
into account by the band structure $\varepsilon_{n\mathbf{q}}$ and Bloch
states $\psi_{n\mathbf{q}}$ and an appropriate sampling of the Brillouin zone.

##  IV Results

###  IV.1 Momentum maps

Figure 2:  Panels a) and e) show LEED images of 5A monolayers on Ag(110) and
Cu(110), respectively. b) and c) show ARUPS momentum maps of 5A/Ag(110) at
binding energies 0.15 and 1.20 eV (_exp_), respectively, compared to simulated
LUMO and HOMO maps of the isolated 5A molecule (_sim_). Panels f) and g) show
the corresponding data for 5A/Cu(110) at binding energies of 0.80 and 1.50 eV.
Panels d) and h) show structural models of 5A/Ag(110) and 5A/Cu(110),
respectively, as deduced from LEED and the momentum maps. In h) the orbital
structure of the LUMO is overlaid.

Monolayers of 5A on Ag(110) and Cu(110) were characterized by low energy
electron diffraction (LEED) and ARUPS. Fig. 2 shows LEED images of the 5A
monolayers on the two substrates Ag(110) and Cu(110) as well as the momentum
maps taken at the binding energies of the prominent adsorption induced
emissions at 0.15 and 1.2 eV on Ag and 0.8 and 1.5 eV on Cu, respectively.
These maps clearly display the character of the LUMO and HOMO of pentacene as
can be seen from the comparison between the measured (_exp_) and the
corresponding simulated maps computed for the isolated molecule (_sim_). The
simulated LUMO map, for instance, is characterized by a major lobe at (0,1.4)
Å-1 and minor lobes at ($\pm 1.8$, $1.2$) Å-1. Indeed the data on Ag has
recently been used to reconstruct these orbitals in real space in excellent
agreement with calculated orbitals for the isolated molecule Lüftner  _et al._
(2014a). The maps also unambiguously reveal the _orientation_ of the
molecules: on Ag(110) the data implies flat lying 5A oriented parallel to the
[001] azimuth, while on Cu the molecules orient along the [1-10] azimuth as
depicted in Fig. 2d and h, respectively. It should be noted that the structure
and density of the two monolayers are very similar, however, on Ag the
molecules lie across while on Cu they lie parallel to the close packed rows of
the substrate.

Thus, the ARUPS momentum maps of Fig. 2 immediately clarify the 5A/Cu system.
In contrast to speculations in the literature, the symmetry of the orbitals is
not modified significantly on hybridization. This allows the emissions to be
unambiguously assigned. The features at 0.8 and 1.5 eV are the fully occupied
LUMO and HOMO and not the HOMO and HOMO-1 emissions as previously and, in the
light of our results, erroneously assigned on the basis of a comparison
between DFT-calculated and experimental binding energies Müller  _et al._
(2012). On closer inspection of the LUMO map on Cu, however, an internal
structure and a change in the $k$-position of the minor lobes becomes visible
that will be argued to be the result of intermolecular dispersion and changes
in the orbital size, respectively.

###  IV.2 DFT results for adsorbed monolayers

Figure 3:  Density of states of a monolayer of pentacene on Ag(110) projected
onto orbitals of the free pentacene molecule using PBE-GGA (a) and the hybrid
functional HSE06 (b). Projections onto the LUMO (blue), the HOMO (red), the
HOMO-1 (orange) and the HOMO-2 (green) are shown. For comparison, panel (c)
shows experimental ARPES data along the [001] azimuth (black solid line) and
an azimuthal direction 45∘ between [001] and [1-10]. The colored curves are
obtained by $k_{x}k_{y}$-integrating experimental momentum maps identifying
the orbital character of the emissions. Panels (d)–(f) show the corresponding
data for a monolayer of pentacene on Cu(110).

Before analyzing this internal structure of the LUMO in more detail, we
discuss DFT results of the pentacene monolayers on the Ag(110) and Cu(110)
surfaces in comparison with ARPES data. Starting from the relaxed structure of
5A/Ag(110), we have computed the density of states projected onto the orbitals
of a freestanding pentacene layer, both, within GGA and self-consistently
within the hybrid functional HSE Heyd and Scuseria (2004); Heyd  _et al._
(2006). It is known that the incorporation of a fraction of exact exchange in
the hybrid HSE mediates self-interaction errors thereby improving the orbital
energies and thus the description of adsorbate systems Körzdörfer  _et al._
(2009); Lüftner  _et al._ (2014b, c). The results are shown in panels (a) and
(b) of Fig. 3. For comparison, panel (c) shows experimental ARPES data
recorded over a large energy window along the two different azimuths [001]
(black solid line) and in a direction 45∘ between [001] and [1-10] (black
dashed line). Panel (c) also includes the $k_{x}k_{y}$-integrated ARPES data
of the momentum maps that have been recorded over four energy windows of 1 eV
centered around the molecular emissions of the LUMO, HOMO, HOMO-1, and HOMO-2,
respectively. When comparing GGA with HSE results, we observe that the HOMO,
HOMO-1 and HOMO-2 features are shifted to somewhat larger binding energies
when using the hybrid functional while the position of the LUMO orbital
slightly below the Fermi level remains unchanged. Compared to experiment, we
note that the HSE calculation clearly improves energy position of the HOMO-1
and HOMO-2 and also leads to a slightly better agreement of the HOMO position.
Regarding the LUMO, both functionals GGA and HSE indicate partial occupation
in good agreement with experiment.

For the case of Cu(110) (Fig. 3, panel d–f), HSE again yields much better
agreement with experiment. However, the LUMO, of prime importance in this
study, is still in poor agreement with experiment. Its computed position,
both, obtained within GGA and within HSE, is in error compared to experiment
as the LUMO is located at the Fermi level and not fully occupied with only a
slight improvement gained by the hybrid functional calculation Hofmann  _et
al._ (2013). Note that various commensurate structures and adsorption sites
did not change this finding. On the other hand, HSE does considerably improve
the energy position of the deeper lying orbitals, for instance, the HOMO. As a
side-note: It was exactly the wrong energy position of the HOMO in the GGA
result which had motivated Müller et al. to erroneously assign the LUMO
emission to the HOMO. Note that the computed DOS of their paper Müller  _et
al._ (2012) agrees with our GGA-DOS and their experimental ARUPS data is in
line with our experiment.

###  IV.3 1D Band dispersion

Let us now concentrate on the energy dispersion of the LUMO state as suggested
by the internal structure of the momentum maps. Fig. 4 shows the $E$ _vs._ $k$
relation (bandmap) along the long molecular axes of 5A for a) a monolayer an
Ag(110), b) a half monolayer on Cu(110) and c) a complete monolayer on Cu(110)
from the Fermi edge ($E_{F}$) down to 2.2 eV binding energy. On Ag, the
bandmap shows a feature located at the Fermi edge at a $k$ value around 1.35
Å-1 whose intensity gradually tales off with increasing binding energy
indicating that the LUMO is half-filled. Note that in this energy range, there
are no apparent changes in the momentum maps. In Fig. 4b, the half monolayer
of 5A on Cu(110) is investigated. By examining such a low coverage, with no
long-range order apparent from LEED, any effects of intermolecular dispersion
are minimized and only effects due to the molecule-substrate interaction are
expected. The bandmap of the 5A submonolayer on Cu(110) in the [1-10]
direction (Fig. 2b), _i.e._ along the long molecular axes, shows, besides the
Cu $d$ and $sp$ bands, emissions from the major 5A LUMO lobe visible around 1
eV binding energy. Comparing with Ag, the LUMO shifts from $E_{F}$ and being
half-occupied to a binding energy of roughly 1 eV and being fully occupied on
Cu. This is indicative of the strong bonding interaction of 5A with Cu,
nevertheless, the band map shows no obvious sign of intermolecular dispersion.

Figure 4: a) Bandmap of a 5A monolayer on Ag(110) in the [001] azimuth. b) and
c) show band maps in the [1-10] azimuth of a _half_ and a _complete_ monolayer
of 5A on Cu(110), respectively. Vertical lines indicate the Brillouin zone
boundaries (yellow) and the $\Gamma$ points (purple) of the 5A overlayer.

Upon the formation of the well ordered monolayer on Cu, the situation is
changed with new structure appearing in both the band and momentum maps. As
evident in Fig. 4c, the principle molecular induced emission now oscillates in
energy from $E_{F}$ down to 1.2 eV in the range of $k_{[1-10]}$ from 1.0 to
1.9 Å-1 suggestive of strong intermolecular dispersion (see Fig. 5b for an
enlarged view). To emphasize the $E$ vs. $k$ relation of the LUMO along the
long molecular axes direction, the Brillouin zone (BZ) boundaries and $\Gamma$
points of the 5A monolayer structure deduced from LEED are indicated in Fig.
4c. Note that without differentiation or other enhancements of the raw data,
the large LUMO dispersion with the periodicity of the overlayer is clearly
evident in the 4th and 5th zones. The high binding energy side of the band is
seen to be at $\Gamma$ and the low binding energy at the zone boundaries.
Given the LUMO orbital’s topology (Fig. 2h), the band dispersion running up
from $\Gamma$ to the BZ boundary is qualitatively consistent with a chain of
pentacene molecules.

###  IV.4 2D Band dispersion

Having qualitatively understood the one-dimensional dispersion along the long
molecular axes, a more quantitative and, above all, an understanding of the
two-dimensional (2D) LUMO dispersion on Cu and its manifestation in the
photoemission experiment is desired. To this end, we first discuss the
computed 2D LUMO dispersion of the _free-standing layer_ which is plotted in
Fig. 5a with the band energy $E(k_{[001]},k_{[1-10]})$ color-coded: green
indicating the center of the band, and red (blue) the top (bottom) of the
band, respectively, where for reasons of clarity just one of the mirror
domains is included. Note that although our 5A/Cu(110) HSE calculation
including the substrate considerably improves the agreement with experiment
over GGA results, we refrained from using it here for the analysis of the
experimentally observed band dispersion of the LUMO for the following reason.
The experimental overlayer structure is non-commensurate while, in order to
enable a DFT calculation including Cu(110), a commensurate surface unit cell
had to be imposed (similar to a previous study Müller  _et al._ (2012)). This
difference would result in distinctly different periodicities in the measured
and simulated momentum maps preventing a one-to-one comparison. In contrast,
the free-standing layer simulation with the correct structure enables a one-
to-one analysis of the experimentally observed LUMO dispersion.

Figure 5:  a) Color-coded two-dimensional dispersion $E(k_{[001]},k_{[1-10]})$
of the LUMO for a freestanding 5A monolayer with the structure of Fig. 2h.
Superimposed in white are the Brillouin zones of Cu(110) (rectangle) and the
5A monolayer (rhombus), and in black an intensity iso-line of the calculated
isolated 5A LUMO momentum distribution. b) Close-up of the dispersing LUMO
band with red, green and blue lines indicating the energy positions at which
momentum maps are extracted for comparison with the calculated photoemission
momentum distributions in c) top d) middle and e) bottom of the band. Note
that in panels c)–e) simulated momentum maps (_sim_) are compared to
experimental maps (_exp_) as detailed in the text.

Inspection of Fig. 5a shows that the character of the dispersion, i.e., the
$k$-positions of the tops and bottoms of the band, may well be understood by
knowing (i) the respective Brillouin zone, and (ii) the nodal structure of the
respective molecular orbital in conjunction with the intermolecular
arrangement of molecules. The first determines the reciprocal periodicities
and the second the direction of the dispersion, i.e., whether the band would
run up or down from $\Gamma$ to the BZ boundary. The main dispersion is along
the long real space unit cell vector, _i.e._ , roughly along the long
molecular axes. The calculated energy spread of 20 meV for the free-standing
layer is, however, very small and more than an order of magnitude lower than
the experimentally observed one.

In the ARUPS experiment, one does not observe the 2D dispersion of the LUMO-
derived band over the entire $k$-range due to matrix element effects. For
isolated molecules, the selection rules are well-described by the FT of the
isolated molecular orbitals. In analogy, in order to simulate the
photoemission intensity distribution in a quantitative manner, we thus need to
consider the wave functions of the extended 2D molecular system and compute
their Fourier transforms according to the plane-wave final state approximation
Puschnig  _et al._ (2009). Evaluating Eq. 4 for a free-standing layer, we
simulate momentum maps (_sim_) at the top, middle and bottom of LUMO-derived
band and compare them with the experimental maps (_exp_) taken at the 0.15 eV,
0.60 eV and 1.0 eV binding energy in panels c, d and e of Fig. 5,
respectively. These energy positions are also indicated as horizontal lines in
the enlarged band map shown in panel b. Note that in the simulations, we have
taken into account that the monolayer structure consists of two mirror
domains.

The experimental and simulated momentum maps are in remarkably good agreement.
For the major lobe, in going from top to bottom of the band, not only is the
general shape in agreement, one also observes an increase in the extension of
the feature in $k_{[001]}$ and a splitting of the emission in $k_{[1-10]}$.
Although weak there is also agreement in the behavior of the minor lobes. For
instance, the shift in $k_{[1-10]}$ from 1.0 to 1.3 Å-1 is seen in both
experiment and simulation. Also, in the middle of the band the splitting of
the minor lobe seen in the experiment, although difficult to discern, is also
observable in the calculated map. We conclude that the photoemission intensity
of the extended system is essentially that of the isolated molecule modulated
by the intermolecular dispersion. This naturally implies that the practice of
searching for dispersion in extended systems outside the $k$-range expected
for the isolated molecules is questionable. This is illustrated in the
calculated 2D dispersion of Fig. 5a by overlaying the isoline of the computed
ARUPS map of the the isolated 5A LUMO.

While giving extremely good agreement in the three momentum maps shown in Fig.
5c–e, the exact $E$ vs. $k$ relation is naturally problematic when analyzed in
terms of a free-standing layer. Indeed, as experience with other systems has
shown, the substrate may or may not (depending of the relative energy position
of molecular and substrate levels) enhance the dispersion as e.g. for
PTCDA/Ag(110) Wießner  _et al._ (2013a) and NTCDA/Ag(110) or NTCDA/Cu(100)
Wießner  _et al._ (2013b). But in all these cases the _character_ of the
dispersion is left unaltered and therefore the appearance of momentum maps is
also unchanged when going from the free-standing layer to the adsorbed
monolayer. The _magnitude_ of the dispersion is generally influenced by the
molecule-substrate interaction. The case of pentacene/Cu(110) stands out
because the observed dispersion is the largest one measured so far and may
thus be even termed substrate-_induced_ since the band width of the free-
standing layer’s LUMO band is almost negligible.

Indeed, our calculation for 5A/Cu(110) using the commensurate structure does
show enhanced dispersion compared to the free-standing layer without the
substrate, however, not as large as the experimental observation presumably as
a result of the underestimation of the LUMO binding energy.

###  IV.5 Molecule-Metal Hybridization

For the 5A/Cu system, in addition to the structure introduced by dispersion,
there is also a distinct change in the $k$-position of the minor lobes in the
$k_{[001]}$ direction. In Fig. 6, simulated momentum maps of the LUMO for the
gas-phase and for the 5A/Cu(110) adsorbate system are compared to experimental
maps of the 5A LUMO on both surfaces. On inspection of Fig. 6d, one sees that
on Cu, the minor lobes (green crosses) have maxima at $k_{y}=\pm 1.5$ Å-1
which is significantly lower than, both, the simulated values for the isolated
molecule (b) and the experimental value on the Ag surface (e) ($\pm 1.8$ Å-1).
This shift is a consequence of the molecule-copper interaction and can be
interpreted as a $\approx 20$% increase of the lateral orbital size upon
adsorption. It can be rationalized by recalling the reciprocal relationships
between the LUMO shape in real space and the corresponding momentum space
patterns as illustrated in Figs. 6a and b. Here the maxima of the major and
minor lobes are marked as red and green crosses in the momentum map,
respectively, while the corresponding real space dimensions are indicated by
scale bars of length $2\pi/k$. We note that the width of the LUMO orbital
along the short molecular axis is reflected by the $k_{y}$ maximum of the
minor lobe. Thus, unlike all molecular adsorbate systems so far reported
Puschnig  _et al._ (2009); Dauth  _et al._ (2011); Lüftner  _et al._ (2014a);
Wießner  _et al._ (2012), the orbital tomography of the 5A LUMO on Cu(110)
shows a distinct modification of the orbital’s shape from the isolated
molecule.

Figure 6: a) Calculated LUMO of an isolated 5A molecule in real space, scale
bars mark characteristic dimensions as detailed in the text. b) and c) are
simulated momentum maps for a free 5A molecule and the monolayer 5A/Cu(110),
while d) and e) show experimental ARPES maps of 5A/Cu(110) and 5A/Ag(110),
respectively. The red and green crosses indicate the position of the maxima of
the LUMO’s major and minor lobe, the horizontal green line is a guide for the
eye.

Our interpretation of the experimental data in terms of a spatial distortion
of the orbital is further supported by our DFT calculations. Starting from the
orbital energies and wave functions of the HSE calculation for 5A/Cu(110), we
simulate an ARPES intensity map of the LUMO as shown in Fig. 6c. These
simulations indeed indicate a shift of the minor lobe peak (green cross) to
lower $k$-values upon adsorption on Cu(110). We note, however, that this
change is smaller than the experimental observation. The discrepancy is likely
due to a shortcoming of the DFT results regarding the LUMO binding energy (too
low) and width (too narrow) at the HSE-DFT level indicating an underestimation
of the interaction strength of the LUMO with the substrate. The pronounced
difference between the minor lobe peak position of 5A/Cu(110) compared to
5A/Ag(110) can thus be seen as direct evidence for a strong molecule-metal
hybridization for the former.

While hybridization is generally observed by indirect means in effects on
molecular emissions such as energy broadening of orbital emissions Wießner
_et al._ (2013a) or reduced emission intensity in gaps of bulk states
Berkebile  _et al._ (2011), we here provide direct evidence for changes in
orbital shape. Moreover, further signatures of hybridization can be observed
as additional Cu-$sp$-like features appearing with the overlayer periodicity
(Fig. 4). These are best seen on the _negative_ $k$ side where the molecular
features are weak as shown in Fig. 4c. These emissions are not simply the bulk
$sp$ bands scattered by the overlayer since mere scattering would replicate
the Cu $sp$ band over the entire energy range. Instead one observes that they
do not extend up to $E_{F}$ but stop at 1.2 eV binding energy, just below the
onset of the LUMO orbital. We therefore suggest they are due to interfacial Cu
$sp$ bands hybridizing with the 5A LUMO such that states appearing in these
bands in the 0–1 eV binding energy range have changed their character from Cu
$sp$ to that of the molecular LUMO.

##  V Conclusion

In summary, a very strong substrate-induced dispersion of the LUMO orbital of
5A/Cu(110) along the Cu-row direction is traced back to a significant
hybridization between the organic and the metallic states, while the LUMO of
5A/Ag(110) exhibits only a minor intermolecular dispersion. By making use of
the reciprocal relation between the structure of real space orbitals and
features in the momentum maps, we are able to deduce a significant geometrical
modification of the LUMO orbital upon adsorption on Cu(110), while it remains
essentially free-molecule like on Ag(110). We believe that the main difference
between the adsorption behavior of pentacene on these two surfaces arises from
the distinct molecular orientation with respect to the close-packed metal
rows. While on Ag(110) pentacene orients perpendicular to the rows, it aligns
parallel to the metal rows on Cu(110). Because the electronic structure, e.g.
the band dispersion of the metallic $sp$-bands or the location of the surface
state in the surface Brillouin zone, is distinctly different along the rows
$[1\overline{1}0]$ and perpendicular to them $[001]$, also the hybridization
with adsorbed elongated species can be expected to be different.

In conclusion, for the case of pentacene monolayers on Ag(110) and Cu(110)
surfaces, we have demonstrated the power of the orbital-tomography method
using extensive angle-resolved photoemission data for revealing the electronic
structure of such two-dimensional organic-metal interfacial layers. By
generalizing the theoretical description of the photoemission process from
isolated molecule systems to extended two-dimensional systems, we are able to
explain the ”fine-structure” in the experimental momentum maps and obtain a
comprehensive understanding of the band dispersion of organic-metal
interfacial layers.

###### Acknowledgements.

We acknowledge financial support from the Austrian Science Fund (FWF)
P21330-N20 and P23190-N16. We further acknowledge the Helmholtz-Zentrum Berlin
- Electron storage ring BESSY II for provision of synchrotron radiation at
beamline U125/2-SGM.

## References

  * Puschnig  _et al._ (2011) P. Puschnig, E.-M. Reinisch, T. Ules, G. Koller, S. Soubatch, M. Ostler, L. Romaner, F. S. Tautz, C. Ambrosch-Draxl, and M. G. Ramsey, Phys. Rev. B 84, 235427 (2011).
  * Ziroff  _et al._ (2010) J. Ziroff, F. Forster, A. Schöll, P. Puschnig, and F. Reinert, Phys. Rev. Lett. 104, 233004 (2010).
  * Stadtmüller  _et al._ (2012) B. Stadtmüller, M. Willenbockel, E. Reinisch, T. Ules, M. Ostler, F. Bocquet, S. Soubatch, P. Puschnig, G. Koller, M. G. Ramsey, F. S. Tautz, and C. Kumpf, European Physics Letters 100, 26008 (2012).
  * Wießner  _et al._ (2012) M. Wießner, N. Rodriguez-Lastra, J. Ziroff, F. Forster, P. Puschnig, L. Dössel, K. Müllen, A. Schöll, and F. Reinert, New J. Phys. 14, 113008 (2012).
  * Reinisch  _et al._ (2014) E. M. Reinisch, T. Ules, P. Puschnig, S. Berkebile, M. Ostler, T. Seyller, M. G. Ramsey, and G. Koller, New J. Phys. 16, 023011 (2014).
  * Lüftner  _et al._ (2014a) D. Lüftner, T. Ules, E. M. Reinisch, G. Koller, S. Soubatch, F. S. Tautz, M. G. Ramsey, and P. Puschnig, Proc. Nat. Acad. Sci. U. S. A. 111, 605 (2014a).
  * Feyer  _et al._ (2014) V. Feyer, M. Graus, P. Nigge, M. Wießner, R. Acres, C. Wiemann, C. Schneider, A. Schöll, and F. Reinert, Surface Science 621, 64 (2014).
  * Yamane  _et al._ (2007) H. Yamane, D. Yoshimura, E. Kawabe, R. Sumii, K. Kanai, Y. Ouchi, N. Ueno, and K. Seki, Phys. Rev. B 76, 165436 (2007).
  * Yamane  _et al._ (2008) H. Yamane, E. Kawabe, D. Yoshimura, R. Sumii, K. Kanai, Y. Ouchi, N. Ueno, and K. Seki, phys. stat. sol. (b) 245, 793 (2008).
  * Baldacchini  _et al._ (2007) C. Baldacchini, C. Mariani, M. G. Betti, I. Vobornik, J. Fujii, E. Annese, G. Rossi, A. Ferretti, A. Calzolari, R. Di Felice, A. Ruini, and E. Molinari, Phys. Rev. B 76, 245430 (2007).
  * Ferretti  _et al._ (2007) A. Ferretti, C. Baldacchini, A. Calzolari, R. Di Felice, A. Ruini, E. Molinari, and M. G. Betti, Phys. Rev. Lett. 99, 046802 (2007).
  * Annese  _et al._ (2008) E. Annese, J. Fujii, C. Baldacchini, B. Zhou, C. E. Viol, I. Vobornik, M. G. Betti, and G. Rossi, Phys. Rev. B 77, 205417 (2008).
  * Müller  _et al._ (2012) K. Müller, A. P. Seitsonen, T. Brugger, J. Westover, T. Greber, T. Jung, and A. Kara, J. Phys. Chem. C 116, 23465 (2012).
  * Puschnig  _et al._ (2009) P. Puschnig, S. Berkebile, A. J. Fleming, G. Koller, K. Emtsev, T. Seyller, J. D. Riley, C. Ambrosch-Draxl, F. P. Netzer, and M. G. Ramsey, Science 326, 702 (2009).
  * Willenbockel  _et al._ (2013) M. Willenbockel, B. Stadtmüller, K. Schönauer, F. Bocquet, D. Lüftner, E. M. Reinisch, T. Ules, G. Koller, C. Kumpf, S. Soubatch, P. Puschnig, M. G. Ramsey, and F. S. Tautz, New J. Phys. 15, 033017 (2013).
  * Wießner  _et al._ (2013a) M. Wießner, J. Ziroff, F. Forster, M. Arita, K. Shimada, P. Puschnig, A. Schöll, and F. Reinert, Nature Communications 4, 1514 (2013a).
  * Wießner  _et al._ (2013b) M. Wießner, J. Kübert, V. Feyer, P. Puschnig, A. Schöll, and F. Reinert, Phys. Rev. B 88, 075437 (2013b).
  * Broekman  _et al._ (2005) L. Broekman, A. Tadich, E. Huwald, J. Riley, R. Leckey, T. Seyller, K. Emtsev, and L. Ley, J. Electron. Spectrosc. Relat. Phenom. 144-147, 1001 (2005).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Methfessel and Paxton (1989) M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989).
  * Neugebauer and Scheffler (1992) J. Neugebauer and M. Scheffler, Phys. Rev. B 46, 16067 (1992).
  * Grimme (2006) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
  * Heyd  _et al._ (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124, 219906 (2006).
  * Feibelman and Eastman (1974) P. J. Feibelman and D. E. Eastman, Phys. Rev. B 10, 4932 (1974).
  * Gadzuk (1974) J. W. Gadzuk, Phys. Rev. B 10, 5030 (1974).
  * Shirley  _et al._ (1995) E. L. Shirley, L. J. Terminello, A. Santoni, and F. J. Himpsel, Phys. Rev. B 51, 13614 (1995).
  * Mugarza  _et al._ (2003) A. Mugarza, J. E. Ortega, F. J. Himpsel, and F. J. García de Abajo, Phys. Rev. B 67, 081404 (2003).
  * Heyd and Scuseria (2004) J. Heyd and G. E. Scuseria, J. Chem. Phys. 121, 1187 (2004).
  * Körzdörfer  _et al._ (2009) T. Körzdörfer, S. Kümmel, N. Marom, and L. Kronik, Phys. Rev. B 79, 201205(R) (2009).
  * Lüftner  _et al._ (2014b) D. Lüftner, M. Milko, S. Huppmann, M. Scholz, N. Ngyuen, M. Wießner, A. Schöll, F. Reinert, and P. Puschnig, J. Electron Spectrosc. Relat. Phenom. available online (2014b), 10.1016/j.elspec.2014.06.002.
  * Lüftner  _et al._ (2014c) D. Lüftner, S. Refaely-Abramson, M. Pachler, R. Resel, M. G. Ramsey, L. Kronik, and P. Puschnig, Phys. Rev. B 90, 075204 (2014c).
  * Hofmann  _et al._ (2013) O. T. Hofmann, V. Atalla, N. Moll, P. Rinke, and M. Scheffler, New Journal of Physics 15, 123028 (2013).
  * Dauth  _et al._ (2011) M. Dauth, T. Körzdörfer, S. Kümmel, J. Ziroff, M. Wiessner, A. Schöll, F. Reinert, M. Arita, and K. Shimada, Phys. Rev. Lett. 107, 193002 (2011).
  * Berkebile  _et al._ (2011) S. Berkebile, T. Ules, P. Puschnig, L. Romaner, G. Koller, A. J. Fleming, K. Emtsev, T. Seyller, C. Ambrosch-Draxl, F. P. Netzer, and M. G. Ramsey, Phys. Chem. Chem. Phys. 13, 3604 (2011).


# Near-Infrared and Low-Rank Adaptation of
Vision Transformers in Remote Sensing

Irem Ulku1, O. Ozgur Tanriover1, Erdem Akagündüz2,
1Dept. of Computer Science, Ankara University, Ankara, Türkiye
2Dept. of Modeling and Simulation, Graduate School of Informatics, METU,
Ankara, Türkiye
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Plant health can be monitored dynamically using multispectral sensors that
measure Near-Infrared reflectance (NIR). Despite this potential, obtaining and
annotating high-resolution NIR images poses a significant challenge for
training deep neural networks. Typically, large networks pre-trained on the
RGB domain are utilized to fine-tune infrared images. This practice introduces
a domain shift issue because of the differing visual traits between RGB and
NIR images. As an alternative to fine-tuning, a method called low-rank
adaptation (LoRA) enables more efficient training by optimizing rank-
decomposition matrices while keeping the original network weights frozen.
However, existing parameter-efficient adaptation strategies for remote sensing
images focus on RGB images and overlook domain shift issues in the NIR domain.
Therefore, this study investigates the potential benefits of using vision
transformer (ViT) backbones pre-trained in the RGB domain, with low-rank
adaptation for downstream tasks in the NIR domain. Extensive experiments
demonstrate that employing LoRA with pre-trained ViT backbones yields the best
performance for downstream tasks applied to NIR images.

_Keywords_ low-rank adaptation, ViT, infrared images, semantic segmentation.

##  1 Introduction

The near-infrared (NIR) spectrum offer the most valuable insights into the
health of vegetation on Earth [7]. For example, when a plant is attacked by a
pest or under stress, there’s a noticeable change in the reflectance of NIR
bands due to alterations in the internal leaf structure. NIR reflectance can
uniquely identify the internal leaf structure and biochemical composition,
while also distinguishing healthy plants. Therefore, multispectral sensors
covering NIR bands have the capability to automatically monitor vegetation
health [1]. However, NIR sensors typically possess lower spatial resolutions
and lower signal-to-noise ratios compared to RGB sensors [2]. Consequently,
acquiring and annotating a sufficient number of NIR images with high
resolution and signal-to-noise ratio to effectively train deep neural networks
is a state-of-the-art challenge [19].

A common practice to address this problem is to train deep neural networks
with RGB images and transferring these pre-trained networks with infrared
images [14]. However, the visual characteristics of RGB images exhibit
distinct distributions compared to infrared images, resulting in the widely
recognized ”domain shift problem”. This issue in deep learning is often
addressed by fine-tuning the initial model using inputs from the shifted
domain. Depending on factors such as model complexity, domain shift scale, and
dataset size, solving this problem can prove achievable to a certain extent.
The reader may refer to [2] for a review of solutions to the domain shift
problem in deep learning.

Currently, fine-tuning is realized mainly by foundation models, which require
large-scale data for training. Among various foundation models, transformer-
based architectures [6, 13] offer the best solution, since they achieve lower
inductive bias than convolutional neural networks with the help of positional
embeddings [14]. Although transformer-based models can improve performance,
the benefit comes with the cost of the highest computational burden. This
drawback primarily stems from the requirement of large-scale pretraining on
supervised datasets. As a still ongoing field, transformer-based foundation
models have a high potential for usage in remote sensing-based tasks to
support domain adaptation. To effectively utilize large foundation models for
multispectral remote-sensing-based agriculture tasks, the models must address
difficult challenges, such as achieving high accuracy in real-time speed. The
difficulty mainly comes from the need for a higher number of trainable
parameters to meet the accuracy requirements, consequently slowing down both
training and inference speeds. These constraints render foundation models
impractical for integration into real-time mobile applications. [15].

There is a recent trend towards adaptation methods instead of fine-tuning, as
they offer improved generalization in out-of-domain scenarios due to avoiding
catastrophic forgetting, especially in cases of limited data availability
[11]. These family of methods are referred to as “adaptation methods”, which
tune only few parameters to mitigate overfitting, and can alleviate the
storage burden and increase the performance of traditional fine-tuning
approaches [8]. A recent example of this family, namely Low-Rank Adaptation
(LoRA) [9] approach allows for indirect training of dense layers while
freezing original weights by only optimizing the rank-decomposition matrices
of the changes in these layers. The underlying hypothesis is that change in
weight matrices during adaptation has a low intrinsic rank. Although simple
and computationally efficient, its usage remains confined to the transformer-
based large language models (LLM) from which it originated, with limited
research on vision transformers [12].

Few studies exist in recent literature that analyze the parameter-efficient
adaptation strategies with remote-sensing images [18, 10, 5]. These studies
solely focus on adapting RGB images and do not delve into the challenges
associated with domain shift when adapting to the NIR domain. To address this
gap, this paper proposes a framework to investigate the impact of low-rank
adaptation on the semantic segmentation performance of vision transformers
(ViT) on NIR images. To the best of our knowledge, near-infrared and low-rank
adaptation of vision transformers has not been used for the semantic
segmentation of multispectral images in remote sensing.

Figure 1: The proposed LoRA-based ViT+segmentation acrhictecture.

##  2 Model Architecture

###  2.1 Encoder

Fig.1 shows an embedding layer and an encoder component based on the ViT
architecture. The input image $\mathbf{x}\in\mathbb{R}^{H\times W\times C}$
with dimensions as height $H$, width $W$ and channel $C$, is divided into a
sequence of $\mathbf{x}=\left[x_{1},...,x_{N}\right]\in\mathbb{R}^{N\times
P^{2}\times C}$ patches, each with dimensions $\left(P,P\right)$, resulting in
a total of $N=HW/P^{2}$ patches. Each patch gets flattened into a one-
dimensional vector, and these flattened patches are mapped onto a sequence of
patch embeddings
$\mathbf{x}_{0}=\left[\mathbf{E}x_{1},...,\mathbf{E}x_{N}\right]\in\mathbb{R}^{N\times
D}$ in a D-dimensional latent space using a linear projection layer
$\mathbf{E}\in\mathbb{R}^{D\times\left(P^{2}C\right)}$. To capture spatial
relationships, learnable positional embeddings
$\mathbf{E}_{pos}\in\mathbb{R}^{N\times D}$ are added to the projected patch
embedding sequence as follows:

$\mathbf{z}_{0}=\left[\mathbf{E}x_{1},...,\mathbf{E}x_{N}\right]+\mathbf{E}_{pos}.$
(1)

The contextualized encodings sequence $\mathbf{z}_{L}\in\mathbb{R}^{N\times
D}$, obtained by applying $L$ transformer blocks to the token array
$\mathbf{z}_{0}$, is formed. Several transformer blocks are stacked to obtain
latent embeddings $z_{i}$, where each transformer block consists of a multi-
headed self-attention (MHA) block followed by a multilayer perceptrons (MLP)
block:

${z}^{\prime}_{i}=MHA\left(Norm\left(z_{i-1}\right)\right)+z_{i-1}$ (2)

${z}_{i}=MLP\left(Norm\left({z}^{\prime}_{i-1}\right)\right)+{z}^{\prime}_{i-1},$
(3)

where $i=1,...,L$ and $Norm\left(.\right)$ represents layer normalization.

The self-attention mechanism involves the query
$\mathbf{Q}_{i}\in\mathbb{R}^{N\times d}$, key
$\mathbf{K}_{i}\in\mathbb{R}^{N\times d}$, and value
$\mathbf{V}_{i}\in\mathbb{R}^{N\times d}$ by multiplying the patch embedding
matrix $\mathbf{z}_{0}\in\mathbb{R}^{N\times D}$ with weight matrices:

$\mathbf{Q}_{i}=\mathbf{z}_{0}\mathbf{W}^{Q}_{i},\mathbf{K}_{i}=\mathbf{z}_{0}\mathbf{W}^{K}_{i},\mathbf{V}_{i}=\mathbf{z}_{0}\mathbf{W}^{V}_{i},$
(4)

where $\mathbf{W}^{Q}_{i},\mathbf{W}^{K}_{i},\mathbf{W}^{V}_{i}$ are learnable
weight matrices. When the MHA block contains $L$ attention heads, the self-
attention operation for the $i^{th}$ head is as follows:

$Atten_{i}\left(Q,K,V\right)=softmax\left(\mathbf{Q}_{i}\mathbf{K}_{i}^{T}/\sqrt{D}\right)\mathbf{V}_{i},$
(5)

The MHA mechanism which identifies the regions to give attention using cosine
similarity then projects these concatenated heads through another learnable
parameter, denoted as $\mathbf{W}\in\mathbb{R}^{D\times D}$:

$MHA\left(\mathbf{z}_{L}\right)=\left[head_{1},...,head_{L}\right]\mathbf{W}.$
(6)

These series of transformer mappings align the sequence activation
$\mathbf{z}_{i-1}=\left[z_{i-1,1},...,z_{i-1,N}\right]$, to a contextualized
encoded sequence $\mathbf{z}_{i}=\left[z_{i,1},...,z_{i,N}\right]$ with
enriched semantic information.

###  2.2 LoRA Adaptation

The LoRA layers are applied to the query and value projection layers in each
transformer block based on the hypothesis that the updates of weights have a
low intrinsic rank during adaptation. Assuming the pre-trained weight matrix
of any projection layer as $\mathbf{W}\in\mathbb{R}^{C_{in}\times C_{out}}$,
the updated version incorporates two linear layers, denoted as
$\mathbf{A}\in\mathbb{R}^{r\times C_{out}}$ and
$\mathbf{B}\in\mathbb{R}^{C_{in}\times r}$, while $W$ remains frozen to avoid
receiving gradient updates according to:

$h=\mathbf{W}x+\bigtriangleup\mathbf{W}x=\mathbf{W}x+\mathbf{BA}x.$ (7)

The value of $r$, such that $r\leqslant min\left(C_{in},C_{out}\right)$,
indicates a substantial reduction in the trainable parameters. The resulting
output vectors are summed element-wise after multiplying $\mathbf{W}$ and
$\bigtriangleup\mathbf{W}$ with the same input. It is reasonable to apply the
LoRA to the query and value projection layers to adjust the attention scores
within the MHA block:

$\mathbf{Q}_{i}=\mathbf{z}_{0}\mathbf{W}^{Q}_{i}+\mathbf{z}_{0}\mathbf{B}^{Q}_{i}\mathbf{A}^{Q}_{i}$
(8)
$\mathbf{V}_{i}=\mathbf{z}_{0}\mathbf{W}^{V}_{i}+\mathbf{z}_{0}\mathbf{B}^{V}_{i}\mathbf{A}^{V}_{i}$
(9) $\mathbf{K}_{i}=\mathbf{z}_{0}\mathbf{W}^{K}_{i}$ (10)

where $\mathbf{W}^{Q}_{i}$, $\mathbf{W}^{V}_{i}$ and $\mathbf{W}^{K}_{i}$
represent frozen weights obtained from ViT, while $\mathbf{B}^{Q}_{i}$,
$\mathbf{A}^{Q}_{i}$, $\mathbf{B}^{V}$ and $\mathbf{A}^{V}_{i}$ denote
trainable LoRA parameters.

###  2.3 Segmentation Head

ASPP (Atrous Spatial Pyramid Pooling) [4] is used to design a lightweight
segmentation head that captures multiple scales of context information. ASPP
consists of five parallel branches so that while providing feature mapping
with $1\times 1$ convolution, three atrous convolutions with
$\left[12,24,36\right]$ atrous rates enable multi-scale feature extraction.
ASPP pooling with an adaptive average pooling layer and $1\times 1$
convolution can capture global information, while bilinear interpolation
recovers feature size. As a result of concatenating all these feature maps
into the channel dimension and fusing them with $1\times 1$ convolution, the
output feature map is passed through the Sigmoid function for binary
segmentation after $3\times 3$ and $1\times 1$ convolutions.

##  3 Experimental Results

###  3.1 Image sets

The DSTL image set [17] comprises 25 satellite images covering a region of
1000 m × 1000 m. RGB images have dimensions of 3348 × 3392 pixels with a
spatial resolution of 0.31 m. Multispectral images, on the other hand, have
dimensions of 837 × 848 pixels and a spatial resolution of 1.24 m, spanning
the wavelength range of 400-1040 nm. The DSTL image set focuses on the crop
target class.

The RIT-18 image set [17] contains a single aerial training image with a
spatial resolution of 0.047 m and dimensions of 9393 × 5642 pixels. RIT-18
provides multispectral images at wavelengths of 715-725, 795-805, and 890-910
nm and assesses the tree target class.

Images from both sets are cropped into 224 × 224 patches and split into $72\%$
for training, 20% for testing, and 8% for validation.

Table 1: Backbone details Backbone | Patch | Layers | Hidden | MLP | Heads |
---|---|---|---|---|---|---
B_16 | 16 | 12 | 768 | 3072 | 12 |
B_32 | 32 | 12 | 768 | 3072 | 12 |
L_16 | 16 | 24 | 1024 | 4096 | 16 |
L_32 | 32 | 24 | 1024 | 4096 | 16 |

###  3.2 Experimental Setup

The experiments are implemented on NVIDIA Quadro RTX 5000 GPU using the
Pytorch framework. Training employs the Adam algorithm for 70 epochs with a
mini-batch size of 8 and an initial learning rate of 1e-4. Learning rate is
reduced by 9% in every five iterations. The loss function for all experiments
is binary cross-entropy with logits.

The chosen metrics for performance evaluation are the Jaccard index (IoU) and
the $F_{1}$ score. The Jaccard index calculation is $IoU=(TP/(TP+FP+FN))$,
where $TP$, $FP$, and $FN$ are the number of true positive, false positive,
and false negative pixels, respectively. The $F_{1}$ score is the harmonic
mean of precision (p) and recall (r) metrics. The calculation of precision is
$p=(TP/(TP+FP))$, and that of the recall is $r=(TP/(TP+FN))$.

###  3.3 Results and Discussions

Table 2: Tree semantic segmentation test results in terms of Jaccard Index
(IoU) and F1 score for the different architectures with RIT-18 image set

|  | Trainable Parameters(M) | RGB | NIR
---|---|---|---|---
Backbone | Model | IoU | F1 | IoU | F1
Xception | DeepLabV3 | 58.80 | 0.837 ± 0.318 | 0.859 ± 0.303 | 0.857 ± 0.298 | 0.878 ± 0.282
B_16_imagenet1k | ViT + DeepLabV3 Head | 91.49 | 0.859 ± 0.283 | 0.885 ± 0.260 | 0.881 ± 0.258 | 0.905 ± 0.236
LoraViT + DeepLabV3 Head | 5.84 | 0.861 ± 0.269 | 0.891 ± 0.240 | 0.882 ± 0.246 | 0.909 ± 0.219
L_16_imagenet1k | ViT + DeepLabV3 Head | 310.66 | 0.862 ± 0.280 | 0.888 ± 0.258 | 0.882 ± 0.263 | 0.903 ± 0.246
LoraViT + DeepLabV3 Head | 7.75 | 0.863 ± 0.269 | 0.891 ± 0.241 | 0.884 ± 0.243 | 0.911 ± 0.215
B_32_imagenet1k | ViT + DeepLabV3 Head | 93.15 | 0.823 ± 0.325 | 0.849 ± 0.303 | 0.839 ± 0.305 | 0.866 ± 0.282
LoraViT + DeepLabV3 Head | 5.84 | 0.839 ± 0.297 | 0.869 ± 0.269 | 0.855 ± 0.285 | 0.882 ± 0.260
L_32_imagenet1k | ViT + DeepLabV3 Head | 312.87 | 0.821 ± 0.327 | 0.847 ± 0.307 | 0.849 ± 0.297 | 0.875 ± 0.276
LoraViT + DeepLabV3 Head | 7.75 | 0.838 ± 0.299 | 0.868 ± 0.272 | 0.858 ± 0.279 | 0.886 ± 0.254

Table 3: Crop semantic segmentation test results in terms of Jaccard Index
(IoU) and F1 score for the different architectures with DSTL image set

|  | Trainable Parameters(M) | RGB | NIR
---|---|---|---|---
Backbone | Model | IoU | F1 | IoU | F1
Xception | DeepLabV3 | 58.80 | 0.858 ± 0.322 | 0.869 ± 0.316 | 0.867 ± 0.302 | 0.882 ± 0.291
B_16_imagenet1k | ViT + DeepLabV3 Head | 91.49 | 0.873 ± 0.298 | 0.891 ± 0.285 | 0.884 ± 0.290 | 0.895 ± 0.281
LoraViT + DeepLabV3 Head | 5.84 | 0.876 ± 0.289 | 0.891 ± 0.278 | 0.896 ± 0.271 | 0.908 ± 0.261
L_16_imagenet1k | ViT + DeepLabV3 Head | 310.66 | 0.878 ± 0.296 | 0.890 ± 0.287 | 0.887 ± 0.282 | 0.899 ± 0.272
LoraViT + DeepLabV3 Head | 7.75 | 0.886 ± 0.283 | 0.898 ± 0.273 | 0.900 ± 0.266 | 0.912 ± 0.254
B_32_imagenet1k | ViT + DeepLabV3 Head | 93.15 | 0.860 ± 0.308 | 0.876 ± 0.296 | 0.871 ± 0.297 | 0.886 ± 0.285
LoraViT + DeepLabV3 Head | 5.84 | 0.874 ± 0.298 | 0.887 ± 0.288 | 0.888 ± 0.273 | 0.903 ± 0.257
L_32_imagenet1k | ViT + DeepLabV3 Head | 312.87 | 0.871 ± 0.304 | 0.883 ± 0.293 | 0.876 ± 0.298 | 0.888 ± 0.288
LoraViT + DeepLabV3 Head | 7.75 | 0.871 ± 0.303 | 0.883 ± 0.296 | 0.887 ± 0.276 | 0.902 ± 0.263

Experiments employ the ViT base (B) and ViT large (L) backbones after pre-
training on ImageNet-21k and fine-tuning on ImageNet-1k, where Table 1 shows
their details. In ImageNet, because the training is performed only on RGB
images, the backbone networks represent only the features of the RGB domain.

Table 2 shows semantic segmentation test results for the tree class of the
RIT-18 image set using RGB and NIR images, represented in terms of the Jaccard
index and $F_{1}$ score. When using Lora with the pre-trained ViT-L/16
encoder, this model outperforms all other models for the downstream task in
the NIR domain. The Jaccard index result using ViT-L/16 with Lora for the NIR
domain demonstrates that even though the number of parameters decreases by
approximately 97.5%, it performs 0.2% better than the ViT-L/16 without Lora.
Using Lora with the ViT-L/16 backbone shows an increase in performance by 2.1%
for the NIR domain compared to the RGB domain.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 2: Crop predictions for DSTL image set in NIR domain. Light green
represents a hit, dark green represents a miss, and red represents a false
alarm.(a) Image. (b) Ground-truth mask. (c) DeepLabV3. (d) ViT (L$-$16). (e)
LoraViT (L$-$16).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 3: Tree predictions for RIT-18 image set in NIR domain. Light green
represents a hit, dark green represents a miss, and red represents a false
alarm. (a) Image. (b) Ground-truth mask. (c) DeepLabV3. (d) ViT (L$-$16). (e)
LoraViT (L$-$16).

An attractive feature of the low-rank adaptation is that despite a significant
distributional difference between the RGB and NIR domains, the performances of
backbones with Lora in the NIR domain are superior to RGB counterparts across
all results. Besides that, the Lora decreases the total number of trainable
parameters by approximately 97%.

Additionally, the Lora-based ViT-L/16 backbone competes with the state-of-the-
art DeepLabV3 and outperforms the Jaccard index performance of this model by
approximately 2.7% in the NIR domain. As expected, larger ViT models like
ViT-L/16 demonstrate significantly stronger intermediate representations when
large pretraining datasets like ImageNet-21k are used [16].

As seen in Table 3, adapting the ViT-L/16 backbone with Lora in the NIR domain
shows the best performance also on the DSTL image set for crop objects. The
experimental results show that using ViT-L/16 with Lora on NIR images improves
the Jaccard index performance by 1.3% compared to employing ViT-L/16 without
Lora. Adapting the ViT-L/16 backbone to the NIR domain with Lora shows a 1.4%
improvement in Jaccard index performance compared to its adaptation to the RGB
domain. As shown in Table 2 and Table 3, a patch size of $16$ generally
performs better than those with a patch size of $32$, which conforms with the
findings in the literature that a smaller patch size performs better [3, 6].

Fig. 2 and Fig. 3 present visual comparisons of the various ViT-L/16 backbone
cases for DSTL and RIT18 image sets, respectively. ViT-L/16 backbone with Lora
in the NIR domain shows a lower false alarm rate than other cases.
Additionally, there are fewer missed detections when capturing small target
objects and complex details.

##  4 Conclusion

This study explores the potential benefits of using pre-trained ViT backbones
in the RGB domain with low-rank adaptation for downstream tasks in the NIR
domain. Due to the domain shift problem, directly using ViT backbones or
traditional semantic segmentation models does not yield sufficient results in
the NIR domain. However, extensive experiments in this study demonstrate that
using LoRA with pre-trained ViT backbones provides the best performance for
downstream tasks applied to NIR images. Lora is advantageous, as it typically
contributes to the overall stability of the domain adaptation process.

Given the scarcity of large amounts of labeled data for NIR bands that contain
discriminative information for vegetation, this study is a promising approach
to addressing the problem of insufficient utilization of NIR bands in the
semantic segmentation process. For future studies, it may be valuable to
investigate the application of LoRA to segmentation architectures for the NIR
domain using different segmentation heads.

## References

  * [1] Adeola M Arogoundade, Onisimo Mutanga, John Odindi, and Omosalewa Odebiri.  Leveraging google earth engine to estimate foliar c: N ratio in an african savannah rangeland using sentinel 2 data.  Remote Sensing Applications: Society and Environment, 30:100981, 2023.
  * [2] Masoomeh Aslahishahri, Kevin G Stanley, Hema Duddu, Steve Shirtliffe, Sally Vail, Kirstin Bett, Curtis Pozniak, and Ian Stavness.  From rgb to nir: Predicting of near infrared reflectance from visible spectrum aerial images of crops.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 1312–1322, 2021.
  * [3] Mathilde Caron, Hugo Touvron, Ishan Misra, Hervé Jégou, Julien Mairal, Piotr Bojanowski, and Armand Joulin.  Emerging properties in self-supervised vision transformers.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 9650–9660, 2021.
  * [4] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  IEEE transactions on pattern analysis and machine intelligence, 40(4):834–848, 2017.
  * [5] Zhe Dong, Yanfeng Gu, and Tianzhu Liu.  Upetu: A unified parameter-efficient fine-tuning framework for remote sensing foundation model.  IEEE Transactions on Geoscience and Remote Sensing, 2024.
  * [6] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929, 2020.
  * [7] Alexandros Gkillas, Dimitrios Kosmopoulos, and Kostas Berberidis.  Cost-efficient coupled learning methods for recovering near-infrared information from rgb signals: Application in precision agriculture.  Computers and Electronics in Agriculture, 209:107833, 2023.
  * [8] Haoyu He, Jianfei Cai, Jing Zhang, Dacheng Tao, and Bohan Zhuang.  Sensitivity-aware visual parameter-efficient fine-tuning.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 11825–11835, 2023.
  * [9] Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen.  Lora: Low-rank adaptation of large language models.  arXiv preprint arXiv:2106.09685, 2021.
  * [10] Leiyi Hu, Hongfeng Yu, Wanxuan Lu, Dongshuo Yin, Xian Sun, and Kun Fu.  Airs: Adapter in remote sensing for parameter-efficient transfer learning.  IEEE Transactions on Geoscience and Remote Sensing, 2024.
  * [11] Yaqin Li, Dandan Wang, Cao Yuan, Hao Li, and Jing Hu.  Enhancing agricultural image segmentation with an agricultural segment anything model adapter.  Sensors, 23(18):7884, 2023.
  * [12] Yen-Cheng Liu, Chih-Yao Ma, Junjiao Tian, Zijian He, and Zsolt Kira.  Polyhistor: Parameter-efficient multi-task adaptation for dense vision tasks.  Advances in Neural Information Processing Systems, 35:36889–36901, 2022.
  * [13] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 10012–10022, 2021.
  * [14] Yuchi Ma, Shuo Chen, Stefano Ermon, and David B Lobell.  Transfer learning in environmental remote sensing.  Remote Sensing of Environment, 301:113924, 2024.
  * [15] Bekhzod Olimov, Jeonghong Kim, and Anand Paul.  Ref-net: robust, efficient, and fast network for semantic segmentation applications using devices with limited computational resources.  IEEE Access, 9:15084–15098, 2021.
  * [16] Maithra Raghu, Thomas Unterthiner, Simon Kornblith, Chiyuan Zhang, and Alexey Dosovitskiy.  Do vision transformers see like convolutional neural networks?  Advances in neural information processing systems, 34:12116–12128, 2021.
  * [17] Irem Ulku, Erdem Akagündüz, and Pedram Ghamisi.  Deep semantic segmentation of trees using multispectral images.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 15:7589–7604, 2022.
  * [18] Bowei Xue, Han Cheng, Qingqing Yang, Yi Wang, and Xiaoning He.  Adapting segment anything model to aerial land cover classification with low rank adaptation.  IEEE Geoscience and Remote Sensing Letters, 2024.
  * [19] Ting Zhang, Haijian Shen, Sadaqat ur Rehman, Zhaoying Liu, Yujian Li, and Obaid ur Rehman.  Two-stage domain adaptation for infrared ship target segmentation.  IEEE Transactions on Geoscience and Remote Sensing, 2023.


$\displaystyle\prod_{i=1}^{l}\,\deg\large\big{(}\big{(}\partial_{T(s_{i})}A,\,\partial_{t_{i}}A\big{)},\,\,\big{(}B_{d}(\xi^{*}_{i})\cap\partial\Omega\big{)}\times\big{(}t_{i}^{*}-d,\,t_{i}^{*}+d\big{)},\,\,0\large\big{)}$
$\displaystyle\qquad\quad\,\times\deg\left(\big{(}\partial_{T(s_{l+1})}a,\ldots,\partial_{T(s_{m})}a\big{)},\,\,\prod_{k=l+1}^{m}\big{(}B_{d}(\xi^{*}_{k})\cap\partial\Omega\big{)},\,\,0\right)$
$\displaystyle\neq$ $\displaystyle\,0.$

This implies that for any $p>1$ large enough, there exists
$(\mathbf{s}^{p},\mathbf{t}^{p})\in\Lambda_{d}$ such that
$\nabla_{\left(T(\mathbf{s}),\mathbf{t}\right)}\widehat{F}_{p}\big{(}\mathbf{s}^{p},\mathbf{t}^{p}\big{)}=0$.
In particular,
$\mathbf{s}^{p}=(s^{p}_{1},\ldots,s^{p}_{m})\rightarrow(\xi^{*}_{1},\ldots,\xi^{*}_{m})$
as $p\rightarrow+\infty$, which completes the proof. $\square$

Proof of Theorem 1.2. According to Proposition 5.1, we need to find a critical
point
$\xi^{p}=(\xi^{p}_{1},\ldots,\xi^{p}_{m})\in\Omega^{l}\times(\partial\Omega)^{m-l}$
of $F_{p}$ such that points $\xi^{p}_{1},\ldots,\xi^{p}_{m}$ accumulate to
$\xi_{*}$. Using (1.9), (2.21), (6.1), Lemma 2.1 and the fact that
$a(\xi_{i})G(\xi_{i},\xi_{k})=a(\xi_{k})G(\xi_{k},\xi_{i})$ for all
$i,k=1,\ldots,m$ with $i\neq k$, we conclude that $F_{p}$ becomes

$\displaystyle
F_{p}(\xi)=\frac{e}{2p^{2}}\left\\{8\pi\sum_{i=1}^{l}a(\xi_{i})\left[p-2\log
p-8\pi H(\xi_{i},\xi_{i})-8\pi\sum_{k=1,\,k\neq
i}^{l}G(\xi_{i},\xi_{k})\right]-64\pi^{2}\sum_{i=1}^{l}\sum_{k=l+1}^{m}a(\xi_{k})G(\xi_{k},\xi_{i})\right.$
$\displaystyle\left.+4\pi\sum_{i=l+1}^{m}a(\xi_{i})\left[p-2\log
p+4\sum_{k=l+1,\,k\neq
i}^{m}\log|\xi_{i}-\xi_{k}|\right]+O\left(1\right)\right\\}\qquad\qquad\qquad\,\,\qquad\qquad\qquad$
(7.7)

$C^{0}$-uniformly in
$\mathcal{O}^{*}_{d,\,p}:=\left\\{\,\xi=(\xi_{1},\ldots,\xi_{m})\in\big{(}B_{d}(\xi_{*})\cap\Omega\big{)}^{l}\times\big{(}B_{d}(\xi_{*})\cap\partial\Omega\big{)}^{m-l}\,\left|\,\,\,\min\limits_{i,k=1,\ldots,m,\,i\neq
k}\,|\xi_{i}-\xi_{k}|>\frac{1}{p^{\kappa}},\right.\right.\\\
\left.\min\limits_{1\leq i\leq
l}\,\text{dist}(\xi_{i},\partial\Omega)>\frac{1}{p^{\kappa}}\right\\}$ for any
$d>0$ sufficiently small, fixed and independent of $p$. Here we claim that for
any $p>1$ large enough, the following maximization problem

$\displaystyle\max\limits_{(\xi_{1},\ldots,\xi_{m})\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi_{1},\ldots,\xi_{m})$

has a solution in the interior of $\mathcal{O}^{*}_{d,p}$. Once this claim is
proven, we can easily get the qualitative properties of solutions of (1.6) as
predicted in Theorem 1.3.

Let $\xi^{p}=(\xi^{p}_{1},\ldots,\xi^{p}_{m})$ be the maximizer of $F_{p}$
over $\overline{\mathcal{O}}^{*}_{d,p}$. We are led to prove that $\xi^{p}$
lies in the interior of $\mathcal{O}^{*}_{d,p}$. First, we obtain a lower
bound for $F_{p}$ over $\overline{\mathcal{O}}^{*}_{d,p}$. Let us consider a
smooth change of variables

$\displaystyle H_{\xi_{*}}^{p}(y)=e^{p/2}H_{\xi_{*}}(e^{-p/2}y),$

where $H_{\xi_{*}}:B_{d}(\xi_{*})\mapsto\mathcal{M}$ is a diffeomorphism and
$\mathcal{M}$ is an open neighborhood of the origin such that
$H_{\xi_{*}}(B_{d}(\xi_{*})\cap\Omega)=\mathcal{M}\cap\mathbb{R}_{+}^{2}$ and
$H_{\xi_{*}}(B_{d}(\xi_{*})\cap\partial\Omega)=\mathcal{M}\cap\partial\mathbb{R}_{+}^{2}$.
Let

$\displaystyle\xi^{0}_{i}=\xi_{*}-\frac{t_{i}}{\sqrt{p}}\nu(\xi_{*}),\quad
i=1,\ldots,l,\qquad\quad\textrm{but}\quad\qquad\xi^{0}_{i}=e^{-p/2}(H_{\xi_{*}}^{p})^{-1}\left(\frac{e^{p/2}}{\sqrt{p}}\hat{\xi}_{i}^{0}\right),\quad
i=l+1,\ldots,m,$

where $t_{i}>0$ and
$\hat{\xi}_{i}^{0}\in\mathcal{M}\cap\partial\mathbb{R}_{+}^{2}$ satisfy
$t_{i+1}-t_{i}=\rho$, $|\hat{\xi}_{i}^{0}-\hat{\xi}^{0}_{i+1}|=\rho$ for all
$\rho>0$ sufficiently small, fixed and independent of $p$. From the expansion
$(H_{\xi_{*}}^{p})^{-1}(z)=e^{p/2}\xi_{*}+z+O(e^{-p/2}|z|)$, we get

$\displaystyle\xi^{0}_{i}=\xi_{*}+\frac{1}{\sqrt{p}}\hat{\xi}_{i}^{0}+O\left(\frac{e^{-p/2}}{\sqrt{p}}|\hat{\xi}_{i}^{0}|\right),\quad
i=l+1,\ldots,m.$

Due to $\kappa>1$, we find
$\xi^{0}=(\xi^{0}_{1},\ldots,\xi^{0}_{m})\in\mathcal{O}^{*}_{d,p}$. Since
$\xi_{*}\in\partial\Omega$ is a strict local maximum point of $a(x)$ over
$\overline{\Omega}$ and satisfies $\partial_{\nu}a(\xi_{*})=\langle\nabla
a(\xi_{*}),\,\nu(\xi_{*})\rangle=0$, there exists a constant $C>0$ independent
of $p$ such that

$\displaystyle a(\xi_{*})-\frac{C}{p}\leq a(\xi^{0}_{i})<a(\xi_{*}),\qquad
i=1,\ldots,m.$

On the other hand, from definition (1.9), Lemma 2.2 and Corollary 2.3, we can
compute that for any $i=1,\ldots,l$ and $k=1,\ldots,m$ with $i\neq k$,

$\displaystyle H(\xi^{0}_{i},\xi^{0}_{i})=\frac{1}{4\pi}\log
p+O\left(1\right),\qquad\quad
G(\xi^{0}_{k},\xi^{0}_{i})=H(\xi^{0}_{k},\xi^{0}_{i})-\frac{1}{2\pi}\log|\xi^{0}_{k}-\xi^{0}_{i}|=\frac{1}{2\pi}\log
p+O\left(1\right).$

Additionally, for any $i,k=l+1,\ldots,m$ with $i\neq k$,

$\displaystyle\log|\xi^{0}_{i}-\xi^{0}_{k}|=-\frac{1}{2}\log
p+O\left(1\right).$

Hence by (7), we find

$\displaystyle\max\limits_{\xi\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi)\geq
F_{p}(\xi^{0})\geq\frac{e}{2p^{2}}\big{\\{}4\pi(m+l)pa(\xi_{*})-8\pi(m+l)^{2}a(\xi_{*})\log
p+O(1)\big{\\}}.$ (7.8)

Next, we suppose
$\xi^{p}=(\xi^{p}_{1},\ldots,\xi^{p}_{m})\in\partial\mathcal{O}^{*}_{d,p}$.
Then there exist four cases:
C1. There exists an $i_{0}\in\\{1,\ldots,l\\}$ such that
$\xi^{p}_{i_{0}}\in\partial B_{d}(\xi_{*})\cap\Omega$, in which case,
$a(\xi^{p}_{i_{0}})<a(\xi_{*})-d_{0}$ for some
$d_{0}>0$ independent of $p$;
C2. There exists an $i_{0}\in\\{l+1,\ldots,m\\}$ such that
$\xi^{p}_{i_{0}}\in\partial B_{d}(\xi_{*})\cap\partial\Omega$, in which case,
$a(\xi^{p}_{i_{0}})<a(\xi_{*})-d_{0}$ for
some $d_{0}>0$ independent of $p$;
C3. There exists an $i_{0}\in\\{1,\ldots,l\\}$ such that
$\text{dist}(\xi_{i_{0}}^{p},\partial\Omega)=p^{-\kappa}$;
C4. There exist indices $i_{0}$, $k_{0}$, $i_{0}\neq k_{0}$ such that
$|\xi_{i_{0}}^{p}-\xi_{k_{0}}^{p}|=p^{-\kappa}$.
Observe that for all $i=1,\ldots,l$ and $k=1,\ldots,m$ with $i\neq k$, by
(1.8), (2.3), (2.5) and the maximum principle,

$\displaystyle
G(\xi^{p}_{k},\xi^{p}_{i})>0,\qquad\,\,H(\xi^{p}_{k},\xi^{p}_{i})>0\qquad\,\,\textrm{and}\qquad\,\,H(\xi^{p}_{i},\xi^{p}_{i})=-\frac{1}{2\pi}\log\big{[}\text{dist}(\xi_{i}^{p},\partial\Omega)\big{]}+O\left(1\right).$
(7.9)

In the first and second cases, by (7) and (7.9) we get

$\displaystyle\max\limits_{\xi\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi)=F_{p}(\xi^{p})\leq\frac{e}{2p^{2}}\left\\{4\pi
p\big{[}(m+l)a(\xi_{*})-d_{0}\big{]}+O\left(\log p\right)\right\\},$ (7.10)

which contradicts to (7.8). This shows that $a(\xi_{i}^{p})\rightarrow
a(\xi_{*})$. Using the assumption of $a(x)$ over $\overline{\Omega}$, we
deduce $\xi_{i}^{p}\rightarrow\xi_{*}$ for all $i=1,\ldots,m$.

In the third case, by (7) and (7.9) we get

$\displaystyle\max\limits_{\xi\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi)=F_{p}(\xi^{p})\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})-64\pi^{2}a(\xi^{p}_{i_{0}})H(\xi^{p}_{i_{0}},\xi^{p}_{i_{0}})+O(1)\right\\}$
$\displaystyle\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})-32\pi\kappa a(\xi^{p}_{i_{0}})\log p+O(1)\right\\}.\quad\,\,\,\,$
(7.11)

In the last case, by (7) and (7.9) we get, if $i_{0}\in\\{1,\ldots,m\\}$ and
$k_{0}\in\\{1,\ldots,l\\}$,

$\displaystyle\max\limits_{\xi\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi)=F_{p}(\xi^{p})\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})+32\pi
a(\xi^{p}_{i_{0}})\log|\xi^{p}_{i_{0}}-\xi^{p}_{k_{0}}|+O(1)\right\\}$
$\displaystyle\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})-32\pi\kappa a(\xi^{p}_{i_{0}})\log
p+O(1)\right\\},\qquad\,\,\,\,\,$ (7.12)

while if $i_{0}\in\\{l+1,\ldots,m\\}$ and $k_{0}\in\\{l+1,\ldots,m\\}$,

$\displaystyle\max\limits_{\xi\in\overline{\mathcal{O}}^{*}_{d,p}}F_{p}(\xi)=F_{p}(\xi^{p})\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})+16\pi
a(\xi^{p}_{i_{0}})\log|\xi^{p}_{i_{0}}-\xi^{p}_{k_{0}}|+O(1)\right\\}$
$\displaystyle\leq\frac{e}{2p^{2}}\left\\{4\pi(m+l)(p-2\log
p)a(\xi_{*})-16\pi\kappa a(\xi^{p}_{i_{0}})\log
p+O(1)\right\\}.\qquad\,\,\,\,\,$ (7.13)

Comparing (7)-(7) with (7.8), we obtain

$\displaystyle 8\pi(m+l)a(\xi_{*})\log p+32\pi\kappa a(\xi^{p}_{i_{0}})\log
p+O(1)\leq 8\pi(m+l)^{2}a(\xi_{*})\log p+O(1),$ (7.14)

which is impossible by the choice of $\kappa$ in (2.10). $\square$

## Acknowledgments

This research is supported by the Fundamental Research Funds for the Central
Universities under Grant No. KYZ201649, and the National Natural Science
Foundation of China under Grant Nos. 11601232, 11671354 and 11775116.

##  8\. Appendix

Proof of lemma 2.4. Observe that for any $\beta\in(0,1)$, by (2.7), (2.12) and
(2.15),

$\displaystyle\left\\{\begin{aligned}
&-\Delta_{a}H_{i}+H_{i}=\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left\\{\left(-4+\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\left[\frac{(x-\xi_{i})\cdot\nabla\log
a(x)}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}-\frac{1}{2}\log\big{(}\delta_{i}^{2}+|x-\xi_{i}|^{2}\big{)}\right]-\log(8\delta_{i}^{2})\right.\\\\[5.69054pt]
&\left.\,\quad\qquad\qquad\qquad+\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{i}+\frac{1}{p}O_{\large
L^{\infty}\big{(}\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{\delta_{i}}{\delta_{i}+|x-\xi_{i}|}+\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\right.\\\\[5.69054pt]
&\left.\,\quad\qquad\qquad\qquad+\frac{1}{p}O_{\large
L^{\infty}\big{(}\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{|(x-\xi_{i})\cdot\nabla\log
a(x)|}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}+\log\frac{\delta_{i}^{2}+|x-\xi_{i}|^{2}}{\delta_{i}^{2}}\right)\right\\}\quad\quad\qquad\,\textrm{in}\,\,\,\,\
\,\,\Omega,\\\\[5.69054pt] &\frac{\partial
H_{i}}{\partial\nu}=-\frac{1}{\gamma\mu_{i}^{2/(p-1)}}\left\\{\left(-4+\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\frac{(x-\xi_{i})\cdot\nu(x)}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}+\frac{1}{p}O_{\large
L^{\infty}\big{(}\partial\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\right.\\\\[5.69054pt]
&\left.\,\quad\qquad+\frac{1}{p}O_{\large
L^{\infty}\big{(}\partial\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{|(x-\xi_{i})\cdot\nu(x)|}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\right\\}\,\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\textrm{on}\,\,\
\ \,\partial\Omega.\end{aligned}\right.$

From (1.8)-(1.9) we have that the regular part of Green’s function,
$H(x,\xi_{i})$, satisfies

$\left\\{\begin{aligned}
&-\Delta_{a}H(x,\xi_{i})+H(x,\xi_{i})=\frac{4}{c_{i}}\log|x-\xi_{i}|-\frac{4}{c_{i}}\frac{(x-\xi_{i})\cdot\nabla\log
a(x)}{|x-\xi_{i}|^{2}}\quad\,\,\textrm{in}\,\,\,\,\,\Omega,\\\\[2.84526pt]
&\frac{\partial
H(x,\xi_{i})}{\partial\nu}=\frac{4}{c_{i}}\frac{(x-\xi_{i})\cdot\nu(x)}{|x-\xi_{i}|^{2}}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\textrm{on}\,\,\,\partial\Omega.\end{aligned}\right.$
(8.1)

So, if we set

$\displaystyle
Z_{i}(x)=\gamma\mu_{i}^{2/(p-1)}H_{i}(x)-\left(4-\frac{C_{1}}{p}-\frac{C_{2}}{p^{2}}\right)\frac{c_{i}}{4}H(x,\xi_{i})+\log(8\delta_{i}^{2})-\left(\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\log\delta_{i},$

then $Z_{i}(x)$ satisfies

$\displaystyle\left\\{\begin{aligned}
&-\Delta_{a}Z_{i}+Z_{i}=\left(-4+\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\left[\,\frac{1}{2}\log\left(\frac{|x-\xi_{i}|^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)-\frac{(x-\xi_{i})\cdot\nabla\log
a(x)}{|x-\xi_{i}|^{2}}\cdot\frac{\delta_{i}^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right]\\\\[5.69054pt]
&\quad\qquad\qquad\qquad+\frac{1}{p}O_{\large L^{\infty}\big{(}\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{\delta_{i}}{\delta_{i}+|x-\xi_{i}|}+\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\\\\[5.69054pt]
&\quad\qquad\qquad\qquad+\frac{1}{p}O_{\large L^{\infty}\big{(}\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{|(x-\xi_{i})\cdot\nabla\log
a(x)|}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}+\log\frac{\delta_{i}^{2}+|x-\xi_{i}|^{2}}{\delta_{i}^{2}}\right)\quad\qquad\,\textrm{in}\,\,\
\,\,\ \,\,\Omega,\\\\[5.69054pt] &\frac{\partial
Z_{i}}{\partial\nu}=\left(-4+\frac{C_{1}}{p}+\frac{C_{2}}{p^{2}}\right)\frac{(x-\xi_{i})\cdot\nu(x)}{|x-\xi_{i}|^{2}}\cdot\frac{\delta_{i}^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}+\frac{1}{p}O_{\large
L^{\infty}\big{(}\partial\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\\\\[5.69054pt]
&\qquad\quad+\frac{1}{p}O_{\large L^{\infty}\big{(}\partial\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})\big{)}}\left(\frac{|(x-\xi_{i})\cdot\nu(x)|}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\,\qquad\qquad\quad\qquad\qquad\quad\qquad\quad\qquad\qquad\,\,\,\textrm{on}\,\
\,\,\,\partial\Omega.\end{aligned}\right.$

Using polar coordinates with center $\xi_{i}$, i.e. $r=|x-\xi_{i}|$, and
changing variables $s=r/\delta_{i}$, we estimate that for any $q>1$,

$\displaystyle\int_{\Omega}\left|\log\left(\frac{|x-\xi_{i}|^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right)\right|^{q}dx$
$\displaystyle\leq
2\pi\int_{0}^{\text{diam}(\Omega)}\left|\log\left(\frac{r^{2}}{\delta_{i}^{2}+r^{2}}\right)\right|^{q}rdr$
$\displaystyle\leq
2\pi\delta_{i}^{2}\int_{0}^{+\infty}\left|\log\left(1+\frac{1}{s^{2}}\right)\right|^{q}sds$
$\displaystyle\leq C\delta_{i}^{2},$

and

$\displaystyle\int_{\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})}\left|\frac{\delta_{i}}{\delta_{i}+|x-\xi_{i}|}+\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|^{q}dx$
$\displaystyle\leq
C\int_{\delta_{i}^{\beta/2}}^{\text{diam}(\Omega)}\left[\frac{\delta_{i}^{q}}{(\delta_{i}+r)^{q}}+\frac{\delta_{i}^{q}}{(\delta_{i}^{2}+r^{2})^{q}}\right]rdr$
$\displaystyle\leq
C\int_{\delta_{i}^{\beta/2-1}}^{+\infty}\left[\frac{\delta_{i}^{2}}{(1+s)^{q}}+\frac{\delta_{i}^{2-q}}{(1+s^{2})^{q}}\right]sds$
$\displaystyle\leq C\delta_{i}^{\beta+q(1-\beta)},$

and for any $1<q<2$,

$\displaystyle\int_{\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})}\left|\frac{|(x-\xi_{i})\cdot\nabla\log
a(x)|}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}+\log\frac{\delta_{i}^{2}+|x-\xi_{i}|^{2}}{\delta_{i}^{2}}\right|^{q}dx\leq$
$\displaystyle
C\int_{0}^{\delta_{i}^{\beta/2}}\left[\left(\frac{r}{\delta_{i}^{2}+r^{2}}\right)^{q}+\left|\log\left(\frac{\delta_{i}^{2}+r^{2}}{\delta_{i}^{2}}\right)\right|^{q}\right]rdr$
$\displaystyle=$ $\displaystyle
C\delta_{i}^{2}\int_{0}^{\delta_{i}^{\beta/2-1}}\left[\delta_{i}^{-q}\left(\frac{s}{1+s^{2}}\right)^{q}+\log^{q}\big{(}1+s^{2}\big{)}\right]sds$
$\displaystyle\leq$ $\displaystyle C\delta_{i}^{\beta(1-q/2)},$

and

$\displaystyle\int_{\Omega}\left|\frac{(x-\xi_{i})\cdot\nabla\log
a(x)}{|x-\xi_{i}|^{2}}\cdot\frac{\delta_{i}^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|^{q}dx$
$\displaystyle\leq
C\int_{0}^{\text{diam}(\Omega)}\left|\frac{\delta_{i}^{2}}{r(\delta_{i}^{2}+r^{2})}\right|^{q}rdr$
$\displaystyle\leq
C\delta_{i}^{2-q}\int_{0}^{+\infty}\frac{s^{1-q}}{(1+s^{2})^{q}}ds$
$\displaystyle\leq C\delta_{i}^{2-q}.$

Thus for any $\xi_{i}\in\overline{\Omega}$ and any $1<q<2$,

$\displaystyle\big{\|}-\Delta_{a}Z_{i}+Z_{i}\big{\|}_{L^{q}(\Omega)}\leq
C\delta_{i}^{\beta(1/q-1/2)}.$ (8.2)

As for the boundary terms, if $\xi_{i}\in\Omega$, by (2.9) we get, for any
$x\in\partial\Omega$,

$\displaystyle\left|\frac{(x-\xi_{i})\cdot\nu(x)}{|x-\xi_{i}|^{2}}\cdot\frac{\delta_{i}^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|\leq\frac{\delta_{i}^{2}}{|x-\xi_{i}|^{3}}\leq\delta_{i}^{2}p^{3\kappa},$
$\displaystyle\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\leq\frac{\delta_{i}}{|x-\xi_{i}|^{2}}\leq\delta_{i}p^{2\kappa},$

and further,

$\displaystyle\left\|\frac{\partial
Z_{i}}{\partial\nu}\right\|_{L^{\infty}(\partial\Omega)}\leq
C\delta_{i}p^{2\kappa-1}.$ (8.3)

While if $\xi_{i}\in\partial\Omega$, using the fact that
$|(x-\xi_{i})\cdot\nu(x)|\leq C|x-\xi_{i}|^{2}$ for any $x\in\partial\Omega$
(see [5]), we estimate that for any $q>1$,

$\displaystyle\int_{\partial\Omega}\left|\frac{(x-\xi_{i})\cdot\nu(x)}{|x-\xi_{i}|^{2}}\cdot\frac{\delta_{i}^{2}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|^{q}dx$
$\displaystyle\leq
C\int_{0}^{+\infty}\frac{\delta_{i}^{2q}}{(\delta_{i}^{2}+r^{2})^{q}}dr=C\delta_{i}\int_{0}^{+\infty}\frac{1}{(1+s^{2})^{q}}ds\leq
C\delta_{i},$ $\displaystyle\int_{\partial\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})}\left|\frac{(x-\xi_{i})\cdot\nu(x)}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|^{q}dx$
$\displaystyle\leq C\int_{\partial\Omega\bigcap
B_{\delta_{i}^{\beta/2}}(\xi_{i})}\frac{|x-\xi_{i}|^{2q}}{(\delta_{i}^{2}+|x-\xi_{i}|^{2})^{q}}dx\leq
C\left|\partial\Omega\cap B_{\delta_{i}^{\beta/2}}(\xi_{i})\right|\leq
C\delta_{i}^{\beta/2},$

and

$\displaystyle\int_{\partial\Omega\setminus
B_{\delta_{i}^{\beta/2}}(\xi_{i})}\left|\frac{\delta_{i}}{\delta_{i}^{2}+|x-\xi_{i}|^{2}}\right|^{q}dx$
$\displaystyle\leq
C\int_{\delta_{i}^{\beta/2}}^{|\partial\Omega|}\frac{\delta_{i}^{q}}{(\delta_{i}^{2}+r^{2})^{q}}dr\leq
C\int_{\delta_{i}^{\beta/2}}^{|\partial\Omega|}\frac{\delta_{i}^{q}}{\,\,r^{2q}\,\,}dr\leq
C\big{(}\delta_{i}^{q}+\delta_{i}^{q(1-\beta)+\beta/2}\big{)}.$

Thus for any $\xi_{i}\in\partial\Omega$ and any $q>1$,

$\displaystyle\left\|\frac{\partial
Z_{i}}{\partial\nu}\right\|_{L^{q}(\partial\Omega)}\leq
C\delta_{i}^{\beta/2q}.$ (8.4)

Hence by elliptic regularity theory, we obtain that for any $1<q<2$ and any
$0<\theta<1/q$,

$\displaystyle\left\|Z_{i}\right\|_{W^{1+\theta,q}(\Omega)}\leq
C\left(\big{\|}-\Delta_{a}Z_{i}+Z_{i}\big{\|}_{L^{q}(\Omega)}+\left\|\frac{\partial
Z_{i}}{\partial\nu}\right\|_{L^{q}(\partial\Omega)}\right)\leq
C\delta_{i}^{\beta(1/q-1/2)}.$

Then by Morrey embedding,

$\displaystyle\left\|Z_{i}\right\|_{C^{\tau}(\overline{\Omega})}\leq
C\delta_{i}^{\beta(1/q-1/2)},$

where $0<\tau<1/2+1/q$, which implies that expansion (2.22) holds with
$\alpha=2\beta(1/q-1/2)$. $\square$

Proof of lemma 2.5. Notice that if we make the change of variables $s=1/p$,
then system (2.25) can be rewritten as the form of the following vector
function

$\displaystyle S_{i}(s,\xi,\mu):=\log\mu_{i}+\frac{\,\log
8+\frac{1}{4}C_{1}+\frac{1}{4}C_{2}s\,}{\,4\big{(}1-\frac{1}{4}C_{1}s-\frac{1}{4}C_{2}s^{2}\big{)}\,}-\frac{1}{4}c_{i}H(\xi_{i},\xi_{i})-\frac{1}{4}\sum_{k=1,\,k\neq
i}^{m}\left(\frac{\,\mu_{i}\,}{\mu_{k}}\right)^{2s/(1-s)}c_{k}G(\xi_{i},\xi_{k})=0.$

Obviously, from the explicit expression (2.18) of the constant $C_{1}$, we
have that for $s=0$,

$\displaystyle\mu_{i}(0,\xi)=e\large^{-\frac{3}{4}+\frac{1}{4}c_{i}H(\xi_{i},\xi_{i})+\frac{1}{4}\sum_{k=1,\,k\neq
i}^{m}c_{k}G(\xi_{i},\xi_{k})}.$ (8.5)

Using the Taylor expansion of exponential function, we can conclude that for
any $s>0$ small enough,

$\displaystyle\left(\frac{\,C^{2}\,}{s^{\kappa}}\right)^{2s/(1-s)}=e\large^{\frac{2s}{1-s}\big{(}2\log
C+\kappa\log\frac{1}{s}\big{)}}=1+O\left(s\log\frac{1}{s}\right),$

and then, by (2.26),

$\displaystyle\left(\frac{\,\mu_{i}\,}{\mu_{k}}\right)^{2s/(1-s)}=1+O\left(s\log\frac{1}{s}\right),\qquad\forall\,\,i\neq
k.$ (8.6)

Moreover, by (1.9), (2.3) and the fact that
$a(\xi_{i})G(\xi_{i},\xi_{k})=a(\xi_{k})G(\xi_{k},\xi_{i})$ for all
$i,k=1,\ldots,m$ with $i\neq k$, we can easily prove that for any points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{1/s}$,

$\displaystyle
G(\xi_{i},\xi_{k})=O\left(\log\frac{1}{s}\right),\,\qquad\,\forall\,\,i\neq
k.$ (8.7)

Consequently, a simple computation shows that

$\displaystyle\frac{\partial
S_{i}(s,\xi,\mu)}{\partial\mu_{i}}=\frac{1}{\mu_{i}}-\frac{s}{\,2\mu_{i}(1-s)\,}\sum_{k=1,\,k\neq
i}^{m}\left(\frac{\,\mu_{i}\,}{\mu_{k}}\right)^{2s/(1-s)}c_{k}G(\xi_{i},\xi_{k})=\frac{1}{\mu_{i}}\left[1+O\left(s\log\frac{1}{s}\right)\right],$

and for any $k\neq i$,

$\displaystyle\frac{\partial
S_{i}(s,\xi,\mu)}{\partial\mu_{k}}=\frac{s}{\,2\mu_{k}(1-s)\,}\left(\frac{\,\mu_{i}\,}{\mu_{k}}\right)^{2s/(1-s)}c_{k}G(\xi_{i},\xi_{k})=\frac{1}{\mu_{k}}O\left(s\log\frac{1}{s}\right).$

Then

$\displaystyle\det\big{(}\nabla_{\mu}S(s,\xi,\mu)\big{)}=\frac{1}{\,\mu_{1}\ldots\mu_{m}\,}\left[1+O\left(s\log\frac{1}{s}\right)\right]\neq
0.$

Hence $\nabla_{\mu}S(s,\xi,\mu)$ is invertible in the range of points and
variables that we are considering. From the Implicit Function Theorem we find
that $S(s,\xi,\mu)=0$ is solvable in some neighborhood of
$\big{(}0,\xi,\mu(0,\xi)\big{)}$, and thus for any points
$\xi=(\xi_{1},\ldots,\xi_{m})\in\mathcal{O}_{p}$ and any $p>1$ large enough,
system (2.25) has a unique solution $\mu=(\mu_{1},\ldots,\mu_{m})$ satisfying
(2.26). This, together with (8.5)-(8.7), implies

$\displaystyle\mu_{i}=e\large^{-\frac{3}{4}+\frac{1}{4}c_{i}H(\xi_{i},\xi_{i})+\frac{1}{4}\sum\large_{k=1,\,k\neq
i}^{m}c_{k}G(\xi_{i},\xi_{k})}\left[\,1+O\left(\frac{\log^{2}p}{p}\right)\right],\quad\,\,\forall\,\,i=1,\ldots,m.$

Moreover, by (2.3), (2.5), (2.9), (2.25), (8.6) and (8.7) we can conclude that
estimate (2.27) holds. $\square$

## References

  * [1] Adimurthi, G. Mancini, The Neumann problem for elliptic equations with critical nonlinearity, “A tribute in honour of G. Prodi”, Scuola Norm. Sup. Pisa (1991) 9–25.
  * [2] Adimurthi, G. Mancini, Geometry and topology of the boundary in the critical Neumann problem, J. Reine Angew. Math. 456 (1994) 1–18.
  * [3] Adimurthi, G. Mancini, S.L. Yadava, The role of the mean curvature in semilinear Neumann problem involving critical exponent, Comm. Partial Differential Equations 20 (1995) 591–631.
  * [4] Adimurthi, F. Pacella, S.L. Yadava, Interaction between the geometry of the boundary and positive solutions of a semilinear Neumann problem with critical nonlinearity, J. Funct. Anal. 113 (1993) 318–350.
  * [5] O. Agudelo, A. Pistoia, Boundary concentration phenomena for the higher-dimensional Keller-Segel system, Calc. Var. Partial Differential Equations (2016) 55:132.
  * [6] S. Baraket and F. Parcard, Construction of singular limits for a semilinear elliptic equation in dimension 2, Calc. Var. Partial Differential Equations 6 (1998) 1–38.
  * [7] D. Chae, O. Imanuvilov, The existence of non-topological multivortex solutions in the relativistic self-dual Chern-Simons theory, Comm. Math. Phys. 215 (2000) 119–142.
  * [8] C.C. Chen, C.S. Lin, Sharp estimates for solutions of multi-bubbles in compact Riemann surfaces, Comm. Pure Appl. Math. 55 (2002) 728–771.
  * [9] D. Cao, E.S. Noussair, S. Yan, Existence and nonexistence of interior-peaked solution for a nonlinear Neumann problem, Pacific J. Math. 200 (2001) 19–42.
  * [10] D. Cao, E.S. Noussair, The effect of geometry of the domain boundary in an elliptic Neumann problem, Adv. Differential Equations 6 (2001) 931–958.
  * [11] G. Cerami, J. Wei, Multiplicity of multiple interior peak solutions for some singularly perturbed Neumann problems, Internat. Math. Res. Notices 12 (1998) 601–626.
  * [12] M. del Pino, P.L. Felmer, J. Wei, On the role of mean curvature in some singularly perturbed Neumann problems, SIAM J. Math. Anal. 31 (1999) 63–79.
  * [13] M. del Pino, P.L. Felmer, J. Wei, On the role of distance function in some singular perturbation problems, Comm. Partial Differential Equations 25 (2000) 155–177.
  * [14] M. del Pino, M. Kowalczyk, M. Musso, Singular limits in Liouville-type equations, Calc. Var. Partial Differential Equations 24 (2005) 47–81.
  * [15] M. del Pino, F. Mahmoudi, M. Musso, Bubbling on boundary submanifolds for the Lin-Ni-Takagi problem at higher critical exponents, J. Eur. Math. Soc. 16 (2014) 1687–1748.
  * [16] M. del Pino, M. Musso, A. Pistoia, Super-critical boundary bubbling in a semilinear Neumann problem, Ann. Inst. H. Poincaré Analyse Non Linéaire 22 (2005) 45–82.
  * [17] M. del Pino, J. Wei, Collapsing steady states of the Keller-Segel system, Nonlinearity 19 (2006) 661–684.
  * [18] S. Deng, F. Mahmoudi, M. Musso, Bubbling on boundary submanifolds for a semiliear Neumann problem near higher critical exponents, Discrete Contin. Dynam. Systems A 36 (2016) 3035–3076.
  * [19] E.N. Dancer, S. Yan, Multipeak solutions for a singularly perturbed Neumann problem, Pacific J. Math. 189 (1999) 241–262.
  * [20] P. Esposito, Interior estimates for some semilinear elliptic problem with critical nonlinearity, Ann. Inst. H. Poincaré Analyse Non Linéaire 24 (2007) 629–644.
  * [21] P. Esposito, M. Musso, A. Pistoia, Concentrating solutions for a planar elliptic problem involving nonlinearities with large exponent, J. Differential Equations 227 (2006) 29–68.
  * [22] A. Gierer, H. Meinhardt, A theory of biological pattern formation, Kybernetik(Berlin) 12 (1972) 30–39.
  * [23] M. Grossi, A class of solutions for the Neumann problem $-\Delta u+\lambda u=u^{\frac{N+2}{N-2}}$, Duke Math. J. 79 (1995) 309-335.
  * [24] M. Grossi, A. Pistoia, J. Wei, Existence of multipeak solutions for a semilinear Neumann problem via nonsmooth critical point theory, Calc. Var. Partial Differential Equations 11 (2000) 143–175.
  * [25] N. Ghoussoub, C. Gui, M. Zhu, On a singularly perturbed Neumann problem with critical exponent, Comm. Partial Differential Equation 26 (2001) 1929–1946.
  * [26] C. Gui, N. Ghoussoub, Multi-peak solutions for a semilinear Neumann problem involving the critical Sobolev exponent, Math. Z. 229 (1998) 443–474.
  * [27] C. Gui, C.S. Lin, Estimates for boundary-bubbling solutions to an elliptic Neumann problem, J. Reine Angew. Math. 546 (2002) 201–235.
  * [28] C. Gui, J. Wei, Multiple interior spike solutions for some singular perturbed Neumann problems, J. Differential Equations, 158 (1999) 1–27.
  * [29] C. Gui, J. Wei, On multiple mixed interior and boundary peak solutions for some singularly perturbed Neumann problems, Canad. J. Math. 52 (2000) 522–538.
  * [30] C. Gui, J. Wei, M. Winter, Multiple boundary peak solutions for some singularly perturbed Neumann problems, Ann. Inst. H. Poincaré Anal. Non Linéaire, 17 (2000) 249–289.
  * [31] M. Kowalczyk, Multiple spike layers in the shadow Gierer-Meinhardt system: existence of equilibria and the quasi-invariant manifold, Duke Math. J. 98 (1999) 59–111.
  * [32] E.F. Keller, L.A. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theoretical Biology 26 (1970) 399–415.
  * [33] C.S. Lin, W.M. Ni, I. Takagi, Large amplitude stationary solutions to a chemotaxis system, J. Differential Equations 72 (1988) 1–27.
  * [34] C.S. Lin, L. Wang, J. Wei, Bubble accumulations in an elliptic Neumann problem with critical Sobolev exponent, Calc. Var. Partial Differential Equations 30 (2007) 153–182.
  * [35] Y.-Y. Li, On a singularly perturbed equation with Neumann boundary condition, Comm. Partial Differential Equations 23 (1998) 487–545.
  * [36] F. Mahmoudi, A. Malchiodi, Concentration on minimal submanifolds for a singularly perturbed Neumann problem, Adv. Math. 209 (2007) 460–525.
  * [37] A. Malchiodi, Concentration at curves for a singularly perturbed Neumann problem in three-dimensional domains, Geom. Funct. Anal. 15 (2005) 1162–1222.
  * [38] A. Malchiodi, M. Montenegro, Boundary concentration phenomena for a singularly perturbed elliptic problem, Comm. Pure Appl. Math. 55 (2002) 1507–1568.
  * [39] A. Malchiodi, M. Montenegro, Multidimensional boundary-layers for a singularly perturbed Neumann problem, Duke Math. J. 124 (2004) 105–143.
  * [40] B. Manna, A. Pistoia, Boundary-layers for a Neumann problem at higher critical exponents, Boll. Unione Mat. Ital. 10 (2017) 355–368.
  * [41] M. Musso, J. Wei, Stationary solutions to a Keller-Segel chemotaxis system, Asymptotic Anal. 49 (2006) 217–247.
  * [42] W.M. Ni, Diffusion, cross-diffusion, and their spike-layer steady states, Notices Amer. Math. Soc. 45 (1998) 9–18.
  * [43] W.N. Ni, X.B. Pan, I. Takagi, Singular behavior of least-energy solutions of a semi-linear Neumann problem involving critical Sobolev exponents, Duke Math. J. 67 (1992) 1–20.
  * [44] W.M. Ni, I. Takagi, On the shape of least-energy solutions to a semilinear Neumann problem, Comm. Pure Appl. Math. 44 (1991) 819–851.
  * [45] W.M. Ni, I. Takagi, Locating the peaks of least-energy solutions to a semilinear Neumann problem, Duke Math. J. 70 (1993) 247–281.
  * [46] M. H. Protter, H. F. Weinberger, Maximum Principles in Differential Equation, Corrected Reprint of the 1967 Original, Springer-Verlag, New York, 1984.
  * [47] O. Rey, Boundary effect for an elliptic Neumann problem with critical nonlinearity, Comm. Partial Differential Equations 22 (1997) 1055–1139.
  * [48] O. Rey, An elliptic Neumann problem with critical nonlinearity in three dimensional domains, Comm. Contemp. Math. 1 (1999) 405–449.
  * [49] O. Rey, The question of interior blow-up points for an elliptic Neumann problem: the critical case, J. Math. Pures Appl. 81 (2002) 655–696.
  * [50] O. Rey, J. Wei, Blow-up solutions for an elliptic Neumann problem with sub-or-supcritical nonlinearity, II: $N=4$, Ann. Inst. H. Poincaré Analyse Non Linéaire 22 (2005) 459–484.
  * [51] O. Rey, J. Wei, Blow-up solutions for an elliptic Neumann problem with sub-or-supcritical nonlinearity, I: $N=3$, J. Funct. Anal. 212 (2004) 472–499.
  * [52] L. Wang, Arbitrarily many solutions for an elliptic Neumann problem with sub- or supercritical nonlinearity, Commun. Pure Appl. Anal. 9 (2010) 761–778.
  * [53] X.J. Wang, Neumann problem of semilinear elliptic equations involving critical Sobolev exponents, J. Differential Equations 93 (1991) 283–310.
  * [54] J. Wei, On the boundary spike layer solutions of singularly perturbed semilinear Neumann problem, J. Differential Equations 134 (1997) 104–133.
  * [55] J. Wei, On the interior spike layer solutions to a singularly perturbed Neumann problem, Tohoku Math. J. 50 (1998) 159–178.
  * [56] J. Wei, M. Winter, Multiple boundary spike solutions for a wide class of singular perturbation problems, J. Lond. Math. Soc. 59 (1999) 585–606.
  * [57] J. Wei, S. Yan, Arbitrary many boundary peak solutions for an elliptic Neumann problem with critical growth, J. Math. Pures Appl. 88 (2007) 350–378.