# Magnetic response and electronic states of well defined Gr/Fe/Ir(111)
heterostructure

Claudia Cardoso  Centro S3, CNR-Istituto Nanoscienze, I-41125 Modena, Italy
Giulia Avvisati  Dipartimento di Fisica, Universit$\grave{a}$ di Roma “La
Sapienza”, I-00185 Roma, Italy  Pierluigi Gargiani  ALBA Synchrotron Light
Source, Carrer de la Llum, 2-26 08290 Barcelona (Spain)  Marco Sbroscia
Dipartimento di Fisica, Universit$\grave{a}$ di Roma “La Sapienza”, I-00185
Roma, Italy  Madan S. Jagadeesh  Dipartimento di Fisica, Universit$\grave{a}$
di Roma “La Sapienza”, I-00185 Roma, Italy  Carlo Mariani  Dipartimento di
Fisica, Universit$\grave{a}$ di Roma “La Sapienza”, I-00185 Roma, Italy  Dario
A. Leon  Universit$\grave{a}$ di Modena e Reggio Emilia I-41125 Modena, Italy
Centro S3, CNR-Istituto Nanoscienze, I-41125 Modena, Italy  Daniele Varsano
Centro S3, CNR-Istituto Nanoscienze, I-41125 Modena, Italy  Andrea Ferretti
Centro S3, CNR-Istituto Nanoscienze, I-41125 Modena, Italy  Maria Grazia Betti
Dipartimento di Fisica, Universit$\grave{a}$ di Roma “La Sapienza”, I-00185
Roma, Italy

###### Abstract

We investigate a well-defined heterostructure constituted by magnetic Fe
layers sandwiched between graphene (Gr) and Ir(111). The challenging task to
avoid Fe-C solubility and Fe-Ir intermixing has been achieved with atomic
controlled Fe-intercalation at moderate temperature below 500 K. Upon
intercalation of a single ordered Fe layer in registry with the Ir substrate,
an intermixing of the Gr bands and Fe $d$-states breaks the symmetry of the
Dirac cone, with a downshift in energy of the apex by about 3 eV, and well-
localized Fe intermixed states induced in the energy region just below the
Fermi level. First principles electronic structure calculations show a large
spin splitting of the Fe states, resulting in a majority spin channel almost
fully occupied and strongly hybridized with Gr $\pi$-states. X-Ray Magnetic
Circular Dichroism on the Gr/Fe/Ir(111) heterostructure reveals an ordered
spin configuration with a ferromagnetic response of Fe layer(s), with enhanced
spin- and orbital-configurations with respect to the bcc-Fe bulk values. The
magnetization switches from a perpendicular easy magnetization axis when the
Fe single layer is lattice matched with the Ir(111) surface, to a parallel one
when the Fe thin film is almost commensurate with graphene.

intercalation; photoemission; graphene; moiré structure; Ir(111)

##  I Introduction

Graphene can be a very promising spin channel material owing to its long spin-
diffusion lengths of several micrometers Tombros2007 , gate-tunable carrier
concentration, and high electronic mobility Han2014 . Graphene coupled with
ferromagnetic systems can open new perspectives when efficient injection of
spin-polarized electrons can be achieved, as observed for a graphene membrane
on Co(0001) Usachov_2015 . It is well known that nearly flat epitaxial
graphene of high structural quality can be formed on several magnetic 3$d$
metal substrates, like the lattice-matched Ni(111) PhysRevB.84.205431 ;
Pacile2013 ; Varykhalov_PRX_2012 ; massimi2014  and Co(0001) Decker_PRB_2013 ;
Usachov_2015 ; Pacile_PRB_2014 ; Vita_PRB_2014  surfaces. It has been shown
that a tiny magnetic moment arises on the carbon atoms, induced by the strong
hybridization of graphene $\pi$-orbitals with Ni or Co 3$d$-states
Usachov_2015 ; Weser2010 ; Decker_JPhys_2014 ; Decker_PRB_2013 . Furthermore,
graphene grown on metals protects highly reactive magnetic surfaces and
stabilizes them against oxidation Coraux_JPCL_2012 ; Cattelan_NanoScale_2015 ;
lodesani2019graphene . Whereas recently a large research effort has been
dedicated to investigate Gr-Ni and Gr-Co heterostructures, only a few
experimental results for graphene grown on Fe surfaces are available
Vinogradov_2012 ; Vinogradov_2017 ; Varykhalov2015 , though iron is the most
widespread transition metal, and the technology of passivated Fe films with a
graphene membrane can be appealing for several reasons. Among them, the
considerable price reduction in comparison with other transition metal
substrates and, most importantly, its strong magnetic response.

The main hurdle for the formation of graphene on top of Fe surfaces is related
to the rich Fe-C phase diagram Okamoto1992 . In fact, the high carbon
diffusion coefficient and solubility in iron are detrimental for chemical
vapor deposition processes, where a high annealing temperature is necessary to
ensure high quality graphene on top of the metal surfaces. Thus, the epitaxial
growth of a graphene membrane on a Fe film/single crystal is made difficult
because of the formation of iron carbide, which is thermodynamically favored
Okamoto1992 . The epitaxial growth of graphene on Fe is also prevented because
none of the bcc-Fe facets is lattice-matched with graphene, at variance with
close-packed surfaces of other 3$d$ metals, like Co and Ni. Recently, a
detailed structural study has demonstrated that the Fe(110) face, with its
distorted hexagonal symmetry, can be a good candidate, because of a partial
match with the graphene lattice Vinogradov_2012 ; Vinogradov_2017 . In that
study, the substrate temperature was kept below 675 K to guarantee the
formation of a graphene membrane on iron reducing the presence of iron
carbides or segregation of diluted carbon from iron.

A more challenging strategy has been proposed exploiting Fe intercalation
beneath graphene Vlaic_JPCL_2018  as a viable route to overcome the hurdles to
realize a direct growth on Fe surfaces. Intercalation of metals below Gr has
proven to efficiently lead to the formation of atomically smooth metallic
layer(s) Decker_PRB_2013 ; Pacile_PRB_2014 ; Vita_PRB_2014 ; Decker_JPhys_2014
; Avvisati_ASS_2020 , in a layer-by-layer growth mode Gargiani_NatCom_2017 ,
where the graphene cover acts as a protective membrane of the confined
epitaxial metallic layer(s) Coraux_JPCL_2012 ; Cattelan_NanoScale_2015 .
Recently, a scanning tunnelling microscopy (STM) experiment demonstrates that
successful Fe intercalation under graphene grown on Ir(111) can be obtained
with the substrate temperature in the range between 450 K and 550 K
Bazarnik_SurfSci_2015 , giving rise to a smooth Fe layer pseudomorphic with
Ir(111) and a highly corrugated graphene membrane on top.

In this work, following this successful strategy for Fe growth, we have
realized well-defined smooth Fe layers intercalated between Gr and Ir(111),
preventing any alloying and Fe-C intermixing. By means of a combined
experimental and theoretical approach, we gain a detailed insight into the
physical properties of the iron-intercalated Gr/Ir(111) heterostructure. The
sandwiched Fe layers present a redistribution of the minority and majority
electronic density of states triggered by the spatial confinement and by the
peculiar strained structural configurations, as predicted by first principles
spin resolved electronic structure calculations and experimentally confirmed
by angular resolved photoemission and X-Ray magnetic dichroism.

##  II Methods

Angular resolved photoelectron spectroscopy (ARPES) and low-energy electron-
diffraction (LEED) experiments were carried out at the Nanostructure at
Surfaces laboratory, Department of Physics, Sapienza University LOTUS . The
ARPES apparatus is equipped with a Scienta SE200 multi-channel-plate electron
analyzer and a monochromatic Gammadata VUV 5000 microwave He source, with main
lines at 20.218 eV (He Iα) and 40.814 eV (He IIα) photon energies. The ARPES
apparatus was set for an experimental energy and angular resolution of 16 meV
and 0.18∘, respectively.

The X-ray Magnetic Circular Dichroism (XMCD) measurements were performed at
the BOREAS beamline of the ALBA synchrotron radiation facility boreas . Data
were taken in the total electron yield (TEY) mode, by measuring the drain
current with respect to a clean gold grid, used for photon flux normalization.
We used two different geometries, grazing incidence (GI) with 70∘ incidence
angle and normal incidence (NI), so to probe the magnetic signal along the
easy and hard magnetic axes; further details are available in
Avvisati_NanoLett_2018 .

In all laboratories, the Ir(111) surface was cleaned by cycles of ion
sputtering (Ar+, 1.5-2.0 keV) and annealing (1 minute at temperatures higher
than 1300 K). The quality of the obtained surface is confirmed by the sharp
LEED pattern. The Gr layer was prepared by thermal decomposition of ethylene,
by exposing the clean Ir surface to 5$\times$10-8-2$\times$10-7 mbar of C2H4
and annealing at around 1300-1320 K. Metallic Fe was deposited, at 0.3 Å/min
on the Gr/Ir(111) surface kept at about 500 K, in order to ensure a high
quality of Fe layer(s) intercalated beneath Gr and to avoid any intermixing,
that occurs at higher annealing temperature (600-800K) Brede2016 . The C $1s$,
Ir $4f$ and Fe $2p_{3/2}$ core-levels as measured at the Superesca beamline of
the Elettra synchrotron radiation facility, are presented in the Supplemental
Material supp-info .

The completion of a single layer of Fe can be identified by following the
evolution of the $\pi$ states of Gr in the ARPES valence band, as discussed in
detail in the next sections. At the BOREAS beamline, the Fe density in the
intercalated film was evaluated via its jump-edge ratio, previously calibrated
with Auger electron spectroscopy, as reported in the Supporting Information of
Ref. Gargiani_NatCom_2017 .

Density functional theory (DFT) simulations were carried out using the Quantum
ESPRESSO package Giannozzi2009 ; Giannozzi2017  where wavefunctions are
expanded in plane-waves and pseudopotentials are used to account for the
electron-ion interaction. We used the local density approximation (LDA) for
the exchange-correlation potential, according to the Perdew-Zunger
parametrization Perdew_PRB_1981 . Similarly to our previous works
Avvisati_NanoLett_2018 ; CalloniJCP2020 , we simulated the Gr/1ML-Fe/Ir(111)
interface considering the complete moiré induced periodicity by using a
9$\times$9 supercell of Ir(111), corresponding to a 10$\times$10 supercell of
pristine Gr. The lattice parameters were obtained by relaxing Ir bulk at the
LDA level using ultrasoft pseudopotentials (USPP), resulting in a Ir–Ir bond
distance of 2.7048 Å (corresponding to a hexagonal cell of 46.001 Bohr radius
for the moiré structure). In all the calculation we included four metallic
layers (3 Ir plus one Fe layer or 4 Ir layers). In order to make the two sides
of the slab inequivalent we added a layer of H atoms in one of the two sides.
Atomic positions were then fully relaxed (except for the two bottom Ir layers
and the H saturation layer) until ionic forces were smaller than 0.001
Ry/Bohr. For all the self-consistent calculations we used a 2$\times$2 grid of
$\mathbf{k}$-points, ultra soft pseudopotentials to model the electron-ion
interaction and a kinetic energy cutoff of 30 and 300 Ry to represent Kohn-
Sham wavefunctions and density, respectively.

Figure 1:  Low energy electron diffraction (LEED) patterns (electron beam
energy 140 eV) for (a) Gr/Ir(111) and (b) Gr/1 ML Fe/Ir(111); (c) atomic
structure as deduced by DFT, top view of the moiré pattern of Gr/Fe/Ir(111)
with C atoms represented in gray, Fe in red and Ir in cream; (d-e) charge
difference with respect to the free standing Gr and the pristine Fe/Ir(111)
slab, computed with DFT: (d) top view of the charge difference isosurfaces
(positive shown in purple, negative in yellow) and (e) side view of the
Gr/Fe/Ir(111) structure showing the graphene corrugation and charge difference
isosurfaces, using the same colors described for the previous panels.

Since LDA is known to underestimate the values of the orbital magnetic moments
in transition metals Truhlar09 , we have adopted a DFT+U scheme
cococcioni2005linear , with a Hubbard U parameter of 2 eV. In order to compare
calculated band structures with ARPES data, we applied an unfolding procedure
Popescu2012 ; unfold-x  to the computed bands of Gr/Fe/Ir(111) and Gr/Ir(111).
With this procedure, the band structure computed for the 10$\times$10
supercell is mapped into the graphene 1$\times$1 Brillouin zone by using the
unfold-x code unfold-x . In this way we obtain an effective band structure
corresponding to the graphene unit cell. In this picture, the
$\mathbf{k}$-dispersion is broadened by the break of the 1$\times$1
translational symmetry induced by the 10$\times$10 moiré pattern.

##  III Results and Discussion

###  III.1 Intercalation and structural properties

Figure 2:  LEED patterns (electron beam energy 140 eV) and zoomed pattern
around the (10) diffraction spot for Gr/7ML Fe/Ir(111) (a), Gr/0.5ML
Fe/Ir(111) (b) and Gr/Ir(111) (c); below: sketch of the main symmetry
directions and of the diffraction pattern of the Ir(111) surface (red dots);
the Surface Brillouin Zone (SBZ) is shown in black. Photoemission Energy
Distribution curves (EDC) taken at 40.814 eV (He IIα); angular integrated
spectra around $\Gamma$ (d) and $K$ (e) points of the SBZ for different
thickness of the intercalated Fe layer between Gr and Ir

A single layer of iron, intercalated under Gr grown on Ir(111), induces a
corrugation of the Gr membrane and preserves the periodicity of the moiré
superstructure, superimposed to the hexagonal Gr lattice, as unraveled by the
LEED patterns reported in Fig. 1(a,b). The diffraction pattern of Gr/Ir(111)
(Fig. 1a) shows bright spots in a hexagonal pattern surrounded by satellites,
consistent with the moiré superstructure caused by the lattice mismatch
between Gr and Ir NDiaye08 ; Scardamaglia2011 . After Fe intercalation, the
pattern is only slightly attenuated (Fig. 1b), thus the first Fe layer appears
commensurate to the Ir(111) surface lattice. This evidence is corroborated by
a recent work Decker_JPhys_2014  where STM measurements show that a few
monolayers (ML) of Fe intercalated under graphene arrange in registry with the
hexagonal Ir(111) surface, with a corrugation of 1.3 Å, larger than the one
measured for Gr/Ir(111) PhysRevLett.107.036101 .

We have further investigated the structural properties of this system by means
of DFT at the LDA+U level. The Gr/Fe/Ir(111) heterostructure was modeled as
described in Sec. II. The calculations confirm the large corrugation of the
graphene layer (1.44 Å), with a minimum graphene-Fe interplanar distance of
1.85 Å and a maximum of 3.29 Å on the crests of the hills, as shown in Fig.
1(c). Previous DFT calculations Decker_JPhys_2014  done at the PBE level
including van der Waals interactions by using the DFT–D2 method Grimme2006JCC
; Decker_D2_note  found a similar scenario, with slightly larger Gr-Fe
distances (2.05 and 3.33 Å for valleys and hills, respectively).

The topography of the Gr/Fe/Ir(111) heterostructure is similar to the
corrugated moiré superstructure observed for the Gr/Co/Ir(111) system
Avvisati_NanoLett_2018 , but very different from the structure reported for Gr
directly grown on bcc-Fe(110) Vinogradov_2012 . Previous STM measurements and
DFT calculations performed at the LDA level, for Gr/bcc-Fe(110) point out the
formation of a periodic corrugated pattern of the graphene layer parallel to
the [001] direction of the substrate, consisting in a supercell of 7$\times$17
graphene hexagons with a smaller corrugation of 0.6/0.9 Å, and only a small
fraction of the C atoms considerably elevated over the Fe surface, thus making
the entire graphene membrane strongly interacting with the metal substrate.

In the present case, the larger graphene corrugation modulates the Gr/Fe
interaction. This is illustrated by the charge difference computed for
Gr/Fe/Ir(111) shown in Fig. 1(d,e). The excess of negative charge (yellow
isosurface) is accumulated in the graphene membrane, donated by the Fe
intercalated layer. The redistribution of charge is more pronounced in the
valley regions and milder in the areas corresponding to the crests of the
hills, corroborating a different strength of hybridization between graphene
and Fe going from the valleys to the crests due to the increasing graphene–Fe
distance. The periodicity of strongly and weakly bound regions in which
covalent and van der Waals bonding dominate, respectively, induces a different
balance in the charge transfer determined also by the registry with respect to
the underlying Fe atoms.

Another signature of the similarity of the registry and graphene corrugation
resulting from Co and Fe intercalated layers can be deduced by comparing the C
1$s$ core level photoemission data for Gr/Fe/Ir(111) (as reported in
Supplemental Material supp-info ) and Gr/Co/Ir(111). In fact, there is a
multi-component C 1$s$ line-shape for both systems, in which the two main
features are assigned to C atoms in the Gr layer weakly (on the crests) and
strongly (in the valleys) bound to the Co and Fe layers Pacile_PRB_2014 ;
Avvisati_JPCC_2017 ; Avvisati_ASS_2020 . In contrast, when graphene is grown
directly on the bcc-Fe(110) surface, a single C 1$s$ core level component is
present, explained by the presence of extended regions of graphene coupled to
the Fe(110) substrate, even at the crests of the moiré superstructure
Vinogradov_2012 ; Vinogradov_2017 . For the graphene membrane directly grown
on Fe(110), a small feature in the C 1s core level has been detected at lower
binding energy (BE), linked to the presence of Fe carbide resulting from the
Fe-C solubility due to the high temperature reached during the chemical vapor
deposition growth procedure Vinogradov_2012 ; Vinogradov_2017 . The similarity
of the Fe intercalation process under Gr/Ir(111) with the Gr/Co/Ir(111)
heterostructure and the absence of a C 1$s$ component due to the formation of
iron carbide (see details in Supplemental Material supp-info ) demonstrate the
successful strategy to intercalate Fe under graphene at low temperature,
preventing any solubility into C.

The corrugation of the Gr membrane discussed so far is attenuated when the
number of Fe intercalated layers increases above a few ML, with a reduction of
the brightness of the superstructure LEED spots, as reported in Fig. 2(a-c).
At 7 ML, we observe the formation of a Fe film with a lattice constant smaller
than Ir(111) and therefore closer to the graphene structural parameters.
However, even when the strain induced by the Ir(111) surface is released, the
Fe film remains slightly incommensurate with graphene, as shown in Fig. 2(a).
This is in contrast with the case of Co intercalation beneath graphene on
Ir(111): when the Co film is formed, it recovers its hpc lattice constant and
becomes commensurate with the graphene layer Pacile_PRB_2014 ;
Avvisati_JPCC_2017 . Indeed, while the hcp Co(0001) surface lattice parameter
(2.50 Å), is comparable with the 2.46 Å of Gr, none of the bcc-Fe faces is
lattice-matched with graphene. In fact, graphene directly grown on bcc-Fe(110)
presents a distorted hexagonal symmetry with only a partial match with the
graphene lattice, as deduced by STM Vinogradov_2012 . In our case, the Fe
intercalated multilayer has a novel strained structural configuration with a
8% mismatch with respect to the graphene lattice constant (see details of the
LEED spot intensity analysis in the Supplemental Material supp-info ).

Importantly, the magnetic response of Fe films is extremely sensitive to tiny
variations of the structural configuration Moruzzi1986 , and the strained
lattice of the Fe film intercalated under graphene can give rise to an altered
distribution of the electronic majority and minority states and to different
spin and orbital configurations with respect to the bulk reference
Moulas_PRB_2008 .

Gr/Ir(111) Gr/Fe/Ir(111)

Figure 3:  Angular resolved photoelectron spectroscopy h$\nu$=40.814 eV
spectra of Gr/Ir(111) (left panels) and Gr/1 ML Fe/Ir(111) (right panels),
around the K point of the SBZ. Electronic band structure vs. k// as intensity
plot (upper panels); spectral density plotted as sequential EDC curves from
the K point towards $\Gamma$ (bottom panels).

###  III.2 Electronic structure

Gr/Ir(111) Gr/Fe/Ir(111)

Figure 4:  Band Structure and projected Density of States (pDOS) computed
within DFT for Gr/Ir(111) and Gr/1 ML Fe/ Ir(111): a) Gr/Ir(111) DOS projected
on C and Ir atomic orbitals and b) Gr/Ir(111) band structure, unfolded on the
1$\times$1 graphene unit cell as described in the main text; c) Gr/Fe/Ir(111)
majority spin DOS projected on C, Fe and Ir atomic orbitals; d) majority and
e) minority spin band structure, unfolded on the 1$\times$1 graphene unit
cell; f) Gr/Fe/Ir(111) minority spin pDOS.

A deeper insight into the interaction between Gr and the Fe intercalated
layer(s) can be unveiled by the electronic spectral density deduced by angular
resolved photoelectron spectroscopy of the Gr/Fe/Ir(111) heterostructure,
compared with ab initio theoretical predictions of the electronic density of
states (DOS) and band structure calculated by DFT.

Experimental data. The photoemission spectral density as a function of Fe
intercalation thickness is reported in Fig. 2 (d) and (e), at the $\Gamma$ and
K points of the surface Brillouin zone (SBZ), respectively. The Gr/Ir(111)
valence band at the $\Gamma$ and K points presents the expected electronic
spectral density PhysRevLett.102.056808 ; PhysRevB.84.075427  whereas the
intercalated system shows three new features when compared with the bare
Gr/Ir(111) : ($i$) the appearance of a peak at 10 eV and the disappearance of
the one related to the $\pi$\- bands of bare Gr/Ir(111) at 8 eV at the
$\Gamma$ point; ($ii$) the quenching of the intensity from the Ir surface
states close to the Fermi level; and ($iii$) the emergence of extra spectral
density in the low binding energy region.

Concerning ($i$), the intensity of the peak at 8 eV BE at $\Gamma$,
corresponding to the bottom of the $\pi$-band for bare Gr/Ir(111), decreases
with increasing Fe thickness. The peak emerging at 10.0 eV BE can be
associated to the shifted $\pi$-band, due to the interaction of Gr with the
intercalated Fe layer. An interacting $\pi$-band has also been observed for Co
and Ni layers sandwiched between Gr and Ir Pacile_PRB_2014 ; Vita_PRB_2014 ;
Pacile2013 ; massimi2014 . The presence of two features associated with the
$\pi$-bands, with opposite intensity behavior before the completion of the
first Fe ML, demonstrates the coexistence of bare Gr/Ir(111) regions and
regions of intercalated Fe atoms, up to the formation of a smooth Fe single
layer, when the peak at 8 eV (bare Gr/Ir(111)) is eventually quenched.

Furthermore, ($ii$) the progressive intensity lowering of the Ir(111) surface
states in the energy region 0-3 eV BE, upon Fe intercalation, without any
energy shift, suggests the absence of Fe-Ir intermixing observed at higher
intercalation temperatures Bazarnik_SurfSci_2015 ; Brede2016 . The emergence
of a spectral density ($iii$) at low binding energy close to the Fermi level
can be ascribed to electronic states mainly localized in the Fe layer, since
their intensity definitely raises at increasing Fe thickness, as can be
clearly observed in the spectral density of 7 ML Fe intercalated under
graphene (orange curve in Fig. 2e).

A clearer assignment of the electronic states due to the Fe-Gr interaction
process can be derived from the electronic state dispersion. The ARPES data
around the $K$ point of the SBZ for Gr/Ir(111) and Gr/1 ML Fe/Ir(111), are
shown in Fig. 3 (left and right panels, respectively). Considering the energy
region closer to the Fermi level, three localized Fe states can be identified
at about 0.2 eV, 0.5 eV and 1.5 eV BE, and found to be only slightly
dispersing over the whole SBZ (Fig. 3 right panel), as expected for a confined
Fe film. Inspecting the spectra at higher binding energies, other Fe-related
states are observed in the energy region of 3-4 eV BE, resonant with the
graphene $\pi$-band states. Besides the presence of these localized Fe states
close to the Fermi level, the most evident consequence of Fe intercalation
under graphene is the downshift of the Dirac cone at the $K$ point of the SBZ,
similarly to the case of Gr/Co/Ir(111) Pacile_PRB_2014 . In fact, the
$\pi$-band, as deduced by the EDCs (Fig. 3), is broadened and less defined
after intercalation and there is a coexistence of other Fe related states in
this energy region (3-4 eV). The position of the Dirac cone vertex (at K) is
downshifted by 3 eV, while the bottom of the $\pi$ band is less shifted at
$\Gamma$. It is worth noting that the Dirac cone is shifted by 2 eV towards
higher binding energy also for graphene directly grown on Fe(110)
Varykhalov2015 , suggesting a weaker interaction in that case.

Theory. The measured ARPES data are complemented by DFT calculations, as shown
in Fig. 4, where the electronic structures computed for Gr/Ir(111) and
Gr/Fe/Ir(111) are unfolded and mapped into the 1$\times$1 graphene Brillouin
zone along the $\Gamma$ – $K$ direction, as described in Sec. II. Projected
DOS (and projected bands in Supplemental Material supp-info ) are also
provided to complement the band information. The bands obtained for Gr/Ir(111)
are in good agreement with existing literature, as from e.g. Ref.
Voloshina2015JPCM . Concerning the Fe intercalated system, as also observed in
the experiments, the Gr $\pi$\- and $\sigma$-bands are shifted to higher
binding energies by the effect of Fe intercalation ($i$). While for Gr/Ir(111)
the $\pi$ bands are clearly recognizable, in the case of Gr/Fe/Ir(111) they
are strongly hybridized with Fe states above -6 eV and the upper part of the
cone is identifiable only for the minority spin states. The bottom of the
$\pi$ bands at $\Gamma$ moves from about -7.5 eV in Gr/Ir(111) to about -10.0
eV in Gr/Fe/Ir(111) (see Fig. S4 in Supplemental Material supp-info ), with an
overall downshift of $\sim$2.5 eV, in quite good agreement with the
experimental findings. Similarly, the $\sigma$ bands of graphene undergo a
downshift by about -1.5 eV at both $\Gamma$ (going from -2.5 to -4 eV) and $K$
(from -10.5 eV to -12 eV). While the experimental offset of these features is
found at larger binding energies (namely 4.0 eV and 5.5 eV as marked by
$\sigma$ and $\sigma^{\prime}$ in Fig. 2d), the downshift is quite consistent
with the calculations note_ARPES_DFT .

In order to identify Fe-related states, we consider the projected DOS, shown
in Fig. 4(c,f), where a strong spin splitting of the Fe components is observed
(see discussion below). For instance, a large peak of the Fe pDOS is found in
the energy region of 3-4 eV BE, in good agreement with the experimental data,
and mainly located in the majority spin channel. Such peak then overlaps with
the $\pi$ bands of Gr disrupting the Dirac cone for majority states. In
contrast, the Dirac cone is still faintly visible in the minority spin bands
at about 2.5-3 eV. In both spin channels the pDOS reveals overlapping peaks
due to Fe and graphene states, suggesting that the shadowing of the vertex of
the Dirac cone observed in the experimental photoemission data is clearly
induced by the hybridization of the graphene $\pi$-states with the Fe $d$
majority states. Furthermore, the spin resolved DOS in the energy region
closer to the Fermi level shows that the main Fe-related peak of the minority
spin states lies above the Fermi level, with smaller peaks between -2 eV and
the Fermi energy, in good agreement with the experimental observation ($iii$).
Overall, this picture is further confirmed by the projected unfolded bands
provided in Fig. S4 of Supplemental Material supp-info .

###  III.3 Magnetic properties and spin and orbital configuration of the Fe
layer

Figure 5:  XMCD spectra of Fe L2,3 absorption edges for Gr/1.4 ML Fe/Ir(111)
(top) and Gr/7ML Fe/Ir(111) (bottom) acquired in remanence at RT in normal
(left) and grazing (right) incidence geometries (see sketches for the
experimental geometry). The easy magnetization direction switching from
perpendicular (top) to parallel (bottom) to the surface plane is indicated by
the blue arrows.

In contrast with Fe grown on Ir(111), for which the loss of inversion symmetry
at the interface of the magnetic layer and substrate stabilizes different
skyrmion lattices depending on the Fe/Ir(111) stacking Heinze_NatPhys_2011 ;
vonBergmann2007 ; vonBergmann2006 ; vonBergmann2015 ; Grenz_PRL_2017 , when a
Fe monolayer is intercalated between Gr and Ir(111), the heterostructure
exhibits a ferromagnetic order with an out-of-plane easy magnetization axis
Decker_JPhys_2014 . For this reason, and based on the experimental evidence
described below, in the present DFT calculations we have considered only
collinear magnetic configurations, excluding in this way skyrmions or spin
spiral textures.

The pDOS computed for the Gr/Fe/Ir(111) heterostructure shows a large spin
split of the Fe states, with the maxima of the two spin pDOS almost 5 eV
apart. In particular, the majority spin states are almost fully occupied while
the largest peak of the minority spin states is empty, which results in an
average Fe-magnetic moment of 2.2 $\mu_{B}/$atom. The computed ground state
magnetic configuration is ferromagnetic even if the Fe magnetic moment shows a
modulation over different sites ranging from 2.0 to 2.7 $\mu_{B}/$atom, in
good agreement with previous calculation Brede2014 . The modulations is
determined by the graphene layer, with the Fe atoms located below the graphene
crests, i.e. with larger Fe-C distances, having the largest values. Concerning
graphene, the C magnetic moments are non-zero but quite small (at the LDA+U
level ranging from -0.02 to 0.01 $\mu_{B}/$atom), similarly to what reported
for Gr/Co/Ir(111) Vita_PRB_2014 ; Avvisati_NanoLett_2018 , with a distance-
and sublattice-dependent modulation.

Overall, we find that the Fe monolayer displays a clear spin splitting, as an
effect of the band narrowing and subsequent increased number of majority
electrons. This is confirmed by the strong localization of the Fe-induced
electronic states observed in the photoemission spectral density. It is
interesting to point out that a different scenario is found for the
Gr/Co/Ir(111) interface Vita_PRB_2014  (see Supplemental Material supp-info
for a detailed comparison with Gr/Co/Ir(111) pDOS) where Co spectral weight in
the minority channel is shifted at lower energies leading to a spin split
significantly smaller than in the present case.

The strong ferromagnetic behavior is confirmed by the X ray absorption spectra
and the XMCD at the L2,3 edges for 1.4 ML and 7 ML of Fe sandwiched between Gr
and Ir(111), shown in Fig. 5. XAS are acquired in remanence conditions with
circularly polarized radiation and with photons impinging normal (left panel)
and grazing (right panel) to the surface. The XMCD spectra are obtained as the
difference between the absorption edges acquired with left- and right-
circularly polarized radiation. The higher dichroic response with photon
impinging at normal incidence (left panel) reveals an out-of-plane magnetic
anisotropy of the Fe layer when it is pseudomorphic to the Ir(111) surface,
with a stretched Fe-Fe distance with respect the bulk bcc(110) or (111)
surfaces. The higher dichroic response for Fe L2,3 XMCD at grazing incidence
(right panel) for the thicker intercalated Fe film, unveils a switch of the
easy magnetization axis from perpendicular to the surface to in-plane. This is
consistent with the magnetic response of heterostructures of single Co or
FexCo1-x layers on Ir(111) covered with a Gr membrane Avvisati_ASS_2020 ;
Avvisati_NanoLett_2018 .

The spin and the orbital moments at the Fe site, as deduced via XMCD sum
rules, show a doubled $L/S_{\text{eff}}$ ratio ($S_{\text{eff}}=S+7D$ and D is
the dipolar moment Thole_PRL_1992 ; Carra_PRL_1993 ), with respect to the bulk
element (0.11$\pm 0.01$ for a Fe single layer and 0.09 $\pm 0.01$ for the 7
ML). These experimental results, related only to the intensity of the dichroic
signal and independent on the number of $3d$ holes, suggest an enhanced
magnetic response due to the redistribution of the majority spin states in the
conduction Fe bands. The evaluated average total moment is 2.1$\pm 0.2$
$\mu_{B}/\text{atom}$ for 7 ML of Fe intercalated under Gr, in fair agreement
with 2.2 $\mu_{B}/$atom (spin moment), as deduced by the DFT predictions. The
orbital and spin moment for a single layer of Fe, where interface effects can
play a role, shows a comparable magnetic response, as deduced by the similar
dichroic signals reported in Fig. 5. As mentioned above, the pseudomorphic
hexagonal Fe single layer grown directly on Ir(111) presents a complex
magnetic structure with skyrmionic spin spiral textures stabilized by the
$3d–5d$ hybridization between Fe and Ir. vonBergmann2006 ; vonBergmann2007 ;
vonBergmann2015 ; Grenz_PRL_2017 . On the other hand, other magnetic materials
like Europium have been successfully intercalated underneath graphene on the
same Ir(111) surface doi:10.1021/nl402797j  and, depending on its coverage, Eu
displays either a paramagnetic or a ferromagnetic behavior PhysRevB.90.235437
. The clear magnetic dichroism response of the single Fe layer intercalated
under graphene suggests a different spin configuration. Furthermore, this
confined single Fe intercalated layer, with stretched Fe-Fe distances,
presents spin and orbital configurations similar to those of the thin Fe
intercalated film, where the influence of the Fe-Ir interface is completely
released.

##  IV Conclusions

The strategy to sandwich the Fe layer beneath the graphene membrane with an
intercalation process at moderate temperature (500 K), prevents any alloying,
and the absence of any hallmark of Fe-Ir and Fe-C intermixing proves the
formation of a well-defined homogeneous Fe layer in registry with the Ir(111)
surface. Such single layer of Fe, protected by the graphene membrane, induces
a downshift in energy and a symmetry breaking of the Dirac cone due to the
interaction between Gr and Fe majority states resonant in the energy region of
the vertex of the cone. The redistribution of the spin resolved Fe pDOS with a
narrowing of the Fe bands and a larger spin splitting between majority (fully
occupied) and minority states (almost empty) of the Fe states associated with
increased total Fe magnetic moment influences the magnetic response of the Fe
intercalated layer(s).

In contrast to the case of direct growth of Fe on Ir(111), where, despite the
large Fe magnetic moment, a non-collinear magnetic order has been observed
Heinze_NatPhys_2011 ; vonBergmann2007 ; vonBergmann2006 ; vonBergmann2015 ,
the Fe layer with graphene on top reveals a ferromagnetic long range order
with spin and orbital moments that are higher than the ones found for the bulk
phases. The graphene top layer acts not only as a protective membrane, but
also allows for a stable ferromagnetic configuration, counteracting the effect
of Ir substrate. The concomitant dimensionality reduction with a narrowing of
the $d$ bands and a reduced superimposition between the spin-split majority
and minority Fe bands, further contribute to the transition of a single Fe
layer (or few layers of Fe), from weak to strong ferromagnet, when
intercalated beneath graphene. These 3$d$ confined layers protected by a
graphene membrane with a novel structural configuration with respect to the
bulk lattice arrangement offer a powerful playground to tune their electronic
structure and magnetic state for magnetic/spintronic devices.

##  V Acknowledgements

The authors thank the experimental assistance of BOREAS beamline staff at
ALBA, of SuperESCA beamline staff at Elettra, and of Alessio Vegliante in
LOTUS lab in Rome. G.A. and MG.B. acknowledge support from translational
access CalipsoPlus funding. C.C., D.V., A.F., D.A.L. acknowledge partially
support from the MaX – MAterials design at the eXascale – Centre of
Excellence, funded by the European Union program H2020-INFRAEDI-2018-1 (Grant
No. 824143). They also thank SUPER (Supercomputing Unified Platform – Emilia-
Romagna) from Emilia-Romagna POR-FESR 2014-2020 regional funds. The work was
partially supported by PRIN FERMAT (2017KFY7XF) from Italian Ministry MIUR and
by Sapienza Ateneo funds. Computational time on the Marconi100 machine at
CINECA was provided by the Italian ISCRA program.

## References

  * (1) N. Tombros, C. Jozsa, M. Popinciuc, H. T. Jonkman, and B. J. van Wees,  Nature 448, 571 (2007).
  * (2) W. Han, R. Kawakami, M. Gmitra, and J. Fabian,  Nature Nanotechnology 9, 794 (2014).
  * (3) D. Usachov, A. Fedorov, M. M. Otrokov, A. Chikina, O. Vilkov, A. Petukhov, A. G. Rybkin, Y. M. Koroteev, E. V. Chulkov, V. K. Adamchuk, A. Grüneis, C. Laubschat, and D. V. Vyalikh,  Nano Lett. 15, 2396 (2015),  PMID: 25734657.
  * (4) L. V. Dzemiantsova, M. Karolak, F. Lofink, A. Kubetzka, B. Sachs, K. von Bergmann, S. Hankemeier, T. O. Wehling, R. Frömter, H. P. Oepen, A. I. Lichtenstein, and R. Wiesendanger,  Phys. Rev. B 84, 205431 (2011).
  * (5) D. Pacilè, P. Leicht, M. Papagno, P. Sheverdyaeva, P. Moras, C. Carbone, K. Krausert, L. Zielke, M. Fonin, Y. Dedkov, F. Mittendorfer, J. Doppler, A. Garhofer, and J. Redinger,  Phys. Rev. B 87, 035420 (2013).
  * (6) A. Varykhalov, D. Marchenko, J. Sánchez-Barriga, M. R. Scholz, B. Verberck, B. Trauzettel, T. O. Wehling, C. Carbone, and O. Rader,  Phys. Rev. X 2, 041017 (2012).
  * (7) L. Massimi, S. Lisi, D. Pacilè, C. Mariani, and M. G. Betti,  Beilstein J. Nanotech. 5, 308 (2014).
  * (8) R. Decker, J. Brede, N. Atodiresei, V. Caciuc, S. Blügel, and R. Wiesendanger,  Phys. Rev. B 87, 041403 (2013).
  * (9) D. Pacilè, S. Lisi, I. Di Bernardo, M. Papagno, L. Ferrari, M. Pisarra, M. Caputo, S. K. Mahatha, P. M. Sheverdyaeva, P. Moras, P. Lacovig, S. Lizzit, A. Baraldi, M. G. Betti, and C. Carbone,  Phys. Rev. B 90, 195446 (2014).
  * (10) H. Vita, S. Böttcher, P. Leicht, K. Horn, A. B. Shick, and F. Máca,  Phys. Rev. B 90, 165432 (2014).
  * (11) M. Weser, Y. Rehder, K. Horn, M. Sicot, M. Fonin, A. B. Preobrajenski, E. N. Voloshina, E. Goering, and Y. S. Dedkov,  Appl. Phys. Lett. 96, 012504 (2010).
  * (12) R. Decker, M. Bazarnik, N. Atodiresei, V. Caciuc, S. Blügel, and R. Wiesendanger,  J. Phys.: Condens. Mat. 26, 394004 (2014).
  * (13) J. Coraux, A. T. N’Diaye, N. Rougemaille, C. Vo-Van, A. Kimouche, H.-X. Yang, M. Chshiev, N. Bendiab, O. Fruchart, and A. K. Schmid,  J. Phys. Chem. Lett. 3, 2059 (2012).
  * (14) M. Cattelan, G. W. Peng, E. Cavaliere, L. Artiglia, A. Barinov, L. T. Roling, M. Favaro, I. Píš, S. Nappini, E. Magnano, F. Bondino, L. Gavioli, S. Agnoli, M. Mavrikakis, and G. Granozzi,  Nanoscale 7, 2450 (2015).
  * (15) A. Lodesani, A. Picone, A. Brambilla, D. Giannotti, M. S. Jagadeesh, A. Calloni, G. Bussetti, G. Berti, M. Zani, M. Finazzi, et al.,  ACS nano 13, 4361 (2019).
  * (16) N. A. Vinogradov, A. A. Zakharov, V. Kocevski, J. Rusz, K. A. Simonov, O. Eriksson, A. Mikkelsen, E. Lundgren, A. S. Vinogradov, N. Mårtensson, and A. B. Preobrajenski,  Phys. Rev. Lett. 109, 026101 (2012).
  * (17) N. Vinogradov, K. Simonov, A. Generalov, J. Drnec, F. Carlà, A. Vinogradov, A. Preobrajenski, N. Mårtensson, and R. Felici,  Carbon 111, 113 (2017).
  * (18) A. Varykhalov, J. Sánchez-Barriga, D. Marchenko, P. Hlawenka, P. S. Mandal, and O. Rader,  Nature Commun. 6, 7610 (2015).
  * (19) H. Okamoto,  Journal of Phase Equilibria 13, 543 (1992).
  * (20) S. Vlaic, N. Rougemaille, A. Artaud, V. Renard, L. Huder, J.-L. Rouviére, A. Kimouche, B. Santos, A. Locatelli, V. Guisset, et al.,  J. Phys. Chem. Lett. 9, 2523 (2018).
  * (21) G. Avvisati, P. Gargiani, D. Lizzit, M. Valvidares, P. Lacovig, C. Petrillo, F. Sacchetti, and M. G. Betti,  Applied Surface Science 527, 146599 (2020).
  * (22) P. Gargiani, R. Cuadrado, H. B. Vasili, M. Pruneda, and M. Valvidares,  Nat. Commun. 8, 699 (2017).
  * (23) M. Bazarnik, R. Decker, J. Brede, and R. Wiesendanger,  Surf. Sci. 639, 70 (2015).
  * (24) https://sites.google.com/uniroma1.it/nano-surface-physics/.
  * (25) A. Barla, J. Nicolás, D. Cocco, S. M. Valvidares, J. Herrero-Martín, P. Gargiani, J. Moldes, C. Ruget, E. Pellegrin, and S. Ferrer,  Journal of Synchrotron Radiation 23, 1507 (2016).
  * (26) G. Avvisati, C. Cardoso, D. Varsano, A. Ferretti, P. Gargiani, and M. G. Betti,  Nano Lett. 18, 2268 (2018).
  * (27) J. Brede, J. Sławińska, M. Abadia, C. Rogero, J. E. Ortega, I. Piquero-Zulaica, J. Lobo-Checa, A. Arnau, and J. I. Cerdá,  2D Materials 4, 015016 (2016).
  * (28) See Supplemental Material for the detailed description of the C $1s$, Ir $4f_{7/2}$, and Fe $2p_{3/2}$ core levels; the LEED spot intensity analysis for Gr/Ir (111), Gr/1ML Fe/Ir(111) and Gr/7ML Fe/Ir(111); the unfolded bands of Gr/1ML-Fe/Ir(111); the spin-resolved pDOS for Gr/1ML-Fe/Ir(111) and Gr/1ML-Co/Ir(111).
  * (29) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al.,  J. Phys.: Condens. Matter 21, 395502 (2009).
  * (30) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio Jr., A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni,  J. Phys.: Condens. Matter 29, 465901 (2017).
  * (31) J. P. Perdew and A. Zunger,  Phys. Rev. B 23, 5048 (1981).
  * (32) A. Calloni, G. Bussetti, G. Avvisati, M. S. Jagadeesh, D. Pacilè, A. Ferretti, D. Varsano, C. Cardoso, L. Duò, F. Ciccacci, and M. G. Betti,  J. Chem. Phys. 153, 214703 (2020).
  * (33) C. J. Cramer and D. G. Truhlar,  Phys. Chem. Chem. Phys. 11, 10757 (2009).
  * (34) M. Cococcioni and S. De Gironcoli,  Phys. Rev. B 71, 035105 (2005).
  * (35) V. Popescu and A. Zunger,  Phys. Rev. B 85, 085201 (2012).
  * (36) P. Bonfà,  https://bitbucket.org/bonfus/unfold-x.
  * (37) A. N’Diaye, J. Coraux, T. Plasa, C. Busse, and T. Michely,  New J. Phys. 10, 043033 (2008).
  * (38) M. Scardamaglia, G. Forte, S. Lizzit, A. Baraldi, P. Lacovig, R. Larciprete, C. Mariani, and M. G. Betti,  J. Nanopart. Res. 13, 6013 (2011).
  * (39) C. Busse, P. Lazić, R. Djemour, J. Coraux, T. Gerber, N. Atodiresei, V. Caciuc, R. Brako, A. T. N’Diaye, S. Blügel, J. Zegenhagen, and T. Michely,  Phys. Rev. Lett. 107, 036101 (2011).
  * (40) S. Grimme,  J. Comput. Chem. 27, 1787 (2006).
  * (41) Note that the DFT-D2 correction for van der Walls interactions has been reparametrized in these works.
  * (42) G. Avvisati, S. Lisi, P. Gargiani, A. Della Pia, O. De Luca, D. Pacilè, C. Cardoso, D. Varsano, D. Prezzi, A. Ferretti, and M. G. Betti,  J. Phys. Chem. C 121, 1639 (2017).
  * (43) V. L. Moruzzi, P. M. Marcus, K. Schwarz, and P. Mohn,  Phys. Rev. B 34, 1784 (1986).
  * (44) G. Moulas, A. Lehnert, S. Rusponi, J. Zabloudil, C. Etz, S. Ouazi, M. Etzkorn, P. Bencok, P. Gambardella, P. Weinberger, and H. Brune,  Phys. Rev. B 78, 214424 (2008).
  * (45) I. Pletikosić, M. Kralj, P. Pervan, R. Brako, J. Coraux, A. T. N’Diaye, C. Busse, and T. Michely,  Phys. Rev. Lett. 102, 056808 (2009).
  * (46) M. Kralj, I. Pletikosić, M. Petrović, P. Pervan, M. Milun, A. T. N’Diaye, C. Busse, T. Michely, J. Fujii, and I. Vobornik,  Phys. Rev. B 84, 075427 (2011).
  * (47) Y. Dedkov and E. Voloshina,  J. Phys.: Condens. Matter 27, 303002 (2015).
  * (48) Note that comparing ARPES data with KS-DFT calculations absolute quantitative agreement in the position of the features cannot be expected.
  * (49) S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel,  Nat. Phys. 7, 713 (2011).
  * (50) K. von Bergmann, S. Heinze, M. Bode, G. Bihlmayer, S. Blügel, and R. Wiesendanger,  New J. Phys. 9, 396 (2007).
  * (51) K. von Bergmann, S. Heinze, M. Bode, E. Y. Vedmedenko, G. Bihlmayer, S. Blügel, and R. Wiesendanger,  Phys. Rev. Lett. 96, 167203 (2006).
  * (52) K. von Bergmann, M. Menzel, A. Kubetzka, and R. Wiesendanger,  Nano Lett. 15, 3280 (2015).
  * (53) J. Grenz, A. Köhler, A. Schwarz, and R. Wiesendanger,  Phys. Rev. Lett. 119, 047205 (2017).
  * (54) J. Brede, N. Atodiresei, V. Caciuc, M. Bazarnik, A. Al-Zubi, S. Blügel, and R. Wiesendanger,  Nature Nanotechnology 9, 1018 (2014).
  * (55) B. T. Thole, P. Carra, F. Sette, and G. van der Laan,  Phys. Rev. Lett. 68, 1943 (1992).
  * (56) P. Carra, B. T. Thole, M. Altarelli, and X. Wang,  Phys. Rev. Lett. 70, 694 (1993).
  * (57) S. Schumacher, T. O. Wehling, P. Lazić, S. Runte, D. F. Förster, C. Busse, M. Petrović, M. Kralj, S. Blügel, N. Atodiresei, V. Caciuc, and T. Michely,  Nano Letters 13, 5013 (2013).
  * (58) S. Schumacher, F. Huttmann, M. Petrović, C. Witt, D. F. Förster, C. Vo-Van, J. Coraux, A. J. Martínez-Galera, V. Sessi, I. Vergara, R. Rückamp, M. Grüninger, N. Schleheck, F. Meyer zu Heringdorf, P. Ohresser, M. Kralj, T. O. Wehling, and T. Michely,  Phys. Rev. B 90, 235437 (2014).


# Binding of curvature-inducing proteins onto tethered vesicles

Hiroshi Noguchi<EMAIL_ADDRESS>Institute for Solid State Physics,
University of Tokyo, Kashiwa, Chiba 277-8581, Japan

###### Abstract

A tethered vesicle, which consists of a cylindrical membrane tube and a
spherical vesicle, is produced by a mechanical force that is experimentally
imposed by optical tweezers and a micropipette. This tethered vesicle is
employed for examining the curvature sensing of curvature-inducing proteins.
In this study, we clarify how the binding of proteins with a laterally
isotropic spontaneous curvature senses and generates the membrane curvatures
of the tethered vesicle using mean-field theory and meshless membrane
simulation. The force-dependence curves of the protein density in the membrane
tube and the tube curvature are reflection symmetric and point symmetric,
respectively, from the force point, in which the tube has a sensing curvature.
The bending rigidity and spontaneous curvature of the bound proteins can be
estimated from these force-dependence curves. First-order transitions can
occur between low and high protein densities in the tube at both low and high
force amplitudes. The simulation results of the homogeneous phases agree very
well with the theoretical predictions. In addition, beaded-necklace-like tubes
with microphase separation are found in the simulation.

##  I Introduction

In living cells, numerous protein types work together to regulate biomembrane
shapes McMahon and Gallop (2005); Suetsugu  _et al._ (2014); Johannes  _et
al._ (2015); Brandizzi and Barlowe (2013); Hurley  _et al._ (2010); McMahon
and Boucrot (2011). Proteins are also involved in various dynamical processes,
such as end/exocytosis and vesicle transport. The Bin/Amphiphysin/Rvs (BAR)
superfamily proteins bend the membrane along its axis and generate cylindrical
membrane tubes McMahon and Gallop (2005); Suetsugu  _et al._ (2014); Johannes
_et al._ (2015); Itoh and De Camilli (2006); Mim and Unger (2012). Other
proteins, such as clathrin and coat protein complex (COPI and COPII), generate
spherical buds Johannes  _et al._ (2015); Brandizzi and Barlowe (2013); Hurley
_et al._ (2010); McMahon and Boucrot (2011). Thus, understanding the mechanism
of these curvature generations is important.

These curvature-inducing proteins are known to sense the membrane curvature
and are concentrated in membranes which have their preferred curvatures.
Various types of proteins can be examined using a tethered vesicle pulled by
optical tweezers and a micropipette Dimova (2014); Baumgart  _et al._ (2011);
Sorre  _et al._ (2012); Prévost  _et al._ (2015); Roux  _et al._ (2010);
Rosholm  _et al._ (2017); Allard  _et al._ (2020). With increasing force or
length, the vesicle first deforms into a lemon shape and subsequently forms a
narrow membrane tube (tether) protruding from a spherical vesicle Hotani  _et
al._ (1999); Inaba  _et al._ (2005); Wu  _et al._ (2016). Moreover, an
elongational force can be produced by the growth of protein filaments in vitro
Hotani  _et al._ (1999); Gavriljuk  _et al._ (2021) and in vivo Svitkina
(2018); Gallop (2020). Curvature-inducing proteins typically bind more onto
the membrane tube than the remaining spherical component. BAR proteins
Baumgart  _et al._ (2011); Sorre  _et al._ (2012); Prévost  _et al._ (2015),
dynamin Roux  _et al._ (2010), and G-protein coupled receptors Rosholm  _et
al._ (2017) have been reported to exhibit curvature sensing.

The aim of this study is to understand the curvature sensing and generation of
the curvature-inducing proteins with an isotropically spontaneous curvature on
a tethered vesicle. We employ mean-field theory and meshless membrane
simulation. In mean-field theory, a simplified geometry is considered for
vesicles. Previously, we used a vesicle consisting of many spherical
components to study budding Noguchi (2021). Here, we apply the same scheme to
the tethered vesicle. Although mean-field theories have been used to analyze
the experimental results of tethered vesicles Prévost  _et al._ (2015);
Rosholm  _et al._ (2017), they have been applied in narrow ranges of the
parameters, and the curvature of the spherical component has not been
considered. Here, we systematically investigate the protein binding onto the
tethered vesicle over a wide range of parameters.

Several types of membrane models have been developed for coarse-grained
simulations Müller  _et al._ (2006); Venturoli  _et al._ (2006); Noguchi
(2009). For a large-scale simulation, we developed two types of meshless
membrane models Noguchi and Gompper (2006); Shiba and Noguchi (2011), in which
membrane particles self-assemble into a single-layer membrane and the
mechanical properties can be varied over a wide range. Here, we employ a spin
meshless membrane model Shiba and Noguchi (2011), since it can vary the
spontaneous curvature and has been applied to membrane deformation by
curvature-inducing proteins with an isotropic spontaneous curvature Noguchi
(2016a); Goutaland  _et al._ (2021) and with an anisotropic spontaneous
curvature Noguchi (2014, 2015, 2016b); Noguchi and Fournier (2017); Noguchi
(2017, 2019a), as well as topological changes of membranes Noguchi (2019b, c).
In mean-field theory, we assume a uniform distribution of the bound proteins
in each membrane component. However, phase separation has been obtained in
simulations of membrane tubes Noguchi (2016a) and flat membranes Goutaland
_et al._ (2021). We clarify where the phase separation of the bound proteins
occurs in the membrane tubes under a constant external force and where the
assumption of uniform protein density is valid.

The mean-field theory of the tethered vesicle is described in Sec. II.
Simulations of membrane tubes are described in Sec. III. The simulation
results are compared with the theoretical results in Sec. III.2. Finally, a
summary and discussion are presented in Sec. IV.

##  II Tethered vesicle

###  II.1 Mean-field theory

A vesicle consists of a sphere with radius $R_{\rm sp}$ and a cylinder with
radius $R_{\rm cy}$ and length $L_{\rm cy}$ as depicted in Fig. 1. Here, a
thin tube of $R_{\rm cy}\ll R_{\rm sp}$ is considered, such that the end and
foot parts of the cylindrical tube are neglected. This geometry was previously
used in Ref. 38. In experiments, the end shape of the tube can be strongly
affected by force-imposing methods such as optical tweezers, whereas the foot
has a catenoid-like shape with a low mean-curvature Powers  _et al._ (2002).

Figure 1:  Simplified geometry of a tethered vesicle. The vesicle consists of
a sphere with a radius $R_{\rm sp}$ and a cylinder with a radius $R_{\rm cy}$
and length $R_{\rm cy}$.

The total surface area $A$ and volume $V$ are given as follows:

$\displaystyle A$ $\displaystyle=$ $\displaystyle 4\pi R_{0}^{2}=4\pi R_{\rm
sp}^{2}+2\pi R_{\rm cy}L_{\rm cy},$ (1) $\displaystyle V$ $\displaystyle=$
$\displaystyle\frac{4\pi}{3}R_{\rm sp}^{3}+\pi R_{\rm cy}^{2}L_{\rm cy},$ (2)

where $R_{0}$ is the radius of a sphere with the same surface area. From eqn
(1) and (2), the radius and length of the membrane tube can be expressed as
follows:

$\displaystyle\frac{R_{\rm cy}L_{\rm cy}}{R_{0}^{2}}$ $\displaystyle=$
$\displaystyle 2(1-r^{2}),$ (3) $\displaystyle\frac{R_{\rm cy}}{R_{0}}$
$\displaystyle=$ $\displaystyle\frac{2(v_{\rm r}-r^{3})}{3(1-r^{2})},$ (4)

where $r=R_{\rm sp}/R_{0}$ and reduced volume $v_{\rm r}=V/(4\pi
R_{0}^{3}/3)$.

Curvature-inducing proteins bind to the membrane, depending on the local
membrane curvature. The bending free energy $F_{\rm cv}$ is given by

$F_{\rm cv}=\int dA\ \Big{\\{}2\kappa_{\rm d}H^{2}(1-\phi)+\frac{\kappa_{\rm
p}}{2}(2H-C_{0})^{2}\phi\Big{\\}},$ (5)

where $H$ is the mean curvature of the membrane ($H=(C_{1}+C_{2})/2$, where
$C_{1}$ and $C_{2}$ are the principal curvatures). The first and second terms
represent the bending energy of the bare (unbound) and protein-bound
membranes, respectively, and $\phi$ is the protein density ($\phi=1$ at the
maximum coverage). The unbound membrane has bending rigidity $\kappa_{\rm d}$
and zero spontaneous curvature; the bound membrane has a larger bending
rigidity $\kappa_{\rm p}$ and spontaneous curvature $C_{0}$. In this study,
the effects of the saddle-splay modulus $\bar{\kappa}$ are not considered.
Although $\bar{\kappa}$ may depend on the protein density, the cylindrical
membrane has zero Gaussian curvature ($C_{1}C_{2}$) such that it has no
influence in the middle of the tube. Small effects may appear through the
neglected regions (the end and foot of the membrane tube). More general
aspects of the protein-binding model are discussed in Ref. 22.

To induce a tether (thin membrane tube), an external force $f_{\rm ex}$ is
imposed. Here, it is assumed that the proteins are homogeneously distributed
in each membrane component (the spherical or cylindrical component). The
subscripts sp and cy represent the quantities of spherical and cylindrical
components of the vesicle, respectively. The free energy $F$ of the tethered
vesicle is given by

$\displaystyle F$ $\displaystyle=$ $\displaystyle F_{\rm sp}+F_{\rm cy}-f_{\rm
ex}(L_{\rm cy}+2R_{\rm sp}),$ (6) $\displaystyle\frac{F_{\rm sp}}{4\pi}$
$\displaystyle=$ $\displaystyle 2(\kappa_{\rm dif}\phi_{\rm sp}+\kappa_{\rm
d})-2\kappa_{\rm p}C_{0}R_{0}r\phi_{\rm sp}$
$\displaystyle+R_{0}^{2}r^{2}(\sigma_{\rm p}\phi_{\rm sp}+b\phi_{\rm sp}^{2})$
$\displaystyle+\frac{k_{\rm B}TR_{0}^{2}r^{2}}{a_{\rm p}}[\phi_{\rm
sp}\ln(\phi_{\rm sp})+(1-\phi_{\rm sp})\ln(1-\phi_{\rm sp})],$
$\displaystyle\frac{F_{\rm cy}}{\pi}$ $\displaystyle=$
$\displaystyle(\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm d})\frac{L_{\rm
cy}}{R_{\rm cy}}-2\kappa_{\rm p}C_{0}L_{\rm cy}\phi_{\rm cy}$
$\displaystyle+2R_{\rm cy}L_{\rm cy}\Big{\\{}(\sigma_{\rm p}\phi_{\rm
cy}+b\phi_{\rm cy}^{2})$ $\displaystyle+\frac{k_{\rm B}T}{a_{\rm p}}[\phi_{\rm
cy}\ln(\phi_{\rm cy})+(1-\phi_{\rm cy})\ln(1-\phi_{\rm cy})]\Big{\\}}$
$\displaystyle=$ $\displaystyle(\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm
d})\frac{9(1-r^{2})^{3}}{2(v_{\rm r}-r^{3})^{2}}-6\kappa_{\rm
p}C_{0}R_{0}\frac{(1-r^{2})^{2}\phi_{\rm cy}}{v_{\rm r}-r^{3}}$
$\displaystyle+4(1-r^{2})R_{0}^{2}\Big{\\{}(\sigma_{\rm p}\phi_{\rm
cy}+b\phi_{\rm cy}^{2})$ $\displaystyle+\frac{k_{\rm B}T}{a_{\rm p}}[\phi_{\rm
cy}\ln(\phi_{\rm cy})+(1-\phi_{\rm cy})\ln(1-\phi_{\rm cy})]\Big{\\}},$

where $\kappa_{\rm dif}=\kappa_{\rm p}-\kappa_{\rm d}$ and $\sigma_{\rm
p}=-\mu/a_{\rm p}+\kappa_{\rm p}C_{0}^{2}/2$. The chemical potential of the
protein binding is $\mu$, and $a_{\rm p}$ is the membrane area bound by one
protein. The last terms in Eqs. (II.1)–(II.1) represent the mixing entropy of
the bound proteins. The inter-protein interactions are taken into account as
squared density terms of $b\phi^{2}$. Proteins have repulsive or attractive
interactions at $b>0$ and $b<0$, respectively.

The densities $\phi_{\rm sp}$ and $\phi_{\rm cy}$ are obtained from $\partial
F/\partial\phi|_{H}=0$:

$\phi=\frac{1}{1+\exp\big{[}\frac{a_{\rm p}}{k_{\rm B}T}(2\kappa_{\rm
dif}H^{2}-2\kappa_{\rm p}C_{0}H+\sigma_{\rm p}+2b\phi)\big{]}},$ (10)

where $H=1/R_{\rm sp}$ and $H=1/2R_{\rm cy}$ for the spherical and cylindrical
components, respectively Noguchi (2021). For $b=0$, this is a sigmoid function
of $\mu$. For $b\neq 0$, Eq. (10) is iteratively solved with an updated $\phi$
value in the right hand side. Thus, the free energy $F$ is expressed as a
function of one variable $r$ using Eqs. (II.1) and (II.1) with $\phi_{\rm sp}$
and $\phi_{\rm sp}$ obtained using Eq. (10). Hence, the free-energy minimum is
calculated by $\partial F/\partial r=0$.

The maximum binding (sensing) of the proteins occurs at a higher density than
the curvature generation Noguchi (2021). The curvatures of sensing and
generation are obtained by $\partial\phi/\partial H=0$ and $\partial
F/\partial H=0$ in the absence of any constraints and external forces,
respectively; The former and latter are the preferred curvatures for protein
binding and for the entire membrane including bare membrane parts,
respectively. The sensing curvature is $H_{\rm s}=\kappa_{\rm
p}C_{0}/2\kappa_{\rm dif}$, and the curvature of the generation is $H_{\rm
g}=\kappa_{\rm p}\phi C_{0}/2(\kappa_{\rm dif}\phi+\kappa_{\rm d})$. Hence,
the maximum protein density $\phi_{\rm cy}^{\rm max}$ is obtained at $1/R_{\rm
cy}=C_{\rm s}=\kappa_{\rm p}C_{0}/\kappa_{\rm dif}$ for the cylindrical tube:
$\phi_{\rm cy}^{\rm max}=1/\\{1+\exp[-(\mu-\mu_{0})/k_{\rm B}T]\\}$ at $b=0$,
where $\mu_{0}=-a_{\rm p}\kappa_{\rm p}{C_{0}}^{2}(\kappa_{\rm p}/\kappa_{\rm
dif}-1)/2$. At $\mu=\mu_{0}$, $\phi_{\rm cy}^{\rm max}=1/2$ is obtained. In
the dilute limit $\phi_{\rm sp}\ll 1$ and $\phi_{\rm sp}\ll 1$ with $R_{\rm
sp}\gg R_{\rm cy}$ and $b=0$, the density ratio is given by

$\frac{\phi_{\rm cy}}{\phi_{\rm sp}}\simeq\exp\Big{\\{}-\frac{\kappa_{\rm
dif}a_{\rm p}}{2k_{\rm B}T}\Big{[}\Big{(}\frac{1}{R_{\rm cy}}-C_{\rm
s}\Big{)}^{2}-C_{\rm s}^{2}\Big{]}\Big{\\}},$ (11)

as reported in Ref. 12.

The surface tension $\sigma_{\rm d}$ and osmotic pressure $\Pi$ can be
expressed as Lagrange multipliers to maintain the area and volume,
respectively: $\breve{F}=F+\sigma_{\rm d}A-\Pi V$. Then, $\sigma_{\rm d}$,
$\Pi$, and $f_{\rm ex}$ satisfy $\partial\breve{F}/\partial R_{\rm
sp}|_{R_{\rm cy},L_{\rm cy}}=0$, $\partial\breve{F}/\partial R_{\rm
cy}|_{R_{\rm sp},L_{\rm cy}}=0$, and $\partial\breve{F}/\partial L_{\rm
cy}|_{R_{\rm sp},R_{\rm cy}}=0$:

$\displaystyle\Pi$ $\displaystyle=$ $\displaystyle\frac{f_{\rm ex}}{\pi R_{\rm
cy}^{2}}-\frac{2(\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm d})}{R_{\rm
cy}^{3}}+\frac{2\kappa_{\rm p}C_{0}\phi_{\rm cy}}{R_{\rm cy}^{2}},$ (12)
$\displaystyle\sigma_{\rm d}$ $\displaystyle=$ $\displaystyle\Pi R_{\rm
cy}+\frac{\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm d}}{2R_{\rm
cy}^{2}}-\sigma_{\rm p}\phi_{\rm cy}-b\phi_{\rm cy}^{2}$
$\displaystyle-\frac{k_{\rm B}T}{a_{\rm p}}[\phi_{\rm cy}\ln(\phi_{\rm
cy})+(1-\phi_{\rm cy})\ln(1-\phi_{\rm cy})]$ $\displaystyle=$
$\displaystyle\frac{\Pi R_{\rm sp}}{2}+\frac{f_{\rm ex}}{4\pi R_{\rm
sp}}+\frac{\kappa_{\rm p}C_{0}\phi_{\rm sp}}{R_{\rm sp}}-\sigma_{\rm
p}\phi_{\rm sp}-b\phi_{\rm sp}^{2}$ $\displaystyle-\frac{k_{\rm B}T}{a_{\rm
p}}[\phi_{\rm sp}\ln(\phi_{\rm sp})+(1-\phi_{\rm sp})\ln(1-\phi_{\rm sp})].$

The first terms in Eqs. (II.1) and (II.1) represent the Laplace tensions.
Here, the spherical and cylindrical components share the same surface tension,
which is different from the analysis in Ref. 12. Because lipid molecules can
freely move between these components, their surface tension is balanced. Note
that proteins can freely move between these components through binding and
unbinding with bulk diffusion in addition to the surface diffusion.

Figure 2:  Force $f_{\rm ex}$ dependence of (a) the protein density $\phi_{\rm
cy}$ and (b) curvature $1/R_{\rm cy}$ of the cylindrical membrane, (c) the
surface tension $\sigma_{\rm d}$, and (d) the osmotic pressure $\Pi$ for
$\mu/k_{\rm B}T=-6$, $-3$, $-2$, $-1$, and $-0.4$ at $C_{0}R_{0}=400$, $v_{\rm
r}=0.9$, $\kappa_{\rm p}/\kappa_{\rm d}=3$, and $b=0$. The solid lines
represent thermal equilibrium states. The dashed lines represent the
metastable and free-energy-barrier states. From the bottom to top, $\mu/k_{\rm
B}T=-6$, $-3$, $-2$, $-1$, and $-0.4$ in (a), (c), and (d).  Figure 3:  Free
energy profile at $C_{0}R_{0}=400$, $v_{\rm r}=0.9$, $\kappa_{\rm
p}/\kappa_{\rm d}=3$, and $b=0$. (a) $\mu/k_{\rm B}T=-0.4$. (b) $\mu/k_{\rm
B}T=-2$. The free energy has a double minima at $f_{\rm ex}/f_{\rm s}=0.02$
and $0.06$ in (a).

When the volume of the cylindrical tube is negligibly small as $R_{\rm
cy}^{2}L_{\rm cy}/R_{0}^{3}\ll 1$, the spherical component can be approximated
to have a maximum volume of $V=4\pi{R_{\rm sp}}^{3}/3$. In this limit
condition Smith  _et al._ (2004), the following relation is obtained:

$\displaystyle r$ $\displaystyle\simeq$ $\displaystyle v_{\rm r}^{1/3},$ (15)
$\displaystyle\frac{R_{\rm cy}L_{\rm cy}}{R_{0}^{2}}$ $\displaystyle\simeq$
$\displaystyle 2(1-v_{\rm r}^{2/3}).$ (16)

Smith et al. reported that this approximation reproduces the force-length
relation of tethered vesicles with zero spontaneous curvature Smith  _et al._
(2004). In this approximation, $L_{\rm cy}$ is a variable, whereas $r$ is a
constant. We call this approximation limit analysis and use sub- or
superscripts lim to represent the quantities calculated by this method. The
vesicle shape is obtained from $\partial F/\partial L_{\rm cy}=0$, as

$\displaystyle f_{\rm ex}$ $\displaystyle\simeq$
$\displaystyle\frac{\pi(\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm d})L_{\rm
cy}}{(1-v_{\rm r}^{2/3})R_{0}^{2}}-2\pi\kappa_{\rm p}C_{0}\phi_{\rm cy}$ (17)
$\displaystyle=$ $\displaystyle\frac{2\pi(\kappa_{\rm dif}\phi_{\rm
cy}+\kappa_{\rm d})}{R_{\rm cy}}-2\pi\kappa_{\rm p}C_{0}\phi_{\rm cy}$ (18)
$\displaystyle=$ $\displaystyle\Big{[}\Big{(}\frac{\kappa_{\rm
dif}}{\kappa_{\rm d}}\phi_{\rm cy}+1\Big{)}\Big{(}\frac{1}{R_{\rm cy}C_{\rm
s}}-1\Big{)}+1\Big{]}f_{\rm s},$ (19)

where

$f_{\rm s}=2\pi\kappa_{\rm d}C_{\rm s}.$ (20)

The maximum protein density is obtained at $f_{\rm ex}=f_{\rm s}$, since
$1/R_{\rm cy}=C_{\rm s}$ at the maximum. The force $f_{\rm ex}$ linearly
increases with the tube curvature $1/R_{\rm cy}$ for a constant $\phi_{\rm
cy}$. In particular, for completely unbound tubes ($\phi_{\rm cy}=0$), $f_{\rm
ex}$ is proportional to $1/R_{\rm cy}$ as $f_{\rm ex}/f_{\rm s}\simeq 1/R_{\rm
cy}C_{\rm s}$. For any value of $\phi_{\rm cy}$, $f_{\rm ex}=f_{\rm s}$ at
$R_{\rm cy}C_{\rm s}=1$. The surface tension $\sigma_{\rm d}$ is obtained from
$\partial\breve{F}/\partial R_{\rm cy}|_{L_{\rm cy}}=0$ with
$\breve{F}=F+\sigma_{\rm d}A$:

$\displaystyle\sigma_{\rm d}$ $\displaystyle\simeq$
$\displaystyle\frac{\kappa_{\rm dif}\phi_{\rm cy}+\kappa_{\rm d}}{2R_{\rm
cy}^{2}}-\sigma_{\rm p}\phi_{\rm cy}-b\phi_{\rm cy}^{2}$
$\displaystyle-\frac{k_{\rm B}T}{a_{\rm p}}[\phi_{\rm cy}\ln(\phi_{\rm
cy})+(1-\phi_{\rm cy})\ln(1-\phi_{\rm cy})].$

Hence, the first term in Eq. (II.1) is neglected in the limit analysis. For
completely unbound membrane tubes of $\phi_{\rm cy}=0$, a well-known relation
$\sigma_{\rm d}=\kappa_{\rm d}/2R_{\rm cy}^{2}=f_{\rm ex}/4\pi R_{\rm cy}$ is
obtained. This relation has been used to estimate $R_{\rm cy}$ and
$\kappa_{\rm d}$ from $\sigma_{\rm d}$ and $f_{\rm ex}$ experimentally Dimova
(2014); Bo and Waugh (1989); Evans  _et al._ (1996); Cuvelier  _et al._
(2005). However, the surface tension of the bound membranes exhibits a more
complicated dependence, as expressed in Eq. (II.1). As described in Sec. II.2,
this limit method provides a good approximation for thin tubes.

We use $\kappa_{\rm d}=20k_{\rm B}T=8\times 10^{-20}$ J, $a_{\rm p}=100$ nm2,
and $R_{0}=10\,\mu$m, i.e., $a_{\rm p}/R_{0}^{2}=10^{-6}$. External forces are
typically at $f_{\rm ex}\sim 10$ pN, which can be imposed using optical
tweezers. For $C_{\rm s}=100/R_{0}=10^{7}{\rm m}^{-1}$, $f_{\rm
s}=2\pi\kappa_{\rm d}C_{\rm s}\simeq 5$ pN. In the plots, the data at $R_{\rm
cy}/R_{0}<0.1$ are displayed (i.e., too wide tubes are excluded). In this
study, the iterations of the density calculation for $b\neq 0$ were repeated
until the difference in $\phi$ was less than $10^{-8}$; fewer than ten
iterations were typically performed.

Figure 4:  Chemical potential $\mu$ dependence of (a) the protein densities
$\phi$ and (b) the curvature $1/R_{\rm cy}$ of the cylindrical membrane for
$f_{\rm ex}/f_{\rm s}=0.25$, $1.1$, and $2$ at $C_{0}R_{0}=400$, $v_{\rm
r}=0.9$, $\kappa_{\rm p}/\kappa_{\rm d}=3$, and $b=0$. The solid lines
represent thermal equilibrium states. The dashed lines represent the
metastable and free-energy-barrier states. The gray line in (a) represents
$\phi_{\rm sp}$ (all data for $f_{\rm ex}/f_{\rm s}=0.25$, $1.1$, and $2$
overlap this single curve).  Figure 5:  Force $f_{\rm ex}$ dependence of (a)
the protein density $\phi_{\rm cy}$ and (b) curvature $1/R_{\rm cy}$ of the
cylindrical membrane for $\mu/k_{\rm B}T=-3$, $-1$, $0$, $1$, and $2$ at
$C_{0}R_{0}=100$, $v_{\rm r}=0.9$, $\kappa_{\rm p}/\kappa_{\rm d}=3$, and
$b=0$. The solid lines represent thermal equilibrium states. The dashed lines
represent the metastable and free-energy-barrier states.  Figure 6:  Chemical
potential $\mu$ dependence of (a) the protein densities $\phi$ and (b) the
curvature $1/R_{\rm cy}$ of the cylindrical membrane for $f_{\rm ex}/f_{\rm
s}=0.9$ and $5$ at $C_{0}R_{0}=100$, $v_{\rm r}=0.9$, $\kappa_{\rm
p}/\kappa_{\rm d}=3$, and $b=0$. The solid lines represent thermal equilibrium
states. The dashed lines represent the metastable and free-energy-barrier
states. The gray line in (a) represents $\phi_{\rm sp}$ (both data for $f_{\rm
ex}/f_{\rm s}=0.9$ and $5$ overlap this single curve).  Figure 7:  Phase
diagrams at $v_{\rm r}=0.9$ and $b=0$. Two states with narrow and wide tubes
coexist in the region between the two solid lines. (a) $C_{0}R_{0}=100$,
$200$, and $400$ at $\kappa_{\rm p}/\kappa_{\rm d}=3$. (b) $\kappa_{\rm
p}/\kappa_{\rm d}=2$, $3$, and $5$ at $C_{0}R_{0}=400$. The origin of the
chemical potential is at $\mu_{0}$ to give $\phi_{\rm cy}=0.5$ at $f_{\rm
ex}/f_{\rm s}=1$.  Figure 8:  Effects of reduced volume $v_{\rm r}$ at
$C_{0}R_{0}=400$, $\mu/k_{\rm B}T=-0.4$, $\kappa_{\rm p}/\kappa_{\rm d}=3$,
and $b=0$. (a) Length of the cylindrical membrane for $v_{\rm r}=0.8$, $0.9$,
and $0.95$. (b),(c) Difference of radius of the spherical membrane from the
maximum value $r_{\rm lim}$ for (b) large and (c) small $f_{\rm ex}$ regions.
The solid lines represent thermal equilibrium states. The dashed lines
represent the metastable and free-energy-barrier states.  Figure 9:
Dependence on pairwise interactions at $C_{0}R_{0}=400$, $v_{\rm r}=0.9$, and
$\kappa_{\rm p}/\kappa_{\rm d}=3$. (a) Protein densities $\phi_{\rm cy}$ and
$\phi_{\rm sp}$ of the cylindrical and spherical membrane components,
respectively, for $ba_{\rm p}/k_{\rm B}T=-1$, $0$, and $1$ at $f_{\rm
ex}/f_{\rm s}=0.15$. (b) The curvature $1/R_{\rm cy}$ of cylindrical membrane
for $ba_{\rm p}/k_{\rm B}T=-1$, $0$, and $1$ at $f_{\rm ex}/f_{\rm s}=0.15$.
(c) Phase diagrams at $\mu/k_{\rm B}T=-1.5$. Two states with narrow and wide
tubes coexist in the region between the two solid lines.  Figure 10:  Accuracy
of the limit analysis using Eqs. (15)–(II.1) for $\mu/k_{\rm B}T=-2$ and
$-0.4$ at $C_{0}R_{0}=400$, $v_{\rm r}=0.9$, $\kappa_{\rm p}/\kappa_{\rm
d}=3$, and $b=0$. (a) Difference of protein density $\Delta\phi_{\rm
cy}=\phi_{\rm cy}-\phi_{\rm cy}^{\rm lim}$. (b) Difference of the curvature of
cylindrical membrane $\Delta(1/R_{\rm cy}C_{\rm s})=1/R_{\rm cy}C_{\rm
s}-1/R_{\rm cy}^{\rm lim}C_{\rm s}$. (c) Difference of the surface tension
$\Delta\sigma_{\rm d}=\sigma_{\rm d}-\sigma_{\rm d}^{\rm lim}$. The
superscript lim represents data obtained by the limit analysis. The thermal
equilibrium and metastable states are used for $\mu/k_{\rm B}T=-0.4$.

###  II.2 Theoretical results

In this subsection, we describe the results of the mean-field theory using the
simplified geometry depicted in Fig. 1. Figures 2–4 and Figs. 5–6 show the
results at $C_{0}R_{0}=400$ and $C_{0}R_{0}=100$, respectively, at $v_{\rm
r}=0.9$, $\kappa_{\rm p}/\kappa_{\rm d}=3$, and $b=0$. When the proteins are
unbound at a low chemical potential $\mu$, the tube curvature is proportional
to the force strength, as shown in Figs. 2(b) and 5(b). As $\mu$ increases,
the protein density $\phi_{\rm cy}$ increases, and the force curve changes, as
predicted in Eqs. (10) and (19). The force dependence curves of $\phi_{\rm
cy}$ are visually reflection symmetric with respect to $f_{\rm ex}/f_{\rm
s}=1$ and take the maxima at $f_{\rm ex}/f_{\rm s}=1$ as shown in Figs. 2(a)
and 5(a). The tube curvature $1/R_{\rm cy}$ is point symmetric with respect to
$f_{\rm ex}/f_{\rm s}=1$ (see Figs. 2(b) and 5(b)). Thus, the results are
reproduced well by the limit analysis using $r=v_{\rm r}^{1/3}$. When $C_{0}$
and $\kappa_{\rm p}$ are unknown, they can be estimated from the tube radius
at the maximum value of $\phi_{\rm cy}$ and the slope of $d(1/R_{\rm
cy})/df_{\rm ex}=1/2\pi(k_{\rm dif}\phi_{\rm cy}^{\rm max}+k_{\rm d})$ at
$f_{\rm ex}/f_{\rm s}=1$, respectively.

For large $\mu$, first-order transitions occur between the unbound and bound
tubes with different tube radii (see Figs. 2 and 5). These transitions appear
symmetrically with both small and large $f_{\rm ex}$ with respect to $f_{\rm
ex}/f_{\rm s}=1$. Around the transition points, the free energy profile has
two minima (see Fig. 3), so that the unbound and bound states coexist (see van
der Waals loop depicted by the dashed lines in Figs. 2 and 5). These two
transitions can be understood using Eq. (19) of the limit analysis; at large
$\mu$, Eq. (19) has two regions of $df_{\rm ex}/d(1/R_{\rm cy})<0$, which
means an unstable solution, i.e., free-energy barrier states between stable
and metastable phases. With increasing $\mu$, the $f_{\rm ex}$ width of the
coexistence increases, and the range of $f_{\rm ex}$ is shifted outward (see
Fig. 7). At small spontaneous curvatures ($C_{0}R_{0}\lesssim 200$), the
transition at a small $f_{\rm ex}$ disappears, since it moves into the
unphysical region ($f_{\rm ex}<0$). Although the transition at $f_{\rm
ex}<f_{\rm s}$ has been previously reported in Ref. 12, the reentrant
transition at $f_{\rm ex}>f_{\rm s}$ has not yet been reported. These
transitions are similar to the budding transition between a small number of
large buds and a large number of small buds Noguchi (2021). However, one of
the specific features of the present case is that the transitions occur twice
with increasing force. This is because tubes with larger curvatures than
$C_{\rm s}$ are generated by large external forces, unlike spontaneous
budding.

The density $\phi_{\rm sp}$ of the spherical component follows the sigmoid
function $1/\\{1+\exp[-(\mu-\mu_{\rm half})/k_{\rm B}T]\\}$, as shown in Figs.
4 and 6. Since changes in the radius $R_{\rm sp}$ are very small (see Figs.
8(b) and (c)), the deviation from this sigmoid function is negligibly small.
In contrast, $\phi_{\rm cy}$–$\mu$ curves can largely deviate from the sigmoid
function accompanied by changes in $R_{\rm cy}$ (see Figs. 4 and 6).
Interestingly, the proteins can bind more onto the spherical component than
onto the cylindrical tube as shown in Fig. 6(a). This occurs in tubes that are
narrower than $1/C_{\rm s}$ at large $f_{\rm ex}$ values.

With increasing force $f_{\rm ex}$, the surface tension $\sigma_{\rm d}$ and
osmotic pressure $\Pi$ increase together. For $\phi_{\rm cy}\simeq 0$ at a low
$\mu$, $\sigma_{\rm d}=f_{\rm ex}/4\pi R_{\rm cy}$ is obtained. At larger
$\phi_{\rm cy}$ (due to large $\mu$), $\sigma_{\rm d}$ and $\Pi$ slightly
increase (see Figs. 2(c) and (d)). The vesicle ruptures when $\sigma_{\rm d}$
overcomes the lysis tension, which is typically $1$–$25$ mN/m, depending on
the membrane composition and conditions Evans and Ludwig (2000); Evans  _et
al._ (2003); Ly and Longo (2004). Since the maximum value in Fig. 2(c) is
$\sigma_{\rm d}\simeq 0.08$ mN/m, the lipid membranes are not yet ruptured in
this range. However, the tubular membranes may become unstable. The minimum
radius of the tubular membranes is $\sim 10$ nm, depending on the membrane
composition. In the right end regions in Figs. 2 and 5, the tubular membrane
can be ruptured or the tubular radius is saturated to a finite value due to
the repulsion between membranes: For $C_{0}R_{0}=400$, $R_{\rm cy}=1/2C_{\rm
s}\simeq 8$ nm and for $C_{0}R_{0}=100$, $R_{\rm cy}=1/5C_{\rm s}\simeq 13$
nm.

As the bending rigidity $\kappa_{\rm p}$ of the bound membrane increases, the
slope of $1/R_{\rm cy}$–$f_{\rm ex}$ curves at $\phi_{\rm cy}=1$ decreases,
and the first-order transitions start at a lower $\phi_{\rm cy}^{\rm max}$
(see Fig. 7(b)). As the reduced volume $v_{\rm r}$ decreases, the tube length
$L_{\rm cy}$ increases owing to the larger available area of the tube (see
Fig. 8(a)). However, the tube radius $R_{\rm cy}$ remains almost unchanged.

Next, we consider the inter-protein interactions ($b\neq 0$). As $b$
decreases, protein binding is promoted at large $\phi$ (Figs. 9(a),(b)), such
that the coexistence region widens (Fig. 9(c)). The reflection and point
symmetries of the $\phi_{\rm cy}$–$f_{\rm ex}$ and $1/R_{\rm cy}$–$f_{\rm ex}$
curves (Figs. 2(a) and (b)) remain unchanged for $b\neq 0$, respectively (data
not shown). Note that the membrane exhibits a phase separation with large and
small $\phi$ within each component at $b<-2k_{\rm B}T/a_{\rm p}$, in which
$F_{\rm sp}$ and $F_{\rm cy}$ can have double minima Noguchi (2021). Thus, we
only consider $b>-2k_{\rm B}T/a_{\rm p}$ in this study.

To confirm the quality of the limit analysis, the deviations of $\phi_{\rm
cy}$, $1/R_{\rm cy}$, and $\sigma_{\rm d}$ are shown in Fig. 10. All of these
are very small. Note that the larger values at the ends of the curves in Figs.
10(a) and (b) are due to slight shifts in the spinodal points along $f_{\rm
ex}$. Thus, this method provides an accurate approximation, so that one can
focus only on the cylindrical component to study the protein binding on the
tethered vesicle. In the next section, we simulate only the membrane tubes.

##  III Simulation of membrane tubes

###  III.1 Simulation model

A fluid membrane is represented by a self-assembled single-layer sheet of $N$
particles. The position and orientational vectors of the $i$-th particle are
${\bm{r}}_{i}$ and ${\bm{u}}_{i}$, respectively. The details of the spin
meshless membrane model are described in Ref. 27, and the combination with the
protein binding is described in Ref. 29; the model is described briefly here.

Figure 11:  Snapshots of the meshless membrane simulations at (a),(b) $N=2400$
and (c),(d) $N=9600$. (a),(b) At $f_{\rm ex}\sigma/k_{\rm B}T=43$, $\mu_{\rm
ms}/k_{\rm B}T=3.5$, and $C_{0}\sigma=0.05$, (a) a wide tube with a high
protein density and (b) narrow tube with a low protein density coexist.
(c),(d) Phase-separated membranes for (c) $C_{0}\sigma=0.05$ and (d)
$C_{0}\sigma=0.075$ at $f_{\rm ex}\sigma/k_{\rm B}T=38$ and $\mu_{\rm
ms}/k_{\rm B}T=3$. Red and transparent gray spheres represent the bound and
unbound membrane particles, respectively.  Figure 12:  Tubular membranes for
$\mu_{\rm ms}/k_{\rm B}T=2$, $2.5$, $3$, $3.5$, and $4$ at $C_{0}\sigma=0.05$.
(a) Protein density $\phi_{\rm cy}$. (b) Tube radius $R_{\rm cy}$. The symbols
with solid lines represent simulation data, and the gray lines represent the
prediction of the mean-field theory at $\mu=\mu_{\rm ms}-3.5k_{\rm B}T$ and
$ba_{\rm p}/k_{\rm B}T=-1$.  Figure 13:  Tubular membranes for $\mu_{\rm
ms}/k_{\rm B}T=2$, $2.5$, $3$, $3.5$, and $4$ at $C_{0}\sigma=0.075$. (a)
Protein density $\phi_{\rm cy}$. (b) Tube radius $R_{\rm cy}$. The symbols
with solid lines represent simulation data, and the gray lines represent the
prediction of the mean-field theory at $\mu=\mu_{\rm ms}-3.8k_{\rm B}T$ and
$ba_{\rm p}/k_{\rm B}T=-1$.

The membrane particles interact with each other via the potential $U=U_{\rm
rep}+U_{\rm att}+U_{\rm bend}+U_{\rm tilt}$. The potential $U_{\rm rep}$ is an
excluded volume interaction with diameter $\sigma$ for all pairs of particles.
The solvent is implicitly accounted for by the effective attractive potential
$U_{\rm att}$. The bending and tilt potentials are given as follows:

$\displaystyle\frac{U_{\rm bend}}{k_{\rm B}T}$ $\displaystyle=$
$\displaystyle\frac{k_{\rm
bend}}{2}\sum_{i<j}({\bm{u}}_{i}-{\bm{u}}_{j}-C_{\rm
bd}\hat{\bm{r}}_{i,j})^{2}w_{\rm cv}(r_{i,j}),$ (22)
$\displaystyle\frac{U_{\rm tilt}}{k_{\rm B}T}$ $\displaystyle=$
$\displaystyle\frac{k_{\rm
tilt}}{2}\sum_{i<j}[({\bm{u}}_{i}\cdot\hat{\bm{r}}_{i,j})^{2}+({\bm{u}}_{j}\cdot\hat{\bm{r}}_{i,j})^{2}]w_{\rm
cv}(r_{i,j}),\hskip 19.91684pt$ (23)

respectively, where $\hat{\bm{r}}_{i,j}={\bm{r}}_{i,j}/r_{i,j}$,
${\bm{r}}_{i,j}={\bm{r}}_{i}-{\bm{r}}_{j}$, and $w_{\rm cv}(r_{i,j})$ is a
weight function. The spontaneous curvature is given by $C_{0}=C_{\rm
bd}/2\sigma$ Shiba and Noguchi (2011).

Each membrane particle is a binding site and takes two states (bound and
unbound). In this study, $C_{0}=0$ and $k_{\rm bend}=k_{\rm tilt}=10$ for the
unbound membrane particles and $k_{\rm bend}=k_{\rm tilt}=40$ for the bound
membrane particles, where $\kappa_{\rm d}/k_{\rm B}T=16\pm 1$ and $\kappa_{\rm
p}/k_{\rm B}T=71\pm 3$. In the bending and tilt potentials, for a pair of
neighboring bound and unbound particles, we use the mean value $k_{\rm
bend}=k_{\rm tilt}=25$. For the bound membrane, $C_{0}\sigma=0.05$ and $0.075$
are used. The ratio of the Gaussian modulus $\bar{\kappa}$ to $\kappa$ is
constant as follows: $\bar{\kappa}/\kappa=-0.9\pm 0.1$ Noguchi (2019b). The
other parameters are the same as those used in Ref. 29. For the unbound and
bound particles, the membrane areas per particle are $1.251\sigma^{2}$ and
$1.262\sigma^{2}$, respectively, for tensionless membranes at $C_{0}=0$.

A tubular membrane consisting of $N$ particles is set along the $x$-axis
connected by the periodic boundary condition. The force $f_{\rm ex}$ is
imposed along the $x$-axis, such that the tube length fluctuates thermally.
Since the solvent is not explicitly taken into account, the volume of the
membrane tube can freely change, as assumed in the limit analysis. We mainly
used $N=2400$. In addition, four-fold longer tubes with $N=9600$ were used at
several parameter sets, to examine finite-size effects. For $f_{\rm
ex}\sigma/k_{\rm B}T\leq 2$, $N=19200$ was also used. Membrane motion is
solved by molecular dynamics with a Langevin thermostat Allen and Tildesley
(1987); Feller  _et al._ (1995); Noguchi (2011, 2012). The bound and unbound
states are stochastically switched by a Metropolis Monte Carlo procedure with
$\Delta H=\Delta U-\mu_{\rm ms}$, where $\Delta U$ is the energy difference
between the bound and unbound states and $\mu_{\rm ms}$ is the binding
chemical potential of the membrane particles Goutaland  _et al._ (2021). Error
bars are estimated from three independent runs.

###  III.2 Simulation results

Simulation results are shown in Figs. 11–13. As $f_{\rm ex}$ increases, the
protein density $\phi_{\rm cy}$ and tube radius $R_{\rm cy}$ decrease. For
$\mu_{\rm ms}/k_{\rm B}T=3$ and $3.5$ at $C_{0}\sigma=0.05$, two states (large
$R_{\rm cy}$ with high $\phi_{\rm cy}$ and small $R_{\rm cy}$ with low
$\phi_{\rm cy}$) coexist (see Figs. 11(a), (b), and 12). These results show
very good agreement with the results of the mean-field theory with the
following parameter sets (see Figs. 12 and 13). Since the binding/unbinding
processes in the simulation involve the energy change in the other potentials
($U_{\rm rep}$ and $U_{\rm att}$) via a small area change, the chemical
potential is shifted from the input value $\mu_{\rm ms}$ Goutaland  _et al._
(2021). Thus, we use $\mu=\mu_{\rm ms}-3.5k_{\rm B}T$ and $\mu_{\rm
ms}-3.8k_{\rm B}T$ with $a_{\rm p}=1.2\sigma^{2}$ for $C_{0}\sigma=0.05$ and
$0.075$, respectively. In the simulation, no direct interactions are
considered between the bound sites. However, the Casimir-like attractive
interactions occur between them owing to the bending rigidity difference
Goutaland  _et al._ (2021). To mimic this, $ba_{\rm p}/k_{\rm B}T=-1$ is used
here. In the theory, $\phi_{\rm cy}$ has a maximum at $f_{\rm ex}=f_{\rm s}$
($f_{\rm s}\sigma/k_{\rm B}T=6.5$ and $9.7$ for $C_{0}\sigma=0.05$ and
$0.075$, respectively). Indeed, $\phi_{\rm cy}^{\rm max}$ is obtained at
$f_{\rm ex}\simeq f_{\rm s}$ in the simulation for $\phi_{\rm cy}^{\rm
max}\lesssim 0.5$ at $\mu_{\rm ms}=2k_{\rm B}T$; however, it slightly
increases with increasing $\mu_{\rm ms}$ to $f_{\rm ex}\simeq 2f_{\rm s}$ at
$\mu_{\rm ms}=4k_{\rm B}T$. Because the Casimir-like forces are not pairwise,
their multibody interactions likely induce this dependency. At a small force
$f_{\rm ex}$ and low $\phi_{\rm cy}$, the radius $R_{\rm cy}$ calculated by
the simulation is slightly larger than the theoretical prediction. This is due
to the thermal undulation of the membrane. Except for these, the simulation
results are well reproduced.

Homogeneous membrane tubes are not maintained in some ranges of the simulation
parameters. First, the membrane tubes become unstable at large and small
limits of $f_{\rm ex}$. The membranes are ruptured for radii that are too
small $R_{\rm cy}\lesssim 2\sigma$ at large $f_{\rm ex}$. At $f_{\rm
ex}\sigma/k_{\rm B}T\lesssim 1$, the cylindrical membrane becomes unstable and
divided into vesicles (see Movie 1 provided in ESI for $f_{\rm ex}=0$,
$\mu_{\rm ms}=4k_{\rm B}T$, $C_{0}\sigma=0.075$, and $N=19200$). At $f_{\rm
ex}=0$, a membrane tube can form an unduloid shape Kenmotsu (2003); Naito  _et
al._ (1995), in which a constant value of the mean curvature $H$ is maintained
everywhere. The cylindrical tube at $f_{\rm ex}=0$ has the curvature of the
curvature generation, as $1/R_{\rm cy}=2H_{\rm g}=\kappa_{\rm
p}\langle\phi_{\rm cy}\rangle C_{0}/(\kappa_{\rm dif}\langle\phi_{\rm
cy}\rangle+\kappa_{\rm d})$ from Eq. (18). However, longer tubes than the
unduloid wavelength $l_{\rm und}=2\pi R_{\rm cy}$ (large $\mu_{\rm ms}$ at
$N=9600$ and $19200$) deform into unduloid shapes and subsequent membrane
fission at the narrow neck of the tube leads to the formation of spherical
vesicles. Cylindrical shapes maintain for shorter tubes than the unduloid
wavelength owing to the finite-size effect.

In addition, phase separation can unstabilize cylindrical shapes. A beaded-
necklace-like tube with phase separation is formed in the middle-density
region ($0.2\lesssim\phi_{\rm cy}\lesssim 0.5$) around the critical points
(see Figs. 11(c) and (d)). The bound membranes form an ellipsoidal shape,
whereas the unbound membranes form a saddle shape between them. In Figs. 12
and 13, the data of the phase-separated tubes are not plotted. The homogeneous
phases disappear for $3.2\lesssim f_{\rm ex}\sigma/k_{\rm B}T\lesssim 3.9$ at
$C_{0}\sigma=0.075$ and $\mu_{\rm ms}/k_{\rm B}T=3$ (see Figs. 11(c) and Movie
2 provided in ESI using the same parameter set with different initial
conformations). This beaded-necklace-like tube is also formed when the
homogeneous tube becomes unstable at the ends of the metastable state of the
coexistence region (see Figs. 11(d)). Membrane is ruptured when a neck region
becomes too narrow. The coexistence region is narrower in the simulation than
in the theory. The phase separation likely causes this reduction of the
coexistence region.

Similar necklace-like membrane tubes with phase separation have been
previously observed in the experiments on three-component membranes Yanagisawa
_et al._ (2008, 2010). Moreover, membrane fission has been observed in the
necked region Allain  _et al._ (2004). In their system, the difference in
bending rigidity causes a neck-like shape and the phase separation occurs even
in a flat membrane. In contrast, for the present system, the difference in the
spontaneous curvature gives more dominant effects.

Unduloid-like deformations have been reported in tethered vesicles after fore
release Tozzi  _et al._ (2019), as well as in tubular vesicles with polymer
anchoring Tsafrir  _et al._ (2001) and rolled membranes during detachment from
a substrate Noguchi (2019c). In contrast, phase separation can make
cylindrical tubes unstable under strong forces. This beaded-necklace-like
shape may be expressed by a periodic combination of two constant-curvature
surfaces discussed in Ref. 57 with an extension to include the deformation by
the external force.

Although Fig. 11 shows only necklace-like membranes of longer tubes, the
shorter tubes with $N=2400$ also show this structure with fewer periodicities
(one or two circular bumps). The simulation results exhibit no notable
differences between these two tube sizes, except for slightly wider
coexistence regions for the shorter tubes.

Here, we used a constant chemical potential and constant external force. The
same membrane shapes can be obtained with the ensemble of a constant number
$N_{\rm p}$ of the proteins and/or constant tube length $L_{\rm cy}$, when the
condition is adjusted. Note that for fixed $N_{\rm p}$ or $L_{\rm cy}$, phase
separation appears more often in the regions of the first-order transitions
owing to the macroscopic phase separation. For example, at $\phi_{\rm
cy}\simeq 0.4$ and $f_{\rm ex}\sigma/k_{\rm B}T>40$, no solution exists in
Figs. 12 and 13. Thus, when $N_{\rm p}/N=0.4$, the membrane is separated into
two regions of high and low protein densities. This is similar to the
gas–liquid coexistence in the $NVT$ ensemble Watanabe  _et al._ (2012).

##  IV Summary and discussions

We have studied the binding of curvature-inducing proteins onto the tethered
vesicle. Proteins exhibit an isotropic spontaneous curvature such that they
sense and generate the curvature of membranes. For a completely unbound
membrane, the tube curvature and tube length are proportional to the force
strength. As the binding chemical potential increases, the protein density
sigmoidally increases in the spherical component of the vesicle. In contrast,
a discrete increase can also occur in the membrane tube, accompanied by a
change in the tube radius. The force–density curve and force–tube-curvature
curve are reflection and point symmetric to the point, where the tube
curvature equals to the sensing curvature, respectively. The approximation
that neglects the tube volume for a small tube radius well reproduces the
results of the analysis for a finite volume. Meshless simulations of membrane
tubes were conducted to confirm these theoretical results. The results of the
simulation and theory for the homogeneous phases show very good agreement.
Additionally, in the simulation, beaded-necklace-like membrane tubes with
phase separation are found around the critical points. Membrane deformation
induces this microphase separation.

Based on these mean-field analyses, we propose a method for estimating the
bending rigidity change by protein binding. The sensing curvature is obtained
as the tube curvature at the maximum protein density $\phi_{\rm cy}^{\rm max}$
under variation in the tube radius. The bending rigidity ratio is obtained
from the slope of $d(1/R_{\rm cy})/df_{\rm ex}$ at this maximum protein
density as expressed in Eq. (19). Our simulation suggest that not large
density, $\phi_{\rm cy}^{\rm max}\lesssim 0.5$, should be used for these
estimations to avoid the influence of inter-protein interactions. The maximum
density can be varied by the binding chemical potential $\mu$, which is a
function of the buffer protein concentration $\rho$. For a dilute solution,
$\mu(\rho)=\mu(1)+k_{\rm B}T\ln(\rho)$.

In this study, we consider only proteins that isotropically bend the membrane
and have no preferred lateral direction. The BAR superfamily proteins exhibit
anisotropic spontaneous curvatures. Previously, we have reported that such
proteins show characteristic behavior in tubular membranes Noguchi (2016a,
2014, 2015, 2019a); the force–protein-curvature curve has a flat region at low
curvature owing to the adjustment of the protein orientation, and protein
assembly induces elliptic and polyhedral tube formations. However, the mean-
field theory of the isotropic spontaneous curvature has been used to analyze
the experimental results on the binding of the BAR proteins to the tethered
membranes Prévost  _et al._ (2015). Recently, a mean-field theory for a
nematic order coupled with protein bending energy was developed for a fixed-
shaped membrane Tozzi  _et al._ (2021) based on Nascimentos’ liquid-crystal
theory Nascimento  _et al._ (2017). The present method can be extended to
anisotropic proteins by including the orientational degree. This is one of the
directions for further studies.

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Number JP21K03481.

## References

  * McMahon and Gallop (2005) H. T. McMahon and J. L. Gallop, Nature 438, 590 (2005).
  * Suetsugu  _et al._ (2014) S. Suetsugu, S. Kurisu, and T. Takenawa, Physiol. Rev. 94, 1219 (2014).
  * Johannes  _et al._ (2015) L. Johannes, R. G. Parton, P. Bassereau, and S. Mayor, Nat. Rev. Mol. Cell. Biol. 16, 311 (2015).
  * Brandizzi and Barlowe (2013) F. Brandizzi and C. Barlowe, Nat. Rev. Mol. Cell Biol. 14, 382 (2013).
  * Hurley  _et al._ (2010) J. H. Hurley, E. Boura, L.-A. Carlson, and B. Różycki, Cell 143, 875 (2010).
  * McMahon and Boucrot (2011) H. T. McMahon and E. Boucrot, Nat. Rev. Mol. Cell. Biol. 12, 517 (2011).
  * Itoh and De Camilli (2006) T. Itoh and P. De Camilli, Biochim. Biophys. Acta 1761, 897 (2006).
  * Mim and Unger (2012) C. Mim and V. M. Unger, Trends Biochem. Sci. 37, 526 (2012).
  * Dimova (2014) R. Dimova, Adv. Colloid Interface Sci. 208, 225 (2014).
  * Baumgart  _et al._ (2011) T. Baumgart, B. R. Capraro, C. Zhu, and S. L. Das, Annu. Rev. Phys. Chem. 62, 483 (2011).
  * Sorre  _et al._ (2012) B. Sorre, A. Callan-Jones, J. Manzi, B. Goud, J. Prost, P. Bassereau, and A. Roux, Proc. Natl. Acad. Sci. USA 109, 173 (2012).
  * Prévost  _et al._ (2015) C. Prévost, H. Zhao, J. Manzi, E. Lemichez, P. Lappalainen, A. Callan-Jones, and P. Bassereau, Nat. Commun. 6, 8529 (2015).
  * Roux  _et al._ (2010) A. Roux, G. Koster, M. Lenz, B. Sorre, J.-B. Manneville, P. Nassoy, and P. Bassereau, Proc. Natl. Acad. Sci. USA 107, 4141 (2010).
  * Rosholm  _et al._ (2017) K. R. Rosholm, N. Leijnse, A. Mantsiou, V. Tkach, S. L. Pedersen, V. F. Wirth, L. B. Oddershede, K. J. Jensen, K. L. Martinez, N. S. Hatzakis, P. M. Bendix, A. Callan-Jones, and D. Stamou, Nat. Chem. Biol. 13, 724 (2017).
  * Allard  _et al._ (2020) A. Allard, M. Bouzid, T. Betz, C. Simon, M. Abou-Ghali, J. Lemière, F. Valentino, J. Manzi, F. Brochard-Wyart, K. Guevorkian, J. Plastino, M. Lenz, C. Campillo, and C. Sykes, Sci. Adv. 6, eaaz3050 (2020).
  * Hotani  _et al._ (1999) H. Hotani, F. Nomura, and Y. Suzuki, Curr. Opin. Coll. Interface Sci. 4, 358 (1999).
  * Inaba  _et al._ (2005) T. Inaba, A. Ishijima, M. Honda, F. Nomura, K. Takiguchi, and H. Hotani, J. Mol. Biol. 348, 325 (2005).
  * Wu  _et al._ (2016) Z. Wu, H. Yuan, X. Zhang, and X. Yi, Soft Matter 11, 2552 (2016).
  * Gavriljuk  _et al._ (2021) K. Gavriljuk, B. Scocozza, F. Ghasemalizadeh, H. Seidel, A. P. Nandan, M. Campos-Medina, M. Schmick, A. Koseska, and P. I. H. Bastiaens, Nat. Commun. 12, 1548 (2021).
  * Svitkina (2018) T. Svitkina, Cold Spring Harb. Perspect. Biol. 10, a018267 (2018).
  * Gallop (2020) J. L. Gallop, Semin. Cell Dev. Biol. 102, 81 (2020).
  * Noguchi (2021) H. Noguchi, Phys. Rev. E 104, 014410 (2021).
  * Müller  _et al._ (2006) M. Müller, K. Katsov, and M. Schick, Phys. Rep. 434, 113 (2006).
  * Venturoli  _et al._ (2006) M. Venturoli, M. M. Sperotto, M. Kranenburg, and B. Smit, Phys. Rep. 437, 1 (2006).
  * Noguchi (2009) H. Noguchi, J. Phys. Soc. Jpn. 78, 041007 (2009).
  * Noguchi and Gompper (2006) H. Noguchi and G. Gompper, Phys. Rev. E 73, 021903 (2006).
  * Shiba and Noguchi (2011) H. Shiba and H. Noguchi, Phys. Rev. E 84, 031926 (2011).
  * Noguchi (2016a) H. Noguchi, Phys. Rev. E 93, 052404 (2016a).
  * Goutaland  _et al._ (2021) Q. Goutaland, F. van Wijland, J.-B. Fournier, and H. Noguchi, Soft Matter 17, 5560 (2021).
  * Noguchi (2014) H. Noguchi, EPL 108, 48001 (2014).
  * Noguchi (2015) H. Noguchi, J. Chem. Phys. 143, 243109 (2015).
  * Noguchi (2016b) H. Noguchi, Sci. Rep. 6, 20935 (2016b).
  * Noguchi and Fournier (2017) H. Noguchi and J.-B. Fournier, Soft Matter 13, 4099 (2017).
  * Noguchi (2017) H. Noguchi, Soft Matter 13, 7771 (2017).
  * Noguchi (2019a) H. Noguchi, Sci. Rep. 9, 11721 (2019a).
  * Noguchi (2019b) H. Noguchi, J. Chem. Phys. 151, 094903 (2019b).
  * Noguchi (2019c) H. Noguchi, Soft Matter 15, 8741 (2019c).
  * Smith  _et al._ (2004) A.-S. Smith, E. Sackmann, and U. Seifert, Phys. Rev. Lett. 92, 208101 (2004).
  * Powers  _et al._ (2002) T. R. Powers, G. Huber, and R. E. Goldstein, Phys. Rev. E 65, 041901 (2002).
  * Bo and Waugh (1989) L. Bo and R. Waugh, Biophys. J. 55, 509 (1989).
  * Evans  _et al._ (1996) E. Evans, H. Bowman, A. Leung, D. Needham, and D. Tirrell, Science 273, 933 (1996).
  * Cuvelier  _et al._ (2005) D. Cuvelier, I. Derényi, P. Bassereau, and P. Nassoy, Biophys. J. 88, 2714 (2005).
  * Evans and Ludwig (2000) E. A. Evans and F. Ludwig, J. Phys. Condens. Matter 12, A315 (2000).
  * Evans  _et al._ (2003) E. A. Evans, V. Heinrich, F. Ludwig, and W. Rawicz, Biophys. J. 85, 2342 (2003).
  * Ly and Longo (2004) H. V. Ly and M. L. Longo, Biophys. J. 87, 1013 (2004).
  * Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_ (Clarendon Press, Oxford, 1987).
  * Feller  _et al._ (1995) S. E. Feller, Y. Zhang, R. W. Pastor, and B. R. Brooks, J. Chem. Phys. 103, 4613 (1995).
  * Noguchi (2011) H. Noguchi, J. Chem. Phys. 134, 055101 (2011).
  * Noguchi (2012) H. Noguchi, Soft Matter 8, 3146 (2012).
  * Kenmotsu (2003) K. Kenmotsu,  _Surfaces with constant mean curvature_ (American Mathematical Society, Providence, R.I., 2003).
  * Naito  _et al._ (1995) H. Naito, M. Okuda, and O.-Y. Zhong-can, Phys. Rev. Lett. 74, 4345 (1995).
  * Yanagisawa  _et al._ (2008) M. Yanagisawa, M. Imai, and T. Taniguchi, Phys. Rev. Lett. 100, 148102 (2008).
  * Yanagisawa  _et al._ (2010) M. Yanagisawa, M. Imai, and T. Taniguchi, Phys. Rev. E 82, 051928 (2010).
  * Allain  _et al._ (2004) J.-M. Allain, C. Storm, A. Roux, M. B. Amar, and J.-F. Joanny, Phys. Rev. Lett. 93, 158104 (2004).
  * Tozzi  _et al._ (2019) C. Tozzi, N. Walani, and M. Arroyo, New J. Phys. 21, 093004 (2019).
  * Tsafrir  _et al._ (2001) I. Tsafrir, D. Sagi, T. Arzi, M.-A. Guedeau-Boudeville, V. Frette, D. Kandel, and J. Stavans, Phys. Rev. Lett. 86, 1138 (2001).
  * Góźdź and Gompper (1999) W. T. Góźdź and G. Gompper, Phys. Rev. E 59, 4305 (1999).
  * Watanabe  _et al._ (2012) H. Watanabe, N. Ito, and C.-K. Hu, J. Chem. Phys. 136, 204102 (2012).
  * Tozzi  _et al._ (2021) C. Tozzi, N. Walani, A.-L. L. Roux, P. Roca-Cusachs, and M. Arroyo, Soft Matter 17, 3367 (2021).
  * Nascimento  _et al._ (2017) E. S. Nascimento, P. Palffy-Muhoray, J. M. Taylor, E. G. Virga, and X. Zheng, Phys. Rev. E 96, 022704 (2017).