# Anomalous magnetic flux via junction twist-angle in a triplet-
superconducting transmon qubit

Sebastián Domínguez-Calderón  School of Engineering and Sciences, Tecnológico
de Monterrey, Monterrey 64849, México  Harley Scammell  School of Mathematical
and Physical Sciences, University of Technology Sydney, Ultimo, NSW 2007,
Australia

###### Abstract

Superconducting transmon qubits with strong anharmonicity and insensitivity to
offset charge are highly desirable for low-error implementation. In this work
we propose a c-axis junction, comprising triplet superconductors, and set at a
relative twist angle. Invoking spin-orbit coupling and spin polarization,
which are known to occur in the material platform of choice, we examine the
resulting transmon Hamiltonian. This junction allows for direct control of the
single and double Cooper pair tunneling strength, and most remarkably, an
anomalous magnetic flux—i.e. a phase offset equivalent to magnetic flux, yet
in zero magnetic field. Having control over these three parameters—single and
double pair tunneling and anomalous flux—allows for optimal design of the
transmon qubit. Interestingly, in this architecture, the anomalous flux is
determined by the twist angle of the junction, thereby offering a novel zero-
field functionality. Our key results rely on symmetry arguments, for
concreteness we demonstrate the implementation of our concept using a model of
moiré graphene-based c-axis junctions.

first keyword, second keyword, third keyword

##  I Introduction

Superconducting quantum computers have undergone significant improvements and
have been utilized to showcase instances of dramatic speedup compared to their
classical counterparts Arute and et al (2019). However, a critical bottleneck
in their advancement is errors Kjaergaard  _et al._ (2020); Gyenis  _et al._
(2021). As the number of qubits increases, so does the number of error
sources, necessitating these errors lie beneath a certain threshold;
otherwise, adding more qubits becomes futile. Since true fault-tolerant
quantum computation is believed to require formidable hardware resources
Fowler  _et al._ (2012); Reiher  _et al._ (2017), the near-term goal is to dip
below the error threshold. At present, superconducting quantum computers are
approaching the error threshold Acharya and et al (2023), underscoring the
crucial importance of further reducing errors in superconducting qubits as a
primary avenue of research.

Recent advancements have pinpointed favorable transmon qubit Hamiltonians that
offer tunable anharmonicity and offset charge insensitivity, thereby
mitigating errors. These Hamiltonians consist of strong double Cooper pair
(CP) tunneling, weak single CP tunneling, and a magnetic flux Patel  _et al._
(2024); Staples  _et al._ (2024). One such architecture is a triple junction
d-mon subject to magnetic flux in one junction, whereby the d-mon refers to a
d-wave/s-wave c-axis superconducting junction. The second is a non-ideal
$\pi$-SQUID, as discussed in various works Kitaev (2006); Paolo  _et al._
(2019); Ciaccia  _et al._ (2024); Valentini  _et al._ (2024), but further
subjected to a magnetic flux Staples  _et al._ (2024).

In contrast to these works, the present study is motivated by two key
questions: (i) Since magnetic fields are known inhibitors to achieving small
scale devices, can we achieve the desired flux without a magnetic field? (ii)
All things being equal, a single-junction qubit is more suited to small scales
than a multi-junction qubit; can we design a highly tunable transmon qubit
within a single junction?

The present work succeeds in designing a transmon architecture with a high
tunability, within a single junction and without magnetic field. Our effective
transmon maps onto two of the aforementioned recent architectures: the triple
junction d-mon with magnetic flux Patel  _et al._ (2024), and the (two
junction) non-ideal $\pi$-SQUID also subject to a magnetic flux Staples  _et
al._ (2024). Again we stress that our construction achieves this desirable
qubit architecture using all-electrical control (zero magnetic field) and
within a single junction. Our key innovation is the use of triplet
superconductivity in the presence of spin-orbit coupling (SOC) and intrinsic
magnetization. We consider a junction of spin-triplet superconductors with
vectors at some non-zero angle $\Phi_{J}$ relative to each other on either
side of the junction. Then, due to SOC, the twist angle of the junction
determines $\Phi_{J}$. Finally, this $\Phi_{J}$ acts, via a non-linear
relationship, as an anomalous magnetic flux. This is our key finding.

Our results are of a general nature, however, for several reasons, we choose
to model the key features of moiré graphene, e.g. twisted bilayer graphene,
because: (i) we require triplet superconductivity, but without nodes in the
order parameter, since nodes introduce decoherence Fominov  _et al._ (2003). A
nice way to achieve this is to include an extra quantum number, in our case we
choose a valley system which admits a valley quantum number. (ii) We require a
valley system that can host superconductivity and moreover, triplet
superconductivity. (iii) We require a superconducting valley system that can
simultaneously host magnetic order. (iv) We require a system into which SOC
can be induced. These criteria naturally direct us to moiré graphene. Notably,
for (iv) one can proximity-induce SOC into moiré graphene via a transition
metal dichalcogenide (TMD) substrate. The form of the induced SOC can be
readily tuned based on the relative twist angle graphene and TMD layers Naimer
_et al._ (2021); Li and Koshino (2019); David  _et al._ (2019); Lee  _et al._
(2022); Péterfalvi  _et al._ (2022); Zollner and Fabian (2021); Zollner  _et
al._ (2022).

##  II Background: Triplet superconductivity in a valley system

We provide a short summary of key details of triplet superconductivity within
a valley system, subject to various particle-hole orders (e.g. spin
polarization). The particular valley system considered here is modeled on
moiré graphene.

First, we consider the low-energy moiré bands, and to be explicit, we consider
doping with chemical potential $\mu>0$. The partially filled moiré band is
denoted $\varepsilon_{\bm{k}}$ and is four-fold degenerate due to spin and
monolayer-valley degrees of freedom; Pauli matrices acting on these quantum
numbers are denoted by $s_{\nu}$ and $\tau_{\nu}$, respectively, with
$\nu=0,x,y,z$. Moreover, we account for trigonal warping within the
dispersion,

$\displaystyle\varepsilon_{\bm{k}}=\varepsilon^{0}_{\bm{k}}(1+\lambda\tau_{z}\cos
3\theta_{\bm{k}})$ (1)

where $\varepsilon^{0}_{\bm{k}}$ is rotationally symmetric, while $\cos
3\theta_{\bm{k}}$ reduces this symmetry to three-fold, where $\theta_{\bm{k}}$
is the quasi-momentum angle about a given valley. Here $\lambda<1$ is a
dimensionless parameter. See Fig. 1(a) for corresponding Fermi surfaces.

Next, our model assumes the Cooper pairs are formed from the partially-filled
bands, consistent with experimental observations Park  _et al._ (2021); Hao
_et al._ (2021); Lin  _et al._ (2022); Siriviboon  _et al._ (2021). Moreover,
we take the Cooper pairs to be extended-s-wave, i.e. being fully gapped and
transforming trivially about the Fermi surface of a given valley. These orders
have been discussed at various points in the literature Xu and Balents (2018);
You and Vishwanath (2019); Scheurer and Samajdar (2020); Lake  _et al._
(2022); Scammell  _et al._ (2022a); Scammell and Scheurer (2023, 2024).
Explicitly, we consider spin-triplet, valley-singlet ordering, with
corresponding gap function

$\displaystyle{\Delta}_{\bm{k}}$
$\displaystyle=\Upsilon_{\bm{k}}({\bm{d}}\cdot\bm{{s}}){s}_{y}{\tau}_{y}.$ (2)

Here
$\Upsilon_{\bm{k}}=\Upsilon_{\bm{k}}^{+}\delta_{\tau,+}+\Upsilon_{\bm{k}}^{-}\delta_{\tau,-}$
is composed of two complex scalar functions $\Upsilon_{\bm{k}}^{\tau}$
containing the momentum dependence of the superconducting order about the
Fermi surfaces (at each valley $\tau$); see the schematic of Fig. 1(a).
Meanwhile, the three-component vector ${\bm{d}}$ accounts for the spin triplet
order parameter.

Finally, in addition to superconductivity, we include in our modeling the
occurrence of different types of particle-hole orders, such as spin or valley
polarization. We now provide a toy-model Hamiltonian for the partially filled
moiré bands of moiré graphene with account of particle-hole orders

$\displaystyle
h_{\bm{k}}=\varepsilon_{\bm{k}}-\mu+M_{\mu\nu}\tau_{\mu}s_{\nu}$ (3)

here $M_{\mu\nu}$ account for different possible particle-hole order, e.g.
$M_{0z}$ is an out-of-plane spin polarization. A schematic of these key
ingredients is provided in Fig. 1(a).

Figure 1: Model details: (a)i&ii Fermi surfaces of
$\varepsilon^{0}_{\tau\bm{k}}$ for $\tau=\pm$; gap function
$\Upsilon_{\bm{k}}^{\tau}$ is superimposed as the shaded region. (a)iii&iv
Fermi surfaces with account of spin-valley locked polarization
$\varepsilon_{\tau\bm{k}}=\varepsilon_{\tau\bm{k}}+M_{zz}\tau_{z}s_{z}$, blue
and red correspond to the two spin species. (b)i Schematic of the junction and
(b)ii top view of the twisted junction, with twist angle $\Phi_{J}$. (b)iii
Shows the corresponding circuit. (c) Qubit subspace of the transmon
spectrum—red eigenvalues $\\{\ket{0}^{\prime},\ket{1}^{\prime}\\}$ corresponds
to zero twist angle $\Phi_{J}^{\prime}=0$ (denoted with a prime for clarity),
blue eigenvalues $\\{\ket{0},\ket{1}\\}$ corresponds to a non-zero twist
angle, $\Phi_{J}\neq 0$.

##  III Methods: effective transmon Hamiltonian

We model a setup of two stacked moiré graphene subsystems, separated by a thin
insulating layer, and rotated to an angle $\Phi_{J}$. An illustration is
provided in Fig. 1(b). In this setup, we treat the systems as being two copies
of moiré graphene, and crucially, the superconducting order from one system is
understood to influence the other via the proximity effect. We refer to this
as the moiré graphene junction — the relation to a standard Josephson junction
will become apparent shortly.

In the Bogoliubov-de Gennes (BdG) framework, i.e. a mean-field Hamiltonian
that captures the electronic normal state as well as the superconducting order
parameters, the moiré graphene junction is modeled as,

${H}_{\bm{k}}=\frac{1}{2}\begin{bmatrix}{h}_{\bm{k}}&{\Delta}^{(1)}_{\bm{k}}+e^{i\varphi}{\Delta}^{(2)}_{\bm{k}}\\\
{\Delta}^{\dagger(1)}_{\bm{k}}+e^{-i\varphi}{\Delta}^{\dagger(2)}_{\bm{k}}&-{h}_{-\bm{k}}\\\
\end{bmatrix},$ (4)

where ${\Delta}^{(1)}_{\bm{k}}$ and ${\Delta}^{(2)}_{\bm{k}}$ are the gap
functions of subsystem 1 and 2, which take the form of (2). Note that the
magnitude of ${\Delta}^{(2)}_{\bm{k}}$ will be reduced since it is proximity
induced into subsystem 1, see e.g. Patel  _et al._ (2024) for an analogous
calculation. We have made explicit the U(1) phase difference between these
order parameters. We stress now that this phase difference, $\varphi$, and its
canonically conjugate angular momentum, ${n}=-i\partial_{\varphi}$, will be
the dynamical phase-space variables of the effective moiré graphene transmon
Hamiltonian.

To arrive at the transmon Hamiltonian, we compute the free energy to quartic
order in the parameter ${\bm{d}}$, of (2). Such an expansion captures both
single and double tunneling processes. Explicitly, the (regularized) free
energy and subsequent expansion are given by

$\displaystyle{\cal U}$
$\displaystyle=\text{Tr}\left[\ln(\omega-{{H}}_{\bm{k}})-\ln(\omega-{{H}}^{(0)}_{\bm{k}})\right]$
$\displaystyle\simeq-\frac{1}{2}\text{Tr}\left[{\cal O}_{\bm{k}}{{\cal
G}}_{\omega,\bm{k}}{\cal O}_{\bm{k}}{{\cal
G}}_{\omega,\bm{k}}\right]-\frac{1}{4}\text{Tr}\left[\left({\cal
O}_{\bm{k}}{{\cal G}}_{\omega,\bm{k}}{\cal O}_{\bm{k}}{{\cal
G}}_{\omega,\bm{k}}\right)^{2}\right]$ $\displaystyle\simeq
a_{1}\cos\varphi+b_{1}\sin\varphi+a_{2}\cos 2\varphi+b_{2}\sin 2\varphi$ (5)

here $H^{(0)}_{\bm{k}}\equiv H_{\bm{k}}\big{|}_{\Delta^{(i)}_{\bm{k}}\to 0}$;
${{\cal
G}}_{\omega,\bm{k}}={G}_{\omega,\bm{k}}(\eta_{0}+\eta_{z})/2-{G}_{-\omega,-\bm{k}}(\eta_{0}-\eta_{z})/2$
and ${G}_{\omega,\bm{k}}=[i\omega-{h}_{\bm{k}}]^{-1}$; and ${{\cal
O}}_{\bm{k}}=\eta_{+}({\Delta}^{(1)}_{\bm{k}}+e^{i\varphi}{\Delta}^{(2)}_{\bm{k}})+\eta_{-}({\Delta}^{\dagger(1)}_{\bm{k}}+e^{-i\varphi}{\Delta}^{\dagger(2)}_{\bm{k}})$,
with $\eta_{\mu}$ Pauli-matrices accounting for the Nambu space of the BdG
Hamiltonian (4). The approximations in passing from the first-to-second and
second-to-third lines involves ignoring higher-order tunneling processes, i.e.
$\cos n\varphi$ for $n>2$, and ignoring an irrelevant $\varphi$-independent
offset.

Meanwhile, the charging energy of the junction takes the standard form,

$\displaystyle{\cal K}=4E_{C}({n}-n_{g})^{2},$ (6)

with ${n}=-i\partial_{\varphi}$ and $n_{g}$ an offset gate charge. Together,
the moiré graphene transmon Hamiltonian is ${\cal H}={\cal K}+{\cal U}$.
Truncating at double CP tunnelling, working in units of $4E_{C}=1$, and
choosing a simpler gauge than (III),

${\cal H}=({n}-n_{g})^{2}+\alpha_{1}\cos(\varphi+\Phi)+\alpha_{2}\cos
2\varphi.$ (7)

The parameters $\\{\alpha_{1},\alpha_{2},\Phi\\}$ control the anharmonicity
and charge offset insensitivity (band flatness) of the qubit subspace. For
more details of a specific metric, see Patel  _et al._ (2024). An ideal regime
corresponds to $\alpha_{2}\gg\alpha_{1}>1$ with $\Phi\sim\pi/2$ — a
demonstration is provided in Fig. 1(c), which compares the different qubit
subspaces formed taking $\\{\alpha_{1},\alpha_{2}\\}=\\{20,2\\}$ and with
$\Phi=\pi/2$ or $0$. Dramatic improvement, with respect to band flatness (and
not shown, also the anharmonicity) can be seen for $\Phi=\pi/2$ (non-zero
flux). The takeaway message is that, the better the control of
$\\{\alpha_{1},\alpha_{2},\Phi\\}$, the better one can design a transmon
qubit. Hence, the objective of the present work is to find setups that allow
control of $\\{\alpha_{1},\alpha_{2},\Phi\\}$, under the additional constraint
of wanting only a single junction and zero magnetic flux. This is to be
contrasted with previous proposals, which include the triple junction d-mon
subject to magnetic flux Patel  _et al._ (2024) and the non-ideal $\pi$-SQUID:
a superconducting loop formed by two $\pi$-periodic circuit elements, also
subject to an external magnetic flux Staples  _et al._ (2024).

##  IV Results: Twisted triplet-triplet junction

As our specific setup, we consider the triplet-triplet junction, i.e. with
triplet vectors $\bm{d}^{(1)}$ and $\bm{d}^{(2)}$, and further impose an easy
axis spin-orbit coupling, via a TMD substrate layer, see e.g. Scammell and
Scheurer (2024) for detailed discussion. The purpose is to reduce the SO(3)
spin-rotational symmetry of the $\bm{d}^{(1,2)}$, and couple it to a given
spatial (crystal) axis, i.e. $\bm{d}^{(1)}=d_{0}(1,0,0)$ in subsystem 1. Next,
spatially rotating subsystem 2 to an angle $\Phi_{J}$ relative to subsystem 1,
implies that $\bm{d}^{(2)}=d_{0}(\cos\Phi_{J},\sin\Phi_{J},0)$; the triplet
angle follows the spatial twist angle, thanks to SOC. The total gap function
entering the BdG Hamiltonian (4) is then

$\displaystyle{\Delta}^{(1)}_{\bm{k}}+e^{i\varphi}{\Delta}^{(2)}_{\bm{k}}=\Upsilon_{\bm{k}}(\bm{d}^{(1)}+e^{i\varphi}\bm{d}^{(2)})\cdot\bm{{s}}({s}_{y}{\tau}_{y}).$
(8)

In what follows, we denote
$\bm{d}\equiv\bm{d}^{(1)}+e^{i\varphi}\bm{d}^{(2)}$.

Finally, allowing for spin polarization — denoted $\bm{M}=M_{z}{z}$ and taken
to be along $z$, i.e. appearing in (3) as $M_{z}s_{z}$ — the free energy
expansion then permits the terms,

$\displaystyle\delta{\cal U}_{1}$ $\displaystyle\sim
i\bm{M}\cdot(\bm{d}^{*}\times\bm{d})=g_{1}M_{z}\sin\Phi_{J}\sin\varphi,$ (9)
$\displaystyle\delta{\cal U}_{2}$ $\displaystyle\sim
i\bm{M}\cdot(\bm{d}^{*}\times\bm{d})(\bm{d}^{*}\cdot\bm{d})=g_{2}M_{z}\sin
2\Phi_{J}\sin 2\varphi.$

Here $g_{1,2}$ are dimensionless parameters, that can be obtained via direct
computation of (III). The $\sin\varphi$ and $\sin 2\varphi$ contributions to
the transmon potential is our key observation. Their existence relies on the
vector product $(\bm{d}^{*}\times\bm{d})\neq\bm{0}$, i.e. non-unitary
superconductivity Sigrist and Ueda (1991), as well as the vector $\bm{M}$
associated here with spin polarization.

Armed with the above setup, and upon expanding to first-order in $M_{z}$, we
arrive at the parameters of the potential (III) (i.e. before re-gauging to
(7)),

$\displaystyle a_{1}$ $\displaystyle\sim a_{1}^{0}\cos\Phi_{J},\quad b_{1}\sim
g_{1}M_{z}\sin\Phi_{J},$ (10) $\displaystyle a_{2}$ $\displaystyle\sim
a_{2}^{0}\cos 2\Phi_{J},\quad b_{2}\sim g_{2}M_{z}\sin 2\Phi_{J}.$

Here $a_{1}^{0},a_{2}^{0}$ are independent of $M_{z}$ to order $O(M_{z}^{2})$.
Upon re-gauging, we arrive at the transmon Hamiltonian (7), with parameters

$\displaystyle\alpha_{1}$
$\displaystyle=\sqrt{a_{1}^{2}+b_{1}^{2}},\quad\alpha_{2}=\sqrt{a_{2}^{2}+b_{2}^{2}},$
(11) $\displaystyle\Phi$
$\displaystyle=\arctan\left(\frac{g_{1}M_{z}}{a^{0}_{1}}\tan\Phi_{J}\right)-\frac{1}{2}\arctan\left(\frac{g_{2}M_{z}}{a^{0}_{2}}\tan
2\Phi_{J}\right).$

This constitutes the title claim: that the junction twist angle, $\Phi_{J}$,
generates the anomalous magnetic flux $\Phi$, via a nonlinear relationship
(11).

##  V Discussion

We examined c-axis tunnel junctions built from triplet superconductors. Our
objective was to demonstrate that such junctions are useful for creating
highly optimal transmon qubits, based on the criteria of anharmonicity and
offset charge insensitivity. Leaning on two recent analyses Patel  _et al._
(2024); Staples  _et al._ (2024), we distilled the problem statement down to:
how can we arrange the triplet superconductors, so as to have control over the
three relevant parameters in the transmon Hamiltonian,
$\\{\alpha_{1},\alpha_{2},\Phi\\}$. In our work, these parameters are
understood to control, respectively, single and double CP tunneling and an
anomalous (zero-field) magnetic flux that acts solely on the single tunneling
processes.

To achieve our goal, we found that a spin-triplet junction, with rotated
triplet vectors on either side, and in the presence of spin polarization, gave
us the desired Hamiltonian. The ability to rotate the spin vector became tied
to the ability to rotate the junction, thanks to spin-orbit coupling.

All of these ingredients: triplet superconductivity, spin polarization, spin-
orbit interaction and a valley degree of freedom, are features of moiré
graphene (with a TMD substrate), which prompted our searches in the first
instance. Moreover, taking moiré graphene as the qubit platform also gives us
hope that the present work will inspire near-term experimental tests, given
the ongoing research activity into moiré graphene. Aside from moiré graphene,
we note that triplet superconductivity in a valley system has been modeled via
modulated gating of a two-dimensional electron (or hole) gas Li  _et al._
(2020, 2021); Scammell  _et al._ (2022b), which can also be designed to induce
spin orbit coupling Sushkov and Castro Neto (2013); Scammell and Sushkov
(2019). This represents another possible avenue.

Finally, we contrast with recent works Díez-Mérida  _et al._ (2023); de Vries
_et al._ (2021); Rodan-Legrain  _et al._ (2021), which have achieved gate-
defined Josephson junctions within moiré graphene (see also recent theory Xie
_et al._ (2023)); such junctions are in-plane, whereas we consider the out-of-
plane (c-axis) junction. We leave it for future work to explore the ideas
presented here within a gate defined, in-plane junction.

## Acknowledgements

We thank Julian Ingham for critical feedback on the manuscript.

## References

  * Arute and et al (2019) F. Arute and et al, “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505 (2019).
  * Kjaergaard  _et al._ (2020) M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, “Superconducting qubits: Current state of play,” Annual Review of Condensed Matter Physics 11, 369 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050605 .
  * Gyenis  _et al._ (2021) A. Gyenis, A. Di Paolo, J. Koch, A. Blais, A. A. Houck, and D. I. Schuster, “Moving beyond the transmon: Noise-protected superconducting quantum circuits,” PRX Quantum 2, 030101 (2021).
  * Fowler  _et al._ (2012) A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, “Surface codes: Towards practical large-scale quantum computation,” Phys. Rev. A 86, 032324 (2012).
  * Reiher  _et al._ (2017) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, “Elucidating reaction mechanisms on quantum computers,” Proceedings of the National Academy of Sciences 114, 7555 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1619152114 .
  * Acharya and et al (2023) R. Acharya and et al, “Suppressing quantum errors by scaling a surface code logical qubit,” Nature 614, 676 (2023).
  * Patel  _et al._ (2024) H. Patel, V. Pathak, O. Can, A. C. Potter, and M. Franz, “$d$-mon: A transmon with strong anharmonicity based on planar $c$-axis tunneling junction between $d$-wave and $s$-wave superconductors,” Phys. Rev. Lett. 132, 017002 (2024).
  * Staples  _et al._ (2024) J. W. Staples, T. B. Smith, and A. C. Doherty, “Universal flux-based control of a $\pi$-squid,” (2024), arXiv:2312.04321 [quant-ph] .
  * Kitaev (2006) A. Kitaev, “Protected qubit based on a superconducting current mirror,” (2006), arXiv:cond-mat/0609441 [cond-mat.mes-hall] .
  * Paolo  _et al._ (2019) A. D. Paolo, A. L. Grimsmo, P. Groszkowski, J. Koch, and A. Blais, “Control and coherence time enhancement of the 0-$\pi$ qubit,” New Journal of Physics 21, 043002 (2019).
  * Ciaccia  _et al._ (2024) C. Ciaccia, R. Haller, A. C. C. Drachmann, T. Lindemann, M. J. Manfra, C. Schrade, and C. Schönenberger, “Charge-4e supercurrent in a two-dimensional inas-al superconductor-semiconductor heterostructure,” Communications Physics 7, 41 (2024).
  * Valentini  _et al._ (2024) M. Valentini, O. Sagi, L. Baghumyan, T. de Gijsel, J. Jung, S. Calcaterra, A. Ballabio, J. Aguilera Servin, K. Aggarwal, M. Janik, T. Adletzberger, R. Seoane Souto, M. Leijnse, J. Danon, C. Schrade, E. Bakkers, D. Chrastina, G. Isella, and G. Katsaros, “Parity-conserving cooper-pair transport and ideal superconducting diode in planar germanium,” Nature Communications 15, 169 (2024).
  * Fominov  _et al._ (2003) Y. V. Fominov, A. A. Golubov, and M. Y. Kupriyanov, “Decoherence due to nodal quasiparticles in d-wave qubits,” Journal of Experimental and Theoretical Physics Letters 77, 587 (2003).
  * Naimer  _et al._ (2021) T. Naimer, K. Zollner, M. Gmitra, and J. Fabian, “Twist-angle dependent proximity induced spin-orbit coupling in graphene/transition metal dichalcogenide heterostructures,” Phys. Rev. B 104, 195156 (2021).
  * Li and Koshino (2019) Y. Li and M. Koshino, “Twist-angle dependence of the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides,” Phys. Rev. B 99, 075438 (2019).
  * David  _et al._ (2019) A. David, P. Rakyta, A. Kormányos, and G. Burkard, “Induced spin-orbit coupling in twisted graphene–transition metal dichalcogenide heterobilayers: Twistronics meets spintronics,” Phys. Rev. B 100, 085412 (2019).
  * Lee  _et al._ (2022) S. Lee, D. J. P. de Sousa, Y.-K. Kwon, F. de Juan, Z. Chi, F. Casanova, and T. Low, “Charge-to-spin conversion in twisted $\mathrm{graphene}/{\mathrm{wse}}_{2}$ heterostructures,” Phys. Rev. B 106, 165420 (2022).
  * Péterfalvi  _et al._ (2022) C. G. Péterfalvi, A. David, P. Rakyta, G. Burkard, and A. Kormányos, “Quantum interference tuning of spin-orbit coupling in twisted van der waals trilayers,” Phys. Rev. Research 4, L022049 (2022).
  * Zollner and Fabian (2021) K. Zollner and J. Fabian, “Bilayer graphene encapsulated within monolayers of ${\mathrm{WS}}_{2}$ or ${\mathrm{Cr}}_{2}{\mathrm{Ge}}_{2}{\mathrm{Te}}_{6}$: Tunable proximity spin-orbit or exchange coupling,” Phys. Rev. B 104, 075126 (2021).
  * Zollner  _et al._ (2022) K. Zollner, M. Gmitra, and J. Fabian, “Proximity spin-orbit and exchange coupling in aba and abc trilayer graphene van der waals heterostructures,” Phys. Rev. B 105, 115126 (2022).
  * Park  _et al._ (2021) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, “Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene,” Nature 590, 249 (2021).
  * Hao  _et al._ (2021) Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, “Electric field-tunable superconductivity in alternating-twist magic-angle trilayer graphene,” Science 371, 1133 (2021).
  * Lin  _et al._ (2022) J.-X. Lin, P. Siriviboon, H. D. Scammell, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, M. S. Scheurer, and J. I. A. Li, “Zero-field superconducting diode effect in small-twist-angle trilayer graphene,” Nature Physics (2022).
  * Siriviboon  _et al._ (2021) P. Siriviboon, J.-X. Lin, X. Liu, H. D. Scammell, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, M. S. Scheurer, and J. I. A. Li, “A new flavor of correlation and superconductivity in small twist-angle trilayer graphene,” arXiv e-prints (2021), arXiv:2112.07127 [cond-mat.mes-hall] .
  * Xu and Balents (2018) C. Xu and L. Balents, “Topological superconductivity in twisted multilayer graphene,” Phys. Rev. Lett. 121, 087001 (2018).
  * You and Vishwanath (2019) Y.-Z. You and A. Vishwanath, “Superconductivity from valley fluctuations and approximate so(4) symmetry in a weak coupling theory of twisted bilayer graphene,” npj Quantum Materials 4, 16 (2019).
  * Scheurer and Samajdar (2020) M. S. Scheurer and R. Samajdar, “Pairing in graphene-based moiré superlattices,” Phys. Rev. Research 2, 033062 (2020).
  * Lake  _et al._ (2022) E. Lake, A. S. Patri, and T. Senthil, “Pairing symmetry of twisted bilayer graphene: A phenomenological synthesis,” Phys. Rev. B 106, 104506 (2022).
  * Scammell  _et al._ (2022a) H. D. Scammell, J. I. A. Li, and M. S. Scheurer, “Theory of zero-field superconducting diode effect in twisted trilayer graphene,” 2D Mater. 9, 025027 (2022a).
  * Scammell and Scheurer (2023) H. D. Scammell and M. S. Scheurer, “Tunable superconductivity and möbius fermi surfaces in an inversion-symmetric twisted van der waals heterostructure,” Phys. Rev. Lett. 130, 066001 (2023).
  * Scammell and Scheurer (2024) H. D. Scammell and M. S. Scheurer, “Displacement field tunable superconductivity in an inversion-symmetric twisted van der waals heterostructure,” Phys. Rev. B 109, 035159 (2024).
  * Sigrist and Ueda (1991) M. Sigrist and K. Ueda, “Phenomenological theory of unconventional superconductivity,” Rev. Mod. Phys. 63, 239 (1991).
  * Li  _et al._ (2020) T. Li, J. Ingham, and H. D. Scammell, “Artificial graphene: Unconventional superconductivity in a honeycomb superlattice,” Phys. Rev. Res. 2, 043155 (2020).
  * Li  _et al._ (2021) T. Li, M. Geier, J. Ingham, and H. Scammell, “Higher-order topological superconductivity from repulsive interactions in kagome and honeycomb systems,” 2D Materials (2021).
  * Scammell  _et al._ (2022b) H. D. Scammell, J. Ingham, M. Geier, and T. Li, “Intrinsic first- and higher-order topological superconductivity in a doped topological insulator,” Phys. Rev. B 105, 195149 (2022b).
  * Sushkov and Castro Neto (2013) O. P. Sushkov and A. H. Castro Neto, “Topological insulating states in laterally patterned ordinary semiconductors,” Phys. Rev. Lett. 110, 186601 (2013).
  * Scammell and Sushkov (2019) H. D. Scammell and O. P. Sushkov, “Tuning the topological insulator states of artificial graphene,” Phys. Rev. B 99, 085419 (2019).
  * Díez-Mérida  _et al._ (2023) J. Díez-Mérida, A. Díez-Carlón, S. Y. Yang, Y.-M. Xie, X.-J. Gao, J. Senior, K. Watanabe, T. Taniguchi, X. Lu, A. P. Higginbotham, K. T. Law, and D. K. Efetov, “Symmetry-broken josephson junctions and superconducting diodes in magic-angle twisted bilayer graphene,” Nature Communications 14, 2396 (2023).
  * de Vries  _et al._ (2021) F. K. de Vries, E. Portolés, G. Zheng, T. Taniguchi, K. Watanabe, T. Ihn, K. Ensslin, and P. Rickhaus, “Gate-defined josephson junctions in magic-angle twisted bilayer graphene,” Nature Nanotechnology 16, 760 (2021).
  * Rodan-Legrain  _et al._ (2021) D. Rodan-Legrain, Y. Cao, J. M. Park, S. C. de la Barrera, M. T. Randeria, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, “Highly tunable junctions and non-local josephson effect in magic-angle graphene tunnelling devices,” Nature Nanotechnology 16, 769 (2021).
  * Xie  _et al._ (2023) Y.-M. Xie, D. K. Efetov, and K. T. Law, “${\varphi}_{0}$-josephson junction in twisted bilayer graphene induced by a valley-polarized state,” Phys. Rev. Res. 5, 023029 (2023).


# Contact lines on stretched soft solids: Modeling anisotropic surface
stresses

Stefanie Heyden<EMAIL_ADDRESS>Department of Materials, ETH
Zürich, 8093 Zürich, Switzerland.  Nicolas Bain  Department of Materials, ETH
Zürich, 8093 Zürich, Switzerland.  Qin Xu  Department of Physics, The Hong
Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong.  Robert W. Style  Department of Materials, ETH Zürich, 8093 Zürich,
Switzerland.  Eric R. Dufresne<EMAIL_ADDRESS>Department of
Materials, ETH Zürich, 8093 Zürich, Switzerland.

###### Abstract

We present fully analytical solutions for the deformation of a stretched soft
substrate due to the static wetting of a large liquid droplet, and compare our
solutions to recently published experiments (Xu _et al_ , _Soft Matter_ 2018).
Following a Green’s function approach, we extend the surface-stress
regularized Flamant-Cerruti problem to account for uniaxial pre-strains of the
substrate. Surface profiles, including the heights and opening angles of
wetting ridges, are provided for linearized and finite kinematics. We fit
experimental wetting ridge shapes as a function of applied strain using two
free parameters, the surface Lamé coefficients. In comparison with
experiments, we find that observed opening angles are more accurately captured
using finite kinematics, especially with increasing levels of applied pre-
strain. These fits qualitatively agree with the results of Xu _et al_ , but
revise values of the surface elastic constants.

##  I Introduction

Capillary phenomena of liquids are familiar to us from everyday life and have
a long history in theoretical mechanics [1]. There, surface stresses dominate
the mechanics of liquids, leading to a range of phenomena from the spherical
shape of small water droplets to capillary action within narrow tubes. Their
effects are visible to the naked eye because their characteristic scale is set
by the capillary length $L_{c}=\sqrt{\gamma/\rho
g}\sim\mathcal{O}(\mathrm{mm})$, which relates the surface energy, $\gamma$,
to density, $\rho$, and gravitational acceleration, $g$. In solids, capillary
phenomena dominate at length scales much shorter than the elastocapillary
length, $l_{e}\sim\Upsilon/E$. Here, $\Upsilon$ is the surface stress and $E$
denotes the material’s Young’s modulus. For stiff engineering materials
($E\sim\mathcal{O}(\mathrm{GPa})$), the elastocapillary length is comparable
to the atomic scale. Thus, capillary phenomena are rarely significant in stiff
materials. Soft materials ($E\sim\mathcal{O}(\mathrm{kPa})$) feature
elastocapillary lengths that can reach $\mathcal{O}(10~{}\mathrm{\mu m})$.
Thus, capillary phenomena can become relevant in a variety of settings, _e.g._
, the wetting behavior on soft substrates [2, 3, 4, 5, 6], flattening of
nearly plane soft solids [7], stiffening effects in composites with liquid
inclusions [8], or soft adhesion [8].

Solid surface stresses have several features that distinguish them from liquid
surface stresses. Due to the ability of liquid molecules to freely rearrange,
the surface energy of simple liquids is independent of deformation.
Consequently, the surface stress takes on a strain-independent isotropic form,
$\Upsilon_{ij}=\gamma^{l}~{}\delta_{ij}$. Thus, the terms ‘surface tension’
and ‘surface energy’ are used interchangeably in that context [9]. Solids on
the other hand exhibit much more complex surface stresses, including, _e.g._ ,
anisotropy [10, 11] and inelasticity [12, 13]. Our understanding of the
underlying physical principles governing surface stresses in soft solids,
however, is only weakly developed and constitutes an active subject of
research [14, 15, 16, 17, 18, 19, 20, 21].

Figure 1: Sketch of the static wetting problem under consideration,
illustrating opening angle $\alpha$ of the wetting ridge, contact angle
$\theta$, as well as in-plane angle $\psi$ between an applied pre-strain
$\epsilon_{xx}$ and the contact line generated by liquid surface tension
$\gamma_{l}$.

To shed light on surface stresses in soft solids, recent experiments have used
static wetting techniques [22, 9, 15, 16]. There, liquid droplets were
deposited on soft substrates, and the microscopic geometry of the contact-line
(where the three interfaces meet) was precisely measured. As illustrated on
the left of Figure 1, surface stresses were inferred from a local force
balance, ignoring any contributions from bulk elasticity. The net force per
unit length of the contact line was assumed to have contributions from the
surface tensions of the liquid-vapor interface, $\gamma^{l}\hat{t}_{i}^{sl}$,
the solid-liquid interface, $\Upsilon_{ij}^{sl}\hat{t}_{j}^{sl}$, and the
solid-vapor interface, $\Upsilon_{ij}^{sv}\hat{t}_{j}^{sv}$. Here,
$\Upsilon_{ij}$ is the surface stress tensor and $\hat{t}_{j}$ is the tangent
vector of each interface. To gain insights into the strain-dependence of
surface stresses, recent experiments [15, 16] have applied the same approach
to stretched substrates as shown in Figure 1 on the right. In these
experiments, the substrate is stretched in both uniaxial [15] as well as
biaxial [16] set-ups preceding droplet deposition.

Data from [16] is shown in Figure 2, which depicts the measured static wetting
profiles within the direction of applied pre-strain and perpendicular to the
direction of applied pre-strain. The slope of the surface converges to a
finite value below the elastocapillary length since the tip behavior of
wetting ridges is dominated by surface stresses [23]. Thus, the tip geometry
is well-described by the opening angle, $\alpha$. As shown in [16], the
opening angles are strain-dependent. A balance of surface stresses at the tip
suggests that the surface stresses are also strain dependent. An alternate
explanation of these data has emerged recently, based on nonlinearities within
the bulk material, propagated to singular points at the surface [18].

Analytical solutions to the static wetting problem trace back to classical
problems of point- and line forces on an elastic half space commonly referred
to as the Boussinesq- and Flamant problem, respectively [24, 25]. In addition,
Cerruti’s problem considers a point force that is tangent to the surface of a
half-space [26]. Solutions were subsequently extended to account for
anisotropic half-spaces [27, 28], as well as isotropic media of finite
thickness and possibly multiple layers [29, 30, 31].

Here, we extend the analytical solution to the static wetting problem in the
large droplet limit to furthermore account for pre-strains of the surface. To
predict the fully three-dimensional profile of the wetting ridge, solutions
account for varying in-plane angles between the attendant line force and the
applied pre-strain. To gain further insight into the characteristics of the
strain-dependence of surface stresses, the theory is laid out in both
frameworks of linearized and finite kinematics and predicted opening angles of
the wetting ridge are compared to available experimental data. We compare our
results to published experiments, finding good agreement for data with revised
elastic constants.

Figure 2: Top: Experimental static wetting profiles within the direction of
applied pre-strain (A) and perpendicular to the direction of applied pre-
strain (B). Bottom: Zoom into the surface stress dominated region of
experimental static wetting profiles and corresponding fittings. The data
shown here is a subset of experimental static wetting profiles first published
in [16]

.

##  II Formulation of the problem

To model substrate deformations upon static wetting in the limit of large
droplets, we consider a substrate of finite thickness, $h$, as shown in Figure
1. The substrate is assumed to be neutrally-wetted, such that
$\Upsilon_{ij}^{sv}=\Upsilon_{ij}^{sl}$. For convenience, both will be
referred to as $\Upsilon_{ij}$ in the following. We further allow for an
emergent anisotropy in surface stresses due to pre-strains applied to the
substrate. Solutions account for varying in-plane angles $\psi$ between the
applied pre-strain and the in-plane direction perpendicular to the contact
line, thus predicting the fully three-dimensional surface profile of a large
sessile droplet.

###  II.1 Bulk elasticity

We follow a Green’s function approach in studying the effect of a line force
acting upon an elastic substrate within the linearized kinematics framework.
Starting from linear momentum balance, insertion of the Cauchy relations for
an isotropic linear elastic material gives the well-known Navier-Lamé
equations

$\displaystyle(\lambda+2\mu)u_{x,xx}+(\lambda+\mu)(u_{y,yx}+u_{z,xz})+\mu(u_{x,yy}+u_{x,zz})$
$\displaystyle=0$
$\displaystyle(\lambda+2\mu)u_{y,yy}+(\lambda+\mu)(u_{x,xy}+u_{z,yz})+\mu(u_{y,xx}+u_{y,zz})$
$\displaystyle=0$ (1)
$\displaystyle(\lambda+2\mu)u_{z,zz}+(\lambda+\mu)(u_{x,xz}+u_{y,zy})+\mu(u_{z,xx}+u_{z,yy})$
$\displaystyle=0,$

whereby gravitational effects are neglected. For the limiting case of large
droplets, the line force generated by the liquid surface tension $\gamma_{l}$
may be approximated as a straight contact line as illustrated in Figure 1.
Under the assumption of plane strain within the in-plane direction parallel to
the contact line, the set of Equations (1) simplifies to

$\displaystyle(\lambda+2\mu)u_{x,xx}+(\lambda+\mu)u_{z,xz}+\mu u_{x,zz}$
$\displaystyle=0$
$\displaystyle(\lambda+2\mu)u_{z,zz}+(\lambda+\mu)u_{x,xz}+\mu u_{z,xx}$
$\displaystyle=0.$ (2)

Application of Fourier transformations in $x$ leads to a system of ordinary
differential equations, which can be recast into matrix form as

$\begin{pmatrix}\hat{u}_{x,z}\\\ \hat{u}_{z,z}\\\ \hat{u}_{x,zz}\\\
\hat{u}_{z,zz}\end{pmatrix}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\
\frac{(\lambda+2\mu)}{\mu}k^{2}&0&0&-\frac{i(\lambda+\mu)}{\mu}k\\\
0&\frac{\mu}{(\lambda+2\mu)}k^{2}&-\frac{i(\lambda+\mu)}{(\lambda+2\mu)}k&0\end{pmatrix}\begin{pmatrix}\hat{u}_{x}\\\
\hat{u}_{z}\\\ \hat{u}_{x,z}\\\ \hat{u}_{z,z}\end{pmatrix}$

with wave number $k$. In shorthand notation, a solution to the system

$\frac{\partial}{\partial z}\hat{u}_{j}=A_{ij}\hat{u}_{j}$ (3)

may be written in terms of matrix exponentials

$\hat{u}_{j}=e^{A_{ij}z}\hat{u}^{0}_{j},$ (4)

whereby we find the Jordan normal form of $A_{ij}$ to compute the matrix
exponential $B_{ij}=e^{A_{ij}z}$ (the full form of $B_{ij}$ is given in the
appendix). Here, $\hat{u}^{0}_{j}$ denotes boundary conditions applied at the
bottom of the substrate in the form $\hat{u}_{x}^{0}=0$ and
$\hat{u}_{z}^{0}=0$. The solution may be reformulated as

$\hat{u}_{i}=B_{ij}\hat{u}^{0}_{j}=\begin{pmatrix}B_{1}&B_{2}\\\
B_{3}&B_{4}\end{pmatrix}\hat{u}^{0}_{j},$ (5)

with $\hat{u}_{i}=(\hat{u}_{x},\hat{u}_{z},\hat{u}_{x,z},\hat{u}_{z,z})$ and
$\hat{u}^{0}_{i}=(0,0,\hat{u}^{0}_{x,z},\hat{u}^{0}_{z,z})$. Each block
$B_{i}$ furthermore constitutes a $2\times 2$ matrix, with

$\begin{pmatrix}\hat{u}_{x}\\\
\hat{u}_{z}\end{pmatrix}=B_{2}\begin{pmatrix}\hat{u}_{x,z}^{0}\\\
\hat{u}_{z,z}^{0}.\end{pmatrix}$ (6)

Using this relation, we can replace $\hat{u}^{0}_{j}$ with
$(0,0,B_{2}^{-1}(\hat{u}_{x},\hat{u}_{z}))$, which gives

$\hat{u}_{i}=\begin{pmatrix}B_{2}\,B_{2}^{-1}\\\
B_{4}\,B_{2}^{-1}\end{pmatrix}\begin{pmatrix}\hat{u}_{x}\\\
\hat{u}_{z}\end{pmatrix}.$ (7)

In the following, we aim to solve for substrate displacements at the surface
$z=h$ only. Note that in the absence of uniaxial pre-strains
$\epsilon_{xx}^{pre}$, $h=h_{0}$, whereas $h=h(1-\epsilon_{xx}^{pre}/2)$)
otherwise. Tractions applied at the surface due to the line force generated by
the liquid surface tension $\gamma_{l}$ take the form

$\begin{pmatrix}\sigma^{h}_{xz}\\\
\sigma^{h}_{zz}\end{pmatrix}=\begin{pmatrix}\gamma_{l}\,\cos({\theta})\delta(x)\\\
\gamma_{l}\,\sin({\theta})\delta(x)\end{pmatrix}\Rightarrow\begin{pmatrix}\hat{\sigma}^{h}_{xz}\\\
\hat{\sigma}^{h}_{zz}\end{pmatrix}=\begin{pmatrix}\gamma_{l}\,\cos({\theta})\\\
\gamma_{l}\,\sin({\theta})\end{pmatrix},$ (8)

where $\theta$ is the contact angle between the applied line force and the
substrate’s surface. Under the assumption of isotropic linear elasticity,
_e.g._ , $\sigma_{ij}=\lambda\epsilon_{kk}\delta{ij}+2\mu\epsilon_{ij}$,
tractions and displacements may be related via a coefficient matrix $S$ such
that

$\displaystyle\begin{pmatrix}\hat{\sigma}^{h}_{xz}\\\
\hat{\sigma}^{h}_{zz}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}S_{11}&S_{12}&S_{13}&S_{13}\\\
S_{21}&S_{22}&S_{23}&S_{23}\end{pmatrix}\hat{u}^{h}_{i}$
$\displaystyle=\begin{pmatrix}0&\frac{iE}{2(1+\nu)}k&\frac{E}{2(1+\nu)}&0\\\
-\frac{iE\nu}{(1+\nu)(2\nu-1)}k&0&0&\frac{E(\nu-1)}{(1+\nu)(2\nu-1)}\end{pmatrix}\begin{pmatrix}B_{2}\,B_{2}^{-1}\\\
B_{4}\,B_{2}^{-1}\end{pmatrix}\begin{pmatrix}\hat{u}^{h}_{x}\\\
\hat{u}^{h}_{z}\end{pmatrix}.$ (9)

In shorthand notation, we write Equation 9 as

$\hat{\tau}^{h}_{i}=\hat{Q}_{ij}\hat{u}^{h}_{j},$ (10)

and refer to $\hat{Q}_{ij}$ as a generic matrix of spring constants relating
tractions $\hat{\tau}^{h}_{i}$ applied at the substrate’s surface to surface
displacements.

###  II.2 Surface stresses

We now proceed to furthermore take into account surface stresses which are
counteracting tractions applied at the substrate’s surface following [32, 33].
Starting from a force balance at the surface, we have

$[\hat{\sigma}_{ij}\hat{n}_{j}]^{+}_{-}=-(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})\hat{\Upsilon}_{ik,j},$
(11)

where $\hat{\sigma}_{ij}$ denotes the Cauchy stress tensor within the bulk
material, $\hat{\Upsilon}_{ik}$ represents the two-dimensional surface stress
tensor and $\hat{n}_{i}$ is the surface normal in the undeformed
configuration. For an isotropic linear elastic surface stress tensor,
$\Upsilon_{ik}=\Upsilon^{0}\delta_{ik}+\lambda^{s}\epsilon_{mm}\delta_{ik}+2\mu^{s}\epsilon_{ik}$,
Equation 11 gives

$\displaystyle[\hat{\sigma}^{h}_{xz}]^{+}_{-}$
$\displaystyle=[\hat{t_{i}}\hat{\sigma}_{ij}\hat{n}_{j}]^{+}_{-}$
$\displaystyle=-(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})(\hat{\Upsilon}_{ik}\hat{t}_{k})_{,j}+(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})\hat{\Upsilon}_{ik}\hat{t}_{k,j}$
$\displaystyle=-(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})(\hat{\Upsilon}_{ik}\hat{t}_{k})_{,j}$
$\displaystyle=k^{2}(\lambda^{s}+2\mu^{s})\hat{u}^{h}_{x}$ (12)

to first order in $\hat{u}_{i}$, with wave number $k$ [34].

To derive surface stresses counteracting normal tractions
$[\hat{\sigma}^{h}_{zz}]^{+}_{-}$ in an equivalent way, we first note that

$\displaystyle[\hat{\sigma}^{h}_{zz}]^{+}_{-}$
$\displaystyle=[\hat{n_{i}}\hat{\sigma}_{ij}\hat{n}_{j}]^{+}_{-}$
$\displaystyle=-(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})(\hat{\Upsilon}_{ik}\hat{n}_{k})_{,j}+(\delta_{ij}-\hat{n}_{i}\hat{n}_{j})\hat{\Upsilon}_{ik}\hat{n}_{k,j}$
$\displaystyle=0.$ (13)

We therefore derive an expression for the surface normal $\hat{\eta}_{i}$ in a
perturbed configuration as

$\displaystyle\hat{\eta}_{i}$ $\displaystyle=\begin{pmatrix}-\hat{u}_{z,x}\\\
-\hat{u}_{z,y}\\\ 1\end{pmatrix}+\mathcal{O}(u^{2}).$ (14)

Using this relation, surface stresses counteracting tractions
$[\hat{\sigma}^{h}_{zz}]^{+}_{-}$ may be derived as

$\displaystyle[\hat{\sigma}^{h}_{zz}]^{+}_{-}$
$\displaystyle=[\hat{\eta_{i}}\hat{\sigma}_{ij}\hat{\eta}_{j}]^{+}_{-}$
$\displaystyle=(\delta_{ij}-\hat{\eta}_{i}\hat{\eta}_{j})\hat{\Upsilon}_{ik}\hat{\eta}_{k,j}$
$\displaystyle=k^{2}\Upsilon^{0}\hat{u}^{h}_{z}+\mathcal{O}(u^{2}).$ (15)

Generalizing Equation 10 to the additional action of surface stresses
therefore results in

$\hat{\tau}^{h}_{i}=(\hat{Q}_{ij}+\hat{T}_{ij})\hat{u}^{h}_{j},\quad\text{with}\quad\hat{T}_{ij}=\begin{pmatrix}k^{2}(\lambda^{s}+2\mu^{s})&0\\\
0&k^{2}\Upsilon^{0}\end{pmatrix}.$ (16)

###  II.3 Anisotropic surface stresses based on pre-strain

For an isotropic substrate subjected to uniaxial tension in $x$ and free
boundaries in $y$ and $z$, strains are related as

$\epsilon^{pre}_{yy}=-\nu\epsilon^{pre}_{xx},$ (17)

with $\nu$ being the bulk’s Poisson ratio. Resultant surface stresses follow
as

$\displaystyle\Upsilon^{pre}_{xx}$
$\displaystyle=(\lambda^{s}+2\mu^{s}-\nu\lambda^{s})\epsilon^{pre}_{xx}\quad\text{and}$
$\displaystyle\Upsilon^{pre}_{yy}$ $\displaystyle=(\lambda^{s}-\nu
2\mu^{s}-\nu\lambda^{s})\epsilon^{pre}_{xx},$ (18)

such that total surface stresses are given by
$\Upsilon_{ik}=\Upsilon^{0}\delta_{ik}+\Upsilon^{pre}_{ik}+\lambda^{s}\epsilon_{mm}\delta_{ik}+2\mu^{s}\epsilon_{ik}$.
To derive the fully three-dimensional profile of the wetting ridge, we
introduce an arbitrary in-plane rotation $\psi$ connecting reference $(x,y)$
and rotated $(x^{\prime},y^{\prime})$ frames. Surface normals in the perturbed
configuration are related as

$\hat{\eta}_{i}=R(\psi)_{ij^{{}^{\prime}}}\hat{\eta}_{i^{{}^{\prime}}}=\begin{pmatrix}-\cos{(\psi)}\hat{u}_{z^{{}^{\prime}},x^{{}^{\prime}}}\\\
-\sin{(\psi)}\hat{u}_{z^{{}^{\prime}},x^{{}^{\prime}}}\\\ 1\end{pmatrix},$
(19)

whereby the corresponding gradient is given by

$\displaystyle\hat{\eta}_{i,j}$
$\displaystyle=\hat{\eta}_{i,k^{{}^{\prime}}}R(\psi)_{jk^{{}^{\prime}}}$ (20)
$\displaystyle=\begin{pmatrix}-\cos^{2}{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}x^{{}^{\prime}}}&-\cos{(\psi)}\sin{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}x^{{}^{\prime}}}&-\cos{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}z^{{}^{\prime}}}\\\
-\sin{(\psi)}\cos{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}x^{{}^{\prime}}}&-\sin^{2}{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}x^{{}^{\prime}}}&-\sin{(\psi)}u_{z^{{}^{\prime}},x^{{}^{\prime}}z^{{}^{\prime}}}\\\
0&0&0\end{pmatrix}.$

Extending Equation 16 to take into account anisotropy in surface stresses
based on an applied pre-strain hence results in

$\displaystyle\hat{\tau}^{h}_{i}$
$\displaystyle=(\hat{Q}_{ij}+\hat{T}_{ij}+\hat{T}^{pre}_{ij})\hat{u}^{h}_{j},\quad\text{with}$
$\displaystyle\hat{T}^{pre}_{ij}$ $\displaystyle=\begin{pmatrix}0&&&0\\\
0&&&k^{2}(\Upsilon_{xx}^{pre}\cos^{2}{(\psi)}+\Upsilon_{yy}^{pre}\sin^{2}{(\psi))}\end{pmatrix}.$
(21)

###  II.4 Extension to finite deformation

In the following, we extend the model to account for finite pre-strains of the
surface and subsequently linearize about the pre-strained configuration
preceding droplet deposition. We choose a two-dimensional compressible Neo-
Hookean strain energy density

$W^{NH}=\frac{\mu^{s}}{2}\left(\frac{I_{1}}{\text{det}(\mathbf{F})}-2\right)+\frac{\kappa^{s}}{2}(\text{det}(\mathbf{F})-1)^{2},$
(22)

where $I_{1}=\text{tr}(\mathbf{F}^{T}\mathbf{F})$ denotes the first tensor
invariant [35]. Linearization about a given finite deformation gives the
tangent matrix in spatial form as

$C_{ijkl}=\text{det}(\mathbf{F})^{-1}F_{jJ}F_{lL}C_{iJkL},$ (23)

with $C_{iJkL}=\partial(\partial W^{NH}/\partial F_{iJ})/\partial F_{kL}$
being the tangent stiffness matrix. Note the dependence on deformation
gradient $\mathbf{F}$, which in turn renders the resulting tangent matrix
anisotropic in the present case of a given uniaxial pre-strain
$\mathbf{F}=((1+\epsilon_{xx}^{pre},0),(0,1-\nu^{s}\epsilon_{xx}^{pre}))$.
Resultant surface stresses upon subsequent droplet deposition follow as
previously derived in Equations (12) and (15), whereby now
$\hat{\Upsilon}_{ij}=C_{ijkl}\epsilon_{kl}$.

###  II.5 Calculated Surface Profiles

Figure 3 illustrates theoretically predicted wetting profiles within the
linearized kinematics framework at different strain levels. The wetting ridge
decreases in height with increasing pre-strain levels due to a stiffer
material response of the surface. This stiffening effect is more pronounced
within the direction of applied pre-strain at $\Psi=0$. Similarly, opening
angles of the wetting ridge increase with increasing pre-strain, showing a
more distinct increase at $\Psi=0$.

Figure 3: Static wetting on substrates exhibiting anisotropy in surface
stresses based on applied pre-strains (linearized kinematics framework). Left:
Wetting ridges along in-plane angle $\Psi$ under varying pre-strains of
$\epsilon_{xx}=0.0$, $\epsilon_{xx}=0.125$ and $\epsilon_{xx}=0.25$. Right:
Corresponding wetting profiles at $\Psi=0$ and $\Psi=\pi/2$ with inset
illustrating the increase in opening angles with increasing pre-strain.
Simulations employ values of $\Upsilon^{0}=0.02$ N/m, $\mu^{s}=0.012$ N/m,
$\lambda^{s}=0.055$ N/m and $E=2.88$ kPa.

##  III Comparison with Experiments

We now compare the the results of these calculations to experiments of Xu _et
al_. To minimize systematic errors due to variations in material composition
between different samples, we use an experimental subset of the data first
published in [16] spanning one sample over a range of pre-strain levels at
positions $A,B$ at $\Psi=0,\pi/2$, respectively (see Figure 1).

Figure 4: Profiles of substrate profiles far from the contact line with
corresponding their corresponding best fits to the theory described here, with
$E=2.88$ kPa. Individual curves are shifted by an offset $u_{\text{off}}$ for
better visibility.

First, we fit the deformation of the substrate far from the contact line to
determine the substrate elastic modulus. In this region, where the curvature
is small, the liquid surface tension, $\gamma^{l}=0.041~{}\mathrm{N/m}$ [15],
is balanced by bulk elastic forces with contribution from solid capillary
forces. Fitting ranges are chosen from $20$ $\mu m$ to the final data point in
each experimental data set. We find that far-field surface profiles for $4$
values of the pre-strain at both locations A and B can be simultaneously best-
fit with a Young’s modulus of $E=2.88$ kPa. Figure 4 depicts the tail behavior
of experimental static wetting profiles and the corresponding fittings.

Figure 5: Opening angles of the wetting ridge as predicted from linearized and
finite kinematics theory in comparison to opening angles measured in static
wetting experiments. Values are depicted for different pre-strain levels
$\epsilon$ and compared both within the direction of applied pre-strain at
$\Psi=0$ (A) as well as perpendicular to it at $\Psi=\pi/2$ (B).

Subsequently, surface elastic moduli $\lambda^{s}$ and $\mu^{s}$ are
determined from a comparison of experimental and analytical opening angles at
the tip of the wetting ridge, as surface stresses dominate the material
behavior below the elastocapillary length. Opening angles are calculated by
fitting surface profiles within $10$ $\mu$m of either side of the ridge tip to
quadratic polynomials. The position $x_{0}$ of the tip of the wetting ridge is
included as a free parameter during the fitting procedure to account for
possible shifts in experimental profiles. Fitted profiles are depicted in
Figure 2, and their opening angles are depicted in Fig. 5. We note that all
tip profiles are well-fit within the linear regime.

To find best-fit values of surface elastic moduli, we seed the solution space
with test points using a grid of varying discretization. Figure 6 illustrates
reduced chi-square statistics $\chi^{2}$ applied to opening angles observed in
experiments and simulations within the tested solution space. Lowest contour
levels $\chi^{2}=0.95$ indicate that differences in opening angles obtained
from experiments and simulations are in accordance with the experimental error
variance. Resultant best-fit values are obtained as $\lambda^{s}=55\pm 1$ mN/m
and $\mu^{s}=12\pm 1$ mN/m (with standard deviations corresponding to a $95\%$
confidence interval). Here, we fixed $\Upsilon^{0}=0.02$ N/m and
$\gamma^{l}=0.041$ N/m as reported in previous investigations [15, 16]. This
gives a surface Poisson’s ratio,
$\nu^{s}=\lambda^{s}/(2(\lambda^{s}+\mu^{s}))=0.41$. In two dimensions,
Poisson’s ratio must lie on the interval $\nu^{s}\in[0,1]$. While our specific
values of $\lambda_{s}$ and $\mu_{s}$ differ from the values in [16], our
value of the biaxial stretching modulus $\Lambda=2(\lambda^{s}+\mu^{s})$ is
consistent from previously reported values [15, 16].

Figure 6: Contour plot showing reduced chi-square statistics $\chi^{2}$
applied to opening angles observed in experiments and simulations for varying
surface elastic moduli. The marker located at
$(\mu_{s}=20\,\text{mN/m},\lambda_{s}=43\,\text{mN/m})$ illustrates $\chi^{2}$
for surface elastic moduli published in [16].

A comparison of the experimental and best-fit opening angles are plotted in
Figure 5 as a function of applied strain. As expected, both theories of
linearized and finite kinematics predict the same opening angle at the tip of
the wetting ridge for $\epsilon_{xx}=0$. In the direction perpendicular to the
applied pre-strain, opening angles are recovered within experimental error
measures in intermediate as well as high deformation regimes. Within the
direction of applied pre-strain, however, linearized kinematics theory
consistently overpredicts opening angles of the wetting ridge due to an
inherent overstiffening of the material response. Using a finite kinematics
model of surface elasticity softens the material response, resulting in
decreased values of opening angles. This effect is more pronounced with
increasing levels of applied pre-strain. Note that the chosen Neo-Hookean
strain energy density allows for direct comparison to linearized kinematics
theory, as both models of surface elasticity may be formulated in terms of
model parameters $\mu^{s}$ and $\lambda^{s}$. Experimental opening angles have
furthermore been compared to calculated values based on a two dimensional Yeoh
strain energy density, for which results only display a marginal difference to
the chosen Neo-Hookean model.

$\epsilon_{xx}^{\text{pre}}$ | 0.0 | 0.13 | 0.22 | 0.25
---|---|---|---|---
$\alpha^{A}_{\text{exp,sub}}$ | 94 $\pm$ 3 | 109 $\pm$ 5 | 119 $\pm$ 6 | 126 $\pm$ 6
$\alpha^{A}_{\text{theo,lin}}$ | 96 | 117 | 126 | 129
$\alpha^{A}_{\text{theo,fin}}$ | 96 | 115 | 123 | 126
$\alpha^{B}_{\text{exp,sub}}$ | 93 $\pm$ 5 | 106 $\pm$ 7 | 113 $\pm$ 7 | 120 $\pm$ 9
$\alpha^{B}_{\text{theo,lin}}$ | 96 | 107 | 114 | 116
$\alpha^{B}_{\text{theo,fin}}$ | 96 | 106 | 113 | 115
Table 1: Opening angles $\alpha^{exp,sub}$ as determined from an experimental
subset of data originally published in [16] and theoretically determined
values $\alpha_{\text{theo,lin}}$ and $\alpha_{\text{theo,fin}}$.

For a detailed comparison, Table 1 lists values of both analytically
determined opening angles within linearized and finite kinematics theory
($\alpha_{theo,lin}$ and $\alpha_{theo,fin}$, respectively) as well as
previously reported experimental results $\alpha_{exp}$ [16]. In addition,
Table 1 lists opening angles calculated from an experimental subset of the
data first published in [16], which is used in determining best-fit surface
elastic moduli $\lambda^{s}$ and $\mu^{s}$.

##  IV Conclusions

We have described a theory of strain-dependent surface stresses in the
framework of static wetting by extending the generalized Flamant-Cerruti
problem. Analytical solutions provide fully three-dimensional surface profiles
near a contact line, with an arbitary orientation relative to the direction of
applied stretch. We found that finite kinematics more accurately capture
experimental opening angles, especially with increasing levels of applied pre-
strain.

To minimize systematic errors, we limited our analysis to a subset of
experiments from [16], where a series of surface profiles are measured on a
single substrate as a function of strain and location. While we recover re-
analysis of these data returns the same value of the surface bi-axial
stretching modulus, we find different values of the surface Lamé coefficients.

The origins of surface elasticity in gels are unknown. One simple model
described in [9] suggests that surface elasticity of gels can arise due to a
thin layer of material near the surface with a modulus much higher than the
bulk.

In this picture, 2D Lamé moduli (denoted as $s$) are related to the 3D Lamé of
the putative stiff layer (denoted as $l$) through
$(\mu^{s},\lambda^{s})=(\mu^{l}h,\frac{2\mu^{l}\lambda^{l}}{\lambda^{l}+2\mu^{l}}h)$,
where $h$ is the thickness of the stiff layer [36]. Substituting our measured
values of $\lambda^{s}$ and $\mu^{s}$, we find
$\frac{\mu^{l}}{\lambda^{l}}=-0.28$. Since elastic moduli must be positive,
this is an unphysical result. This suggests that the surface elastic moduli do
not come from such a simple layer, or that nonlinear effects at the contact
line [18] need to be explicitly accounted for.

Possible extensions of the theory presented in this work include a re-
formulation of the problem to allow for bulk elasticity within finite
deformations. This extension may leverage further insights into the strain-
dependence of surface stresses within the regime of large pre-deformations of
the elastic substrate. An examination of different constitutive models of
surface elasticity may furthermore give insights into the specific
characteristics of strain dependence. As a secondary measure of comparison,
theoretical results could be compared to local strains at the contact line.
Current discrepancies between the predicted and measured in-plane
displacements are based on the influence of the Laplace pressure difference of
the droplet, and could hence be resolved by extending the theory to droplets
of arbitrary size. Future experiments on a wider range of materials with
different interfacial properties are furthermore needed to differentiate
between the two hypotheses of strain-dependent surface stresses versus the
effects of bulk nonlinearities being propagated to surfaces of constant
surface energy [18]. Possible experimental set-ups in future investigations
could include a droplet deposited on a stretched soft solid substrate immersed
in a second liquid phase. Varying second liquid phases would thus allow
testing of a wide range of surface tensions at the contact line.

##  V Appendix

The full form of components $B_{1}$ through $B_{4}$ of the matrix exponential
$B_{ij}$ in equation (5) is given as

$B_{1}=\begin{pmatrix}\text{cosh}(kz)+\frac{kz\,\text{sinh}(kz)}{(2-4\nu)}&\frac{i(kz\,\text{cosh}(kz)-\text{sinh}(kz))}{4(-1+\nu)}\\\
\frac{2i(kz\,\text{cosh}(kz)-\text{sinh}(kz))}{8\nu-4}&\text{cosh}(kz)+\frac{kz\,\text{sinh}(kz)}{(4\nu-4)}\\\
\end{pmatrix}$ (24)
$B_{2}=\begin{pmatrix}\frac{kz\,\text{cosh}(kz)+3\,\text{sinh}(kz)-4\nu\,\text{sinh}(kz)}{4k-4k\nu}&\frac{2iz\,\text{sinh}(kz)}{8\nu-4}\\\
\frac{iz\,\text{sinh}(kz)}{4(\nu-1)}&\frac{kz\,\text{cosh}(kz)-3\,\text{sinh}(kz)+4\nu\,\text{sinh}(kz)}{4k\nu-2k}\\\
\end{pmatrix}$ (25)
$B_{3}=\begin{pmatrix}\frac{k(-kz\,\text{cosh}(kz)+(4\nu-3)\,\text{sinh}(kz))}{4\nu-2}&\frac{ik^{2}z\,\text{sinh}(kz)}{4(\nu-1)}\\\
\frac{2ik^{2}z\,\text{sinh}(kz)}{8\nu-4}&\frac{k(kz\,\text{cosh}(kz)+(4\nu-3)\,\text{sinh}(kz))}{4(\nu-1)}\end{pmatrix}$
(26)
$B_{4}=\begin{pmatrix}\text{cosh}(kz)+\frac{kz\,\text{sinh}(kz)}{(4-4\nu)}&\frac{2i(kz\,\text{cosh}(kz)+\text{sinh}(kz))}{8\nu-4}\\\
\frac{i(kz\,\text{cosh}(kz)+\text{sinh}(kz))}{4(-1+\nu)}&\text{cosh}(kz)+\frac{kz\,\text{sinh}(kz)}{(4\nu-2)}\end{pmatrix}$
(27)

## References

  * Maxwell [1878] J. C. Maxwell, Encylopaedia Britannica 5 (1878).
  * Carré  _et al._ [1996] A. Carré, J.-C. Gastel, and M. E. Shanahan, Nature 379, 432 (1996).
  * Extrand and Kumagai [1996] C. Extrand and Y. Kumagai, Journal of Colloid and Interface Science 184, 191 (1996).
  * de Gennes  _et al._ [2004] P.-G. de Gennes, F. Brochard-Wyart, and D. Quéré,  _Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves_ (Springer, New York, 2004).
  * Roman and Bico [2010] B. Roman and J. Bico, Journal of Physics: Condensed Matter 22, 493101 (2010).
  * Yu [2012] Y.-s. Yu, Applied Mathematics and Mechanics 33, 1095 (2012).
  * Jagota  _et al._ [2012] A. Jagota, D. Paretkar, and A. Ghatak, Phys. Rev. E 85, 051602 (2012).
  * Style  _et al._ [2015] R. W. Style, J. S. Wettlaufer, and E. R. Dufresne, Soft Matter 11, 672 (2015).
  * Style  _et al._ [2017] R. Style, A. Jagota, C.-Y. Hui, and E. Dufresne, Annu. Rev. Condens. Matter Phys. 8, 99 (2017).
  * Stein  _et al._ [2018] P. Stein, A. Moradabadi, M. Diehm, B.-X. Xu, and K. Albe, Acta Materialia 159, 225 (2018).
  * Dreher  _et al._ [2019] T. Dreher, N. Pineau, E. Bourasseau, P. Malfreyt, L. Soulard, and C. Lemarchand, J Chem Phys 151, 244703 (2019).
  * Altenbach [2012] H. Altenbach, International Journal of Engineering Science 59, 83 (2012).
  * Pepicelli  _et al._ [2019] M. Pepicelli, N. Jaensson, C. Tregouët, B. Schroyen, A. Alicke, T. Tervoort, C. Monteux, and J. Vermant, Journal of Rheology 63, 815 (2019).
  * Jensen  _et al._ [2017] K. E. Jensen, R. W. Style, Q. Xu, and E. R. Dufresne, Phys. Rev. X 7, 041031 (2017).
  * Xu  _et al._ [2017] Q. Xu, K. Jensen, R. Boltyanskiy, R. Safarti, R. Style, and E. Dufresne, Nature Communications 555(8), 2041 (2017).
  * Xu  _et al._ [2018] Q. Xu, R. Style, and E. Dufresne, Soft Matter 14(916), 916 (2018).
  * Schulman  _et al._ [2018] R. Schulman, M. Trejo, T. Salez, E. Raphael, and K. Dalnoki-Veress, Nature Communications 9, 982 (2018).
  * Masurel  _et al._ [2019] R. Masurel, M. Roché, L. Limat, I. Ionescu, and J. Dervaux, Phys. Rev. Lett. 122, 248004 (2019).
  * Bostwick  _et al._ [2014] J. Bostwick, M. Shearer, and K. Daniels, Soft Matter 10, 7361 (2014).
  * Cao  _et al._ [2014] Z. Cao, M. Stevens, and A. Dobrynin, Macromolecules 18, 6515–6521 (2014).
  * Andreotti and Snoeijer [2020] B. Andreotti and J. Snoeijer, Annual Review of Fluid Mechanics 52, 285 (2020).
  * Style  _et al._ [2013] R. W. Style, R. Boltyanskiy, Y. Che, J. S. Wettlaufer, L. A. Wilen, and E. R. Dufresne, Phys. Rev. Lett. 110, 066103 (2013).
  * Wu  _et al._ [2018] H. Wu, Z. Liu, A. Jagota, and C.-Y. Hui, Soft Matter 14, 1847 (2018).
  * Boussinesque [1885] J. Boussinesque,  _Applications des potentiels à l’étude de l’équilibre et mouvement des solides elastiques_ (Gauthier–Villard, Paris, 1885).
  * Flamant [1892] A. Flamant, CR Acad Sci 114, 1465–1468 (1892).
  * Cerruti [1882] V. Cerruti, Rend Accad Lincei 3(13), 81–122 (1882).
  * Sveklo [1964] V. Sveklo, PMM 28(5), 908 (1964).
  * Willis [1967] J. Willis, J Mech Phys Solids 15, 331 (1967).
  * Burmister [1945] D. M. Burmister, Journal of Applied Physics 16, 89 (1945).
  * Yue [1996] Z. Yue, Journal of Mechanics and Applied Mathematics 49, 471 (1996).
  * Merkel  _et al._ [2007] R. Merkel, N. Kirchgeßner, C. M. Cesa, and B. Hoffmann, Biophysical Journal 93, 3314 (2007).
  * Jerison  _et al._ [2011] E. Jerison, Y. Xu, L. Wilen, and E. Dufresne, Physical Review Letters 106, 186103 (2011).
  * Style and Dufresne [2012] R. Style and E. Dufresne, Soft Matter 27(8), 7177 (2012).
  * Style and Xu [2018] R. Style and Q. Xu, Soft Matter 14, 4569 (2018).
  * Holzapfel [2000] G. Holzapfel,  _Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science_ (John Wiley & Sons, Chichester, 2000).
  * Benveniste and Miloh [2007] Y. Benveniste and T. Miloh, J Elasticity 88, 87 (2007).