# Cranked Skyrme–Hartree–Fock–Bogoliubov for a mean-field description of
nuclear rotations near the drip line

Kenichi Yoshida<EMAIL_ADDRESS>Department of Physics,
Kyoto University, Kyoto, 606-8502, Japan

###### Abstract

To describe the yrast states in weakly-bound nuclei, I directly solve the
coordinate-space cranked Skyrme–Hartree–Fock–Bogoliubov equation in the
quasiparticle basis with the continuum states discretized in a box. After the
numerical demonstration for the ground-state band in a medium-mass nucleus, I
apply the newly-developed method to neutron-rich even-$N$ Mg isotopes. I find
that the appearance of the significantly low-energy $I^{\pi}=2^{+}$ state in
40Mg is mainly due to the suppression of pairing. The calculation predicts
that the $2^{+}$ state in 42Mg appears as high in energy as in
${}^{34\text{--}38}$Mg whereas the triaxial deformation is enhanced. The
present numerical framework offers a practical approach for investigating the
near yrast states and revealing structures unique in drip-line nuclei.

††preprint: KUNS-2893

A diverse variety of modes of motion show up in nuclei. To explore unique
structures in exotic nuclei, various spectroscopic studies have been carried
out via decays and making use of direct reactions. Nuclear deformation is a
collective phenomenon that has long attracted interest; a nucleus is deformed
as it moves away from the magic number due to the many-body correlation. It
has been turned out that the shell structure evolves as the conventional magic
numbers disappear and new ones appear instead, depending on the combination of
the particle numbers of protons and neutrons [1]. Energies of excited nuclear
states are often among the first quantities accessible in experiments and have
been used as an indicator of the changing shell structure and the onset of
deformation [2].

Rotational motion is a manifestation of the spontaneous breaking of the
rotational symmetry [3]. The cranking approximation provides an intuitive
picture of the rotation of quantum systems, and the cranked shell model is a
standard tool to investigate the microscopic structure of nuclear rotations
near the ground state and at high spins [4, 5]. The cranked shell model in the
framework of the nuclear energy-density functional (EDF) has provided a
systematic and quantitative description of the (near) yrast states from light
to heavy nuclei [6]. The accumulation of experimental data for low spins in
neutron-rich nuclei has stimulated the application of the EDF-based cranked
shell model to the elucidation of the rotational motion unique in weakly bound
nuclei. Furthermore, it has been pointed out in Ref. [7] that the moments of
inertia are sensitive to the pairing and weak binding.

The island of inversion has been the subject of much experimental and
theoretical interest [8]. A systematic calculation for the Mg isotopes in the
mean-field approximation produces a spherical configuration in 32Mg, a soft
potential energy surface in 34Mg, and a prolate configuration in 36,38,40Mg
[9]. It has been clarified that the shape fluctuation and the correlation
beyond the mean-field approximation is significant in 32Mg [10, 11, 12].
Experimentally, the measurement of not only the first $I^{\pi}=2^{+}$ state
but the $4^{+}$ state and their ratio $R_{4/2}>3$ have revealed a well-
deformed structure in 34,36,38Mg [13]. The significantly low energy of the
$2^{+}$ state in 40Mg is not reproduced by any theoretical models, and it
indicates the emergence of a unique feature associated with weak binding [14].
As a coherent contribution of the pairing in the continuum states induces an
enhanced transition strength to the low-lying vibrational mode in the Mg
isotopes [15, 16, 17, 18, 19], the roles of the weak binding and the continuum
coupling in the rotational motions are interesting to study.

I investigate in this article the low-spin states in the neutron-rich Mg
isotopes close to the drip line. Then, I try to clarify the role of weak
biding in the low-lying excited states. To this end, I develop a new framework
of the cranked shell model within the nuclear EDF approach, which is capable
of handling nuclides with arbitrary mass numbers [6, 20]. To describe the
(near) yrast states with proper account of the pairing in the continuum
states, I directly solve the coordinate-space cranked
Skyrme–Kohn–Sham–Bogoliubov (KSB) or Hartree–Fock–Bogoliubov equation in the
quasiparticle basis:

$\displaystyle\sum_{\sigma^{\prime}}\begin{bmatrix}h^{q\prime}_{\sigma\sigma^{\prime}}(\bm{r})&\tilde{h}^{q}_{\sigma\sigma^{\prime}}(\bm{r})\\\
4\sigma\sigma^{\prime}\tilde{h}^{q*}_{-\sigma-\sigma^{\prime}}(\bm{r})&-4\sigma\sigma^{\prime}h^{q\prime*}_{-\sigma-\sigma^{\prime}}(\bm{r})\end{bmatrix}\begin{bmatrix}\varphi^{q}_{1,\alpha}(\bm{r}\sigma^{\prime})\\\
\varphi^{q}_{2,\alpha}(\bm{r}\sigma^{\prime})\end{bmatrix}$
$\displaystyle=E^{\prime}_{\alpha}\begin{bmatrix}\varphi^{q}_{1,\alpha}(\bm{r}\sigma)\\\
\varphi^{q}_{2,\alpha}(\bm{r}\sigma)\end{bmatrix},$ (1)

which is obtained by extending the formalism developed for describing the
ground-state properties of even-even nuclei near the drip line [21]. Here the
single-particle Routhian and the pair Hamiltonian are defined by using the
Skyrme and paring EDF $E[\rho,\tilde{\rho},\tilde{\rho}^{*}]$ as
$h^{q\prime}_{\sigma\sigma^{\prime}}(\bm{r})=\frac{\delta
E[\rho,\tilde{\rho},\tilde{\rho}^{*}]}{\delta\rho^{q}_{\sigma^{\prime}\sigma}(\bm{r})}-(\lambda^{q}+\omega_{\rm
rot}j_{z})\delta_{\sigma\sigma^{\prime}}$ and
$\tilde{h}^{q}_{\sigma\sigma^{\prime}}(\bm{r})=\frac{\delta
E[\rho,\tilde{\rho},\tilde{\rho}^{*}]}{\delta\tilde{\rho}^{q*}_{\sigma^{\prime}\sigma}(\bm{r})}$.
I define the $z$-axis as a quantization axis of the intrinsic spin and
consider the system rotating uniformly about the $z$-axis. I take the natural
units: $\hbar=c=1$.

The numerical procedure to solve Eq. (1) is basically the same as in solving
the cranked KS equation in Refs. [22, 23]. I impose the reflection symmetry
about the $(x,y)$-, $(y,z)$\- and $(z,x)$-planes. Thus, the parity
$\mathfrak{p}_{k}$ $(=\pm 1)$ and $z$-signature $r_{k}$ ($=\pm\mathrm{i}$) are
a good quantum number. I solve Eq. (1) by diagonalizing the KSB Hamiltonian in
the three-dimensional (3D) Cartesian-mesh representation with the box boundary
condition. Thanks to the reflection symmetries, I have only to consider
explicitly the octant region in space with $x\geq 0$, $y\geq 0$, and $z\geq
0$; see Refs. [24, 25] for details. I use a 3D lattice mesh
$x_{i}=ih-h/2,y_{j}=jh-h/2,z_{k}=kh-h/2\ \ (i,j,k=1,2,\cdots M)$ with a mesh
size $h=1.0$ fm and $M=10$ for each direction. To check the convergence of the
results with respect to the box size and to investigate the effect of the weak
binding, I change $M$ in the discussion below. The differential operators are
represented by the use of the 9-point formula of the finite difference method.
For diagonalization of the matrix of Eq. (1), I use the LAPACK dsyevx
subroutine [26]. A modified Broyden’s method [27] is utilized to calculate new
densities during the selfconsistent iteration. The quasiparticle Routhians are
cut off at 60 MeV.

To see the validity of the present framework, I perform the calculation for
the ground-state rotational band in 172Hf as a typical example of the
collective rotation. Figure 1 shows the calculated kinematic moments of
inertia $\mathcal{J}^{(1)}$ as functions of the rotational frequency. I
employed the SLy4 [28] and the Yamagami–Shimizu–Nakatsukasa (YSN) pairing EDF
in Ref. [29]. The present model describes well the low spin states and the
band crossing. Around $\omega_{\rm rot}=0.25$ MeV, the alignment of neutrons
in the $i_{13/2}$ orbital occurs, whereas this is lower than the measurement
$\sim 0.3$ MeV. One sees that the rotational property beyond the band crossing
is also reasonably described.

Figure 1:  Kinematic moments of inertia $\mathcal{J}^{(1)}$ as functions of
the rotational frequency. The experimental data [30] are denoted by open
symbols.  Figure 2:  Similar to Fig. 1 but for 34Mg. The results obtained by
using several pairing interactions are displayed. The experimental data are
obtained from Ref. [31].  Figure 3:  Absolute value of the pairing gaps
$|\Delta|$ and pairing energies $E_{\rm pair}$ of neutrons in 34Mg.

Then, I investigate the low-spin yrast states in the neutron-rich Mg isotopes.
I use the SkM* functional [32] and the so-called mixed-type pairing
interaction with the strength $V_{0}=-295$ MeV fm3 as in Ref. [18]. The
ground-state properties thus obtained are summarized in Table 1 of Ref. [18].
I found that the triaxiality is negligibly small: $\gamma<1^{\circ}$ in low
spins. Because the chemical potential is a key quantity governing the spatial
structure of quasiparticle wave functions, I put it down:
$-4.17,-3.25,-2.41,-1.60$ MeV for neutrons in 34,36,38,40Mg, respectively.
When the two-basis method, in which the KSB Hamiltonian is diagonalized in a
truncated single-particle basis [6], is employed for such weakly-bound nuclei,
the single-particle continuum states enter the pairing window. Therefore, the
convergence with respect to the number of basis states has to be carefully
examined because the densities should be spatially localized as far as
$\lambda<0$. In the present case, however, the densities are always calculated
to be localized because the full KSB Hamiltonian is directly diagonalized
[21].

Figure 2 shows the calculated $\mathcal{J}^{(1)}$ of 34Mg, and compares with
the experimental data [31]. The measured $R_{4/2}$ value is 3.06, which is
marginally lower than that of the rigid rotor. The calculation reproduces well
the slight increase in $\mathcal{J}^{(1)}$ due to the weakening of pairing.
The role of the density dependence of the pairing interaction has been
discussed in the study of the superdeformed states [33]: the density
dependence results in a retarded alignment. I investigate here the density
dependence of the pairing interaction in the light nuclei. To this end, I use
the volume- and surface-type pairing interactions. I determined the strength
to keep the calculated pairing gap as obtained with the mixed pairing:
$V_{0}=-207$ and $-451$ MeV fm3 for the volume and surface pairing,
respectively. The pairing gap is defined as $\Delta^{q}=\int
d\bm{r}\tilde{h}^{q}(\bm{r})\rho^{q}(\bm{r})/\int d\bm{r}\rho^{q}(\bm{r})$ and
the resultant gap is shown in Fig. 3. Notice that the protons are unpaired at
$\omega_{\rm{rot}}=0$. In low $\omega_{\rm rot}$, the volume pairing also
reproduces the experimental value. However, it gives a faster increase in
$\mathcal{J}^{(1)}$ similarly to the finding in Ref. [33]. The surface pairing
shows a slow increase in $\mathcal{J}^{(1)}$. An interesting feature is that
the magnitude of $\mathcal{J}^{(1)}$ is inversely proportional to the
magnitude of the pairing energy $E^{q}_{\rm pair}=1/2\int
d\bm{r}\tilde{h}^{q}(\bm{r})\tilde{\rho}^{q*}(\bm{r})$, as shown in Fig. 3.

Figure 4:  (a) Similar to Fig. 1 but for 36,38,40Mg. The experimental data,
open symbols, are taken from Refs. [13, 14]. (b) Calculated
$\mathcal{J}^{(1)}$ for 40Mg obtained with the box size of $M=8,10$, and 12.
Results without the pairing are also depicted.

Next, I investigate the rotational property of 36,38,40Mg located close to the
drip line. Figure 4(a) displays the calculated $\mathcal{J}^{(1)}$ together
with the experimental data [13, 14]. The calculated $\mathcal{J}^{(1)}$ for
36,38Mg in low $\omega_{\rm{rot}}$ is similar to the one for 34Mg, and smaller
than that for 40Mg. This is consistent with the calculation in Ref. [18],
where the moments of inertia were evaluated using the Thouless–Valatin
procedure in the framework of the Skyrme EDF-based QRPA [34]. A higher value
for $\mathcal{J}^{(1)}$ of 40Mg is partly because of smaller deformation than
others. Another reason is the weak-binding effect, as discussed below. The
experimental data indicate that the pairing in 36,38Mg would be weaker than
the calculation because the measured $\mathcal{J}^{(1)}$ is slightly larger
than the calculated one. Furthermore, the measurement shows a faster increase
in $\mathcal{J}^{(1)}$ than the calculation for 38Mg; the calculation produces
a stronger pairing correlation. As in the case of 34Mg, I found a good
correspondence between $\mathcal{J}^{(1)}$ and $E_{\rm{pair}}$.

Figure 5:  Similar to Fig. 1 but for 34,36,38,40Mg obtained using the YSN
pairing functional of Ref. [29].

The present calculation fails to describe the rotational motion in 40Mg: the
calculated $\mathcal{J}^{(1)}$ is far below the measured value. One is then
tempted to expect that the weak-binding effect shows up. In Fig. 4(b), I show
the results obtained by varying the box size. With the increase in the box
size, the calculated $\mathcal{J}^{(1)}$ increases. The calculation with $M=8$
does not show the convergence, while the difference between the results with
$M=10$ and 12 becomes small. One cannot expect a further increase in
$\mathcal{J}^{(1)}$ even enlarging the box size. Notice that for the case of
34Mg, the results obtained with $M=8$ and 10 are not very different, as shown
in Fig. 2. In this sense, the weak-binding effect appears in 40Mg; the weak
binding of neutrons reduces the pairing correlation and enhances the moment of
inertia. As an extreme case of weak pairing, I performed the calculation
without the pairing, as displayed in Fig. 4(b). The resultant
$\mathcal{J}^{(1)}$ is much larger than that obtained with pairing. As
increasing the rotational frequency, the results with and without pairing
become closer to each other. The observed $\mathcal{J}^{(1)}$ of 6.0 MeV-1 is
between these results and is relatively closer to the one obtained without
pairing $\mathcal{J}^{(1)}=6.8$ MeV-1. As pointed out in Ref. [9], one has an
oblate minimum in 40Mg. Indeed, I found the oblate solution with
$\beta=0.16,\gamma=60^{\circ}$. With this configuration, the resultant
$\mathcal{J}^{(1)}$ is much lower: $1.6$ MeV-1.

To describe the isotopic dependence of low-spin states in the Mg isotopes, a
key is the isospin dependence of the pairing correlation. Thus, I employ the
optimal pairing EDF, the YSN functional [29], in which the isovector-density
dependence is introduced to describe pairing in neutron-rich nuclei. Figure 5
shows the calculated $\mathcal{J}^{(1)}$ in ${}^{34\text{--}42}$Mg. Since the
coupling constant of the YSN pairing functional was fixed for the rare-earth
nuclei, it gives a weak pairing in light nuclei: the density of single-
particle states near the Fermi energy is low. The calculated
$\mathcal{J}^{(1)}$ are thus all larger than those obtained with the mixed
pairing. The pairing vanishes at as low as $\omega_{\rm rot}=0.8$ MeV and 0.5
MeV in 34Mg and 40Mg, respectively. However, the enhancement in
$\mathcal{J}^{(1)}$ from 38Mg to 40Mg is well described. One sees the
enhancement by 28% at $\omega_{\rm rot}=0.3$ MeV, which is comparable to the
experimental observation of the lowering in the $2^{+}_{1}$ energy from 38Mg
to 40Mg by 23% [14]. When the SLy4 functional is employed together with the
YSN pairing functional, both neutrons and protons are unpaired in 40Mg, while
neutrons are paired in 38Mg. Therefore, a significant lowering of the
excitation energy of $2^{+}_{1}$ state in 40Mg is also described.

When two neutrons are added, I find a further structural change. Comparing
with 40Mg, the deformation gets weaker by about 33%, and the triaxiality
develops. The pairing for protons shows up, and that for neutrons increases.
Accordingly, $\mathcal{J}^{(1)}$ is reduced as shown in Fig. 5. The yrast band
terminates at $\omega_{\rm rot}=0.8$ MeV, reaching $\mathcal{J}^{(1)}=11$
MeV-1. Therefore, the present calculation with SkM*+YSN predicts that the
irregularity in the $2^{+}_{1}$ state appears only at 40Mg. I found that the
calculated $\mathcal{J}^{(1)}$ changes by about 2% when the box size is varied
from $M=10$ to 12. Investigation of the halo structure in such a triaxialy-
deformed rotating nucleus is an interesting future work though it needs a
larger box [35].

To summarize, I have developed a numerical framework for a mean-field
description of near yrast states in nuclei near the drip line in a nuclear EDF
approach. I directly solved the coordinate-space cranked Skyrme-KSB equation
in the quasiparticle basis, with the continuum states being discretized in a
box. The present framework reproduces the low-spin states and band crossing in
a medium-heavy deformed nucleus.

The low-spin states in 34Mg are well described by using the density-dependent
pairing interaction. The calculated moments of inertia show the inverse
correlation with the magnitude of the pairing energy; the moments of inertia
increase as the magnitude of pairing energy decreases. With the increase in
the neutron number, the calculation overestimates the pairing correlation,
thus leading to the underestimation of the moments of inertia. Employing the
pairing EDF constructed to describe the isospin dependence in neutron-rich
nuclei, I have found that the appearance of the significantly low-energy
$I^{\pi}=2^{+}$ state in 40Mg is mainly due to the suppression of pairing. In
42Mg, the $2^{+}_{1}$ state appears higher in energy than in 40Mg due to the
structure change.

###### Acknowledgements.

This work was supported by the JSPS KAKENHI (Grants No. JP19K03824 and No.
JP19K03872). The numerical calculations were performed on the computing
facilities at the Yukawa Institute for Theoretical Physics, Kyoto University,
and at the Research Center for Nuclear Physics, Osaka University.

## References

  * Otsuka  _et al._ [2020] T. Otsuka, A. Gade, O. Sorlin, T. Suzuki, and Y. Utsuno, Evolution of shell structure in exotic nuclei, Rev. Mod. Phys. 92, 015002 (2020).
  * Gade [2015] A. Gade, Excitation energies in neutron-rich rare isotopes as indicators of changing shell structure, Eur. Phys. J. A 51, 118 (2015).
  * Bohr and Mottelson [1975] A. Bohr and B. Mottelson,  _Nuclear Structure: Volume II, Nuclear Deformations_ (Benjamin, Reading, MA, 1975).
  * Ring  _et al._ [1970] P. Ring, R. Beck, and H. J. Mang, On the application of the Hartree-Fock-Bogolyubov-equations to a microscopic theory of nuclear rotations, Z. Phys. A 231, 10 (1970).
  * Bengtsson  _et al._ [1986] R. Bengtsson, S. Frauendorf, and F.-R. May, Quasiparticle levels in rotating rare earth nuclei: A cranked shell-model dictionary, At. Data Nucl. Data Tab. 35, 15 (1986).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Self-consistent mean-field models for nuclear structure, Rev. Mod. Phys. 75, 121 (2003).
  * Yamagami and Shimizu [2008] M. Yamagami and Y. R. Shimizu, Pairing effects for rotational excitations unique to neutron-rich nuclei, Phys. Rev. C 77, 064319 (2008).
  * Sorlin and Porquet [2008] O. Sorlin and M.-G. Porquet, Nuclear magic numbers: New features far from stability, Prog. Part. Nucl. Phys. 61, 602 (2008).
  * Terasaki  _et al._ [1997] J. Terasaki, H. Flocard, P.-H. Heenen, and P. Bonche, Deformation of nuclei close to the two-neutron drip line in the Mg region, Nucl. Phys. A 621, 706 (1997).
  * Kimura and Horiuchi [2002] M. Kimura and H. Horiuchi, Breaking of the Neutron Magic Number $N=20$ in 32Mg and 30Ne and Its Possible Relation to the Cluster Structure, Prog. Theor. Phys. 107, 33 (2002), https://academic.oup.com/ptp/article-pdf/107/1/33/5213174/107-1-33.pdf .
  * Rodríguez-Guzmán  _et al._ [2002] R. Rodríguez-Guzmán, J. Egido, and L. Robledo, Correlations beyond the mean field in magnesium isotopes: angular momentum projection and configuration mixing, Nucl. Phys. A 709, 201 (2002).
  * Hinohara  _et al._ [2011] N. Hinohara, K. Sato, K. Yoshida, T. Nakatsukasa, M. Matsuo, and K. Matsuyanagi, Shape fluctuations in the ground and excited ${0}^{+}$ states of 30,32,34Mg, Phys. Rev. C 84, 061302 (2011).
  * Doornenbal  _et al._ [2013] P. Doornenbal, H. Scheit, S. Takeuchi, N. Aoi, K. Li, M. Matsushita, D. Steppenbeck, H. Wang, H. Baba, H. Crawford, C. R. Hoffman, R. Hughes, E. Ideguchi, N. Kobayashi, Y. Kondo, J. Lee, S. Michimasa, T. Motobayashi, H. Sakurai, M. Takechi, Y. Togano, R. Winkler, and K. Yoneda, In-Beam $\gamma$-Ray Spectroscopy of ${}^{34,36,38}\mathrm{Mg}$: Merging the $N\mathbf{=}20$ and $N\mathbf{=}28$ Shell Quenching, Phys. Rev. Lett. 111, 212502 (2013).
  * Crawford  _et al._ [2019] H. L. Crawford, P. Fallon, A. O. Macchiavelli, P. Doornenbal, N. Aoi, F. Browne, C. M. Campbell, S. Chen, R. M. Clark, M. L. Cortés, M. Cromaz, E. Ideguchi, M. D. Jones, R. Kanungo, M. MacCormick, S. Momiyama, I. Murray, M. Niikura, S. Paschalis, M. Petri, H. Sakurai, M. Salathe, P. Schrock, D. Steppenbeck, S. Takeuchi, Y. K. Tanaka, R. Taniuchi, H. Wang, and K. Wimmer, First Spectroscopy of the Near Drip-line Nucleus ${}^{40}\mathrm{Mg}$, Phys. Rev. Lett. 122, 052501 (2019).
  * Yoshida  _et al._ [2006a] K. Yoshida, M. Yamagami, and K. Matsuyanagi, Pairing and continuum effects on low-frequency quadrupole vibrations in deformed Mg isotopes close to the neutron drip line, Nucl. Phys. A 779, 99 (2006a), arXiv:nucl-th/0605073 .
  * Yoshida  _et al._ [2006b] K. Yoshida, M. Yamagami, and K. Matsuyanagi, Dynamic pairing effects on low-frequency modes of excitation in deformed Mg isotopes close to the neutron drip line, Phys. Scr. T 125, 45 (2006b), arXiv:nucl-th/0507047 .
  * Yoshida and Yamagami [2008] K. Yoshida and M. Yamagami, Low-frequency $K^{\pi}=0^{+}$ modes in deformed neutron-rich nuclei: Pairing- and $\beta$-vibrational modes of neutron, Phys. Rev. C 77, 044312 (2008), arXiv:0802.2341 [nucl-th] .
  * Yoshida [2009] K. Yoshida, Skyrme-QRPA calculations for low-lying excitation modes in deformed neutron-rich nuclei, Eur. Phys. J. A 42, 583 (2009), arXiv:0902.3053 [nucl-th] .
  * Yamagami [2019] M. Yamagami, Pairing effect on $K^{\pi}=0^{+}$ quadrupole excitations in neutron-rich Mg isotopes studied by Skyrme quasiparticle random-phase approximation calculations in wave-number space, Phys. Rev. C 100, 054302 (2019).
  * Nakatsukasa  _et al._ [2016] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Time-dependent density-functional description of nuclear dynamics, Rev. Mod. Phys. 88, 045004 (2016), arXiv:1606.04717 .
  * Dobaczewski  _et al._ [1984] J. Dobaczewski, H. Flocard, and J. Treiner, Hartree-Fock-Bogolyubov description of nuclei near the neutron-drip line, Nucl. Phys. A 422, 103 (1984).
  * Sakai  _et al._ [2020] S. Sakai, K. Yoshida, and M. Matsuo, Signature dependent triaxiality for shape evolution from superdeformation in rapidly rotating 40Ca and 41Ca, Prog. Theor. Exp. Phys. 2020, 063D02 (2020), arXiv:2003.06081 [nucl-th] .
  * Yoshida [2021] K. Yoshida, Competing shell effect of protons and neutrons for superdeformation in 60Zn, 62Zn, and 64Ge, (2021), arXiv:2108.12130 [nucl-th] .
  * Bonche  _et al._ [1987] P. Bonche, H. Flocard, and P. Heenen, Self-consistent calculation of nuclear rotations: The complete yrast line of 24Mg, Nucl. Phys. A 467, 115 (1987).
  * Ogasawara  _et al._ [2009] H. Ogasawara, K. Yoshida, M. Yamagami, S. Mizutori, and K. Matsuyanagi, Rotational Frequency Dependence of Octupole Vibrations on Superdeformed States in 40Ca, Prog. Theor. Phys. 121, 357 (2009).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Baran  _et al._ [2008] A. Baran, A. Bulgac, M. M. Forbes, G. Hagen, W. Nazarewicz, N. Schunck, and M. V. Stoitsov, Broyden’s method in nuclear structure calculations, Phys. Rev. C 78, 014318 (2008).
  * Chabanat  _et al._ [1998] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, A Skyrme parametrization from subnuclear to neutron star densities. 2\. Nuclei far from stablities, Nucl. Phys. A 635, 231 (1998), [Erratum: Nucl.Phys.A 643, 441–441 (1998)].
  * Yamagami  _et al._ [2009] M. Yamagami, Y. R. Shimizu, and T. Nakatsukasa, Optimal pair density functional for description of nuclei with large neutron excess, Phys. Rev. C 80, 064301 (2009), arXiv:0812.3197 [nucl-th] .
  * [30] National Nuclear Data Center, “Evaluated Nuclear Structure Data File”, https://www.nndc.bnl.gov/ensdf/.
  * Michimasa  _et al._ [2014] S. Michimasa, Y. Yanagisawa, K. Inafuku, N. Aoi, Z. Elekes, Z. Fülöp, Y. Ichikawa, N. Iwasa, K. Kurita, M. Kurokawa, T. Machida, T. Motobayashi, T. Nakamura, T. Nakabayashi, M. Notani, H. J. Ong, T. K. Onishi, H. Otsu, H. Sakurai, M. Shinohara, T. Sumikama, S. Takeuchi, K. Tanaka, Y. Togano, K. Yamada, M. Yamaguchi, and K. Yoneda, Quadrupole collectivity in island-of-inversion nuclei ${}^{28,30}\mathrm{Ne}$ and ${}^{34,36}\mathrm{Mg}$, Phys. Rev. C 89, 054307 (2014).
  * Bartel  _et al._ [1982] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. Håkansson, Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force, Nucl. Phys. A 386, 79 (1982).
  * Terasaki  _et al._ [1995] J. Terasaki, P.-H. Heenen, P. Bonche, J. Dobaczewski, and H. Flocard, Superdeformed rotational bands with density dependent pairing interactions, Nucl. Phys. A 593, 1 (1995).
  * Yoshida and Van Giai [2008] K. Yoshida and N. Van Giai, Deformed quasiparticle-random-phase approximation for neutron-rich nuclei using the Skyrme energy density functional, Phys. Rev. C 78, 064316 (2008), arXiv:0809.0169 [nucl-th] .
  * Uzawa  _et al._ [2021] K. Uzawa, K. Hagino, and K. Yoshida, Role of triaxiality in deformed halo nuclei, Phys. Rev. C 104, L011303 (2021), arXiv:2104.02946 [nucl-th] .


# Optimal retirement in presence of stochastic labor income: a free boundary
approach in presence of an incomplete market.

Daniele<EMAIL_ADDRESS>
Politecnico di Milano, Milano, Italy

###### Abstract

In this note, we show how to solve an optimal retirement problem in presence
of a stochastic wage dealing with a free boundary problem. In particular, we
show how to deal with an incomplete market case, where the wage cannot be
fully hedged investing in the risk-free and the risky asset describing the
financial market.

##  1 Introduction

In this note, we extend the work of [Barucci and Marazzina (2011)], where a
stochastic wage perfectly correlated to the risky asset is considered,
assuming that the wage may not be perfectly hedged exploiting the financial
market, and we show how the free boundary problem is modified due to the
incomplete market assumption, exploiting the theoretical framework in [Tepla
(2000)].

More precisely, we study a stochastic control problem involving the
consumption-portfolio-leisure policy and the optimal stopping time of
retirement. Relating to the leisure rate, the agent earns the labour income
with a stochastic wage rate. The sum of labour and leisure rates is assumed to
be constant.. As the complement, the labour rate is lower bounded by a
positive constant for the consideration of retaining the employment state
before retirement, and then the agent realizes the full leisure rate after
declaring retirement.

The problem is solved exploiting the duality method, which can be summarized
into four critical procedures. We first tackle the post-stopping time problem
to deduce a closed form of the corresponding value function. Then we apply the
Legendre-Fenchel transform to the utility function and the value function of
post-stopping time problem obtained in the first step. Afterwards, we
construct the duality between the optimal control problem with the
individual’s shadow prices problem, by the aid of the liquidity and budget
constraints and the dual transforms acquired before. Finally, we cast the dual
shadow price problem as a free boundary problem. We refer to [Barucci and
Marazzina (2011), Choi et al. (2008), Ding and Marazzina (2022), Ding and
Marazzina (2024), Dybvig and Liu (2010)] for related works. The obtained free
boundary problem can not be solved analytically, but only numerically, e.g.,
exploiting a finite difference scheme as in [Barucci and Marazzina (2011)].

The note is organized as follows: in Section 2 we deal with the considered
problem, in Section 3 we show how the duality approach is applied, finally
Section 4 contains the result of the article, that is the partial differential
equation (pde) of the related free boundary problem.

##  2 The Optimal Consumption-Investment Problem

We consider a continuous time economy. The infinite horizon life of the agent
is divided in two parts: before retirement and after retirement. In the first
part of his life the agent consumes, invests in financial markets and chooses
the labor supply rate earning an exogenous stochastic income, after retirement
he does not work, he only consumes and invests his wealth in financial
markets. After retirement the agent cannot go back to work.

The labor income process $Y(t)$ (wage) is exogenous and stochastic. More
specifically, we assume

$dY(t)=\mu_{1}Y(t)dt+\boldsymbol{\mu}^{\top}_{2}Y(t)d\boldsymbol{B}(t),\ \
Y(0)=Y_{0}>0,$ (1)

where $\mu_{1}$ is a scalar function, $\boldsymbol{\mu_{2}}$ is a $2$-vector
of functions, such that $Y(t)>0,\ \forall t>0$ and $\boldsymbol{B}(t)$ is a
$2$-vector of Brownian motions.

We assume that there are 2 assets: the risk-free asset with a constant
instantaneous interest rate $r$ and a risky assets. Asset price $S(t)$ evolve
as

$d{S}(t)={S}(t)\left({b}dt+\sigma d{B}_{1}(t)\right),\ \ S(0)=S_{0},$ (2)

where $\boldsymbol{b}$ is the constant drift of the risky asset price,
$\sigma$ is the volatility matrix. Notice that (2) can be re-written as

$d{S}(t)={S}(t)\left({b}dt+\boldsymbol{\sigma}^{\top}d\boldsymbol{B}(t)\right),\
\ S(0)=S_{0},$ (3)

with $\boldsymbol{\sigma}=[\sigma;0]$. Since the kernel of the vector
$\boldsymbol{\sigma}$ is not empty, the market is incomplete: we are not able
to perfectly hedge/insurance the wage just using the market assets. As in
[Tepla (2000)], let us define the space

$\displaystyle K(\sigma)=\Big{\\{}\boldsymbol{v}\in\mathbb{R}^{2}\ :\
\boldsymbol{\sigma}^{\top}\boldsymbol{v}=0\Big{\\}}.$ (4)

We denote by $c(t)$, $\theta(t)$ and $l(t)$ consumption, risky asset portfolio
and leisure processes, respectively. We make the following assumptions on
these processes:

  * •

$c(t)$ is a non negative progressively process measurable with respect to
${\cal F}$, such that $\int_{0}^{t}c(s)ds<\infty,\ \forall t\geq 0$ a.s.,

  * •

$\theta(t)$ represents the amount of wealth invested in the risky asset at
time $t$, it is a progressively adapted process with respect to ${\cal F}$,
such that $\int_{0}^{t}\theta^{2}(s)ds<\infty,\ \forall t\geq 0$ a.s.,

  * •

$l(t)\in\left[0,1\right]$ is the rate of leisure at time $t$, $l(t)$ is a
process measurable with respect to ${\cal F}$.

Let $\tau\geq 0$ be the retirement date chosen by the agent, i.e.,
$\tau\in{\cal S}$, where ${\cal S}$ denotes the set of all ${\cal F}$-stopping
time. The retirement decision is irreversible: during his working life the
agent chooses the leisure rate, after retirement he fully enjoys leisure. At
time $t$ the agent is endowed with one unit of “time”, before retirement he
allocates his time between labor ($1-l(t)$) and leisure ($l(t)$) with the
constraint that he has to work at least $1-L$, after retirement he only enjoys
leisure:

$0\leq l(t)\leq L<1\ \textrm{if}\ 0\leq t<\tau,\ \ \ l(t)=1\ \textrm{if}\
t\geq\tau.$

Let $W(0)=W_{0}$ be the initial wealth. The consumption-investment-leisure
strategy $(c,l,\theta)$ is admissible if it satisfies the above technical
conditions and the wealth process $W(t)$ satisfies the dynamic budget
constraint

$dW(t)\\!=\\!\left[(1-l(t))Y(t)\\!-\\!c(t)\right]dt+{\theta}(t)\left({b}\,dt+\boldsymbol{\sigma}^{\top}d\boldsymbol{B}(t)\right)+\left(\\!W(t)-\\!\theta(t)\\!\right)\\!rdt$
(5)

where $(1-l(t))Y(t)dt$ is the labor income received before retirement, i.e.,
for $0\leq t<\tau$,
${\theta}(t)\left({b}\,dt+\boldsymbol{\sigma}^{\top}d\boldsymbol{B}(t)\right)$
refers to the financial wealth evolution and $\left(W(t)-\theta(t)\right)rdt$
to the wealth invested in the money market (risk-free bond).

We impose a no borrowing condition

$W(t)\geq 0,\ \ \forall t\geq 0$ (6)

that is the agent cannot borrow money in the risk-free market against future
income.

The agent maximizes the expected utility over the infinite horizon:
preferences are time separable with exponential discounting. The instantaneous
utility at time $t$ is a function of consumption and leisure, i.e.,
$u(c(t),l(t))$. Therefore, given the initial wealth $W_{0}$ and the labor
income process $Y(t)$, we look for an admissible process triplet
$(c,l,{\theta})$ and a stopping time $\tau\in{\cal S}$ that maximize

$E\left[\int_{0}^{+\infty}e^{-\beta t}u(c(t),l(t))dt\right],$ (7)

subject to the budget constraint (5) and the no borrowing condition (6).
$\beta$ is the subjective discount rate and the utility function is given by

$u(c,l)=\frac{1}{\alpha}\frac{(l^{1-\alpha}c^{\alpha})^{1-\gamma}}{1-\gamma},\
\gamma>0,\ 0<\alpha<1.$ (8)

The infinite horizon is a strong assumption that allows us to handle
analytically the problem; to make the setting more plausible, we may allow for
an hazard rate of mortality as in [Dybvig and Liu (2010)] adding a positive
coefficient to the discount factor $\beta$.

Given the initial wealth $W_{0}$ and the labor income $Y_{0}$, the agent
maximizes the expected utility (7):

$\displaystyle{\cal J}(W_{0},Y_{0};c,l,{\theta},\tau)$ $\displaystyle:=$
$\displaystyle E\left[\int_{0}^{+\infty}e^{-\beta t}u(c(t),l(t))dt\right]$
$\displaystyle=$ $\displaystyle E\left[\int_{0}^{\tau}e^{-\beta
t}u(c(t),l(t))dt+\int_{\tau}^{+\infty}e^{-\beta t}u(c(t),1)dt\right]$

acting on $(c,l,{\theta},\tau)$ with $\tau\in{\cal S}$ and subject to the
budget constraint (5) and the no borrowing condition (6). The value function
is defined as

${\cal V}(W_{0},Y_{0}):=\sup_{(c,l,{\theta},\tau)\in{\cal A}}{\cal
J}(W_{0},Y_{0};c,l,{\theta},\tau),$ (9)

where ${\cal A}$ denotes the set of admissible strategies, i.e., processes
($c,l,{\theta},\tau$) that satisfy (5) and (6), $\tau\in{\cal S}$, such that
${\cal J}(W_{0},Y_{0};c,l,{\theta},\tau)<+\infty$.

##  3 A dual approach

As in [Tepla (2000)], considering (4) let us define the market price of risk

$\boldsymbol{\Theta}(\boldsymbol{v})\\!:=\\!\boldsymbol{\sigma}\left(\boldsymbol{\sigma}^{\top}\boldsymbol{\sigma}\right)^{-1}(b-r)+\boldsymbol{v},\
\forall\boldsymbol{v}\in K(\sigma),$

the state-price-density process (or pricing kernel)

$H(t,\boldsymbol{v}):=e^{-\left(r+\frac{1}{2}|\boldsymbol{\Theta}(\boldsymbol{v})|^{2}\right)t-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{B}(t)}$

and, for a given $t>0$, the equivalent risk neutral martingale measure

$\overline{P}(A):=E\left[e^{-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{B}(t)-\frac{1}{2}|\boldsymbol{\Theta}(\boldsymbol{v})|^{2}t}\mathbf{1}_{A}\right]\
\ \forall A\in{\cal F}_{t}.$

Under the (non-unique) risk neutral density we can redefine the budget and the
no borrowing constraint as static constraints. Following [Choi et al. (2008)],
by the optional sampling theorem we have that the budget constraint becomes

$E\left[\int_{0}^{\tau}H(t,\boldsymbol{v})\left(c(t)+Y(t)(l(t)-1)\right)dt+H(\tau,\boldsymbol{v})W(\tau)\right]\leq
W_{0},\ \ \forall\boldsymbol{v}\in K(\sigma)$ (10)

and the no borrowing constraint ($W(t)\geq 0,\ \forall t,\ 0\leq t\leq\tau$)
is equivalent to

$E_{t}\left[\int_{t}^{\tau}\frac{H(s,\boldsymbol{v})}{H(t,\boldsymbol{v})}\left(c(s)+Y(s)(l(s)-1)\right)ds+\frac{H(\tau,\boldsymbol{v})}{H(t,\boldsymbol{v})}W(\tau)\right]\geq
0$ (11)

$\forall\ t\ :\ 0\leq t\leq\tau,\ \forall\boldsymbol{v}\in K(\sigma).$ Note
that the no borrowing constraint is never binding after retirement. The
objective function can be rewritten as

$\displaystyle{\cal J}(W_{0},Y_{0};c,l,{\theta},\tau)\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!E\left[\int_{0}^{\tau}\\!e^{-\beta
t}u(c(t),l(t))dt+e^{-\beta\tau}\int_{\tau}^{+\infty}\\!\\!e^{-\beta(t-\tau)}u(c(t),1)dt\right]$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!E\left[\int_{0}^{\tau}e^{-\beta
t}u(c(t),l(t))dt+e^{-\beta\tau}U(W(\tau))\right],$

where

$\displaystyle
U(W(\tau)):=\sup_{(c,{\theta})}E\left[\int_{\tau}^{+\infty}e^{-\beta(t-\tau)}u(c(t),1)dt\right]$
(12)

is the optimal expected utility attainable at time $\tau$ with wealth
$W(\tau)$: for $t\geq\tau$ the agent solves the classical optimal consumption-
portfolio problem with initial wealth $W(\tau)$ subject to the budget
constraint. $U(W(\tau))$ is the indirect utility function associated with the
maximization problem at time $\tau$.

We define the convex conjugate of the utility function $u(c,l)$:

$\widetilde{u}(z,Y):=\max_{\scriptsize{\begin{matrix}c\geq 0\\\ 0\leq l\leq
L\end{matrix}}}u(c,l)-(c+Yl)z.$ (13)

In a similar way, we also define the convex conjugate of $U$:

$\displaystyle\widetilde{U}(z):=\sup_{w\geq 0}U(w)-wz.$ (14)

If $I$ denotes the inverse of $U^{\prime}$, i.e.,
$I(z)=\inf\\{w\,:\,U^{\prime}(w)=z\\}$, then we also have
$\widetilde{U}(z)=U(I(z))-zI(z)$.

Let $\widehat{c}$ and $\widehat{l}$ be the pair of processes that provides a
solution to (13), then the following proposition holds true.

######  Proposition 3.1.

Let

$\displaystyle\widetilde{z}=\left(\frac{\alpha}{1-\alpha}Y\right)^{\alpha(1-\gamma)-1}L^{-\gamma}.$
(15)

If $z\geq\widetilde{z}$, then

$\widetilde{u}(z,Y)=\frac{1}{\alpha}\frac{(\widehat{l}^{\,1-\alpha}\
\widehat{c}^{\,\alpha})^{1-\gamma}}{1-\gamma}-(\widehat{c}+Y\widehat{l})z,$
(16)

where $\widehat{c}=\frac{\alpha}{1-\alpha}Y\widehat{l}\ \ \textrm{and}\ \
\widehat{l}=z^{-\frac{1}{\gamma}}\left(\frac{\alpha}{1-\alpha}Y\right)^{\frac{\alpha(1-\gamma)-1}{\gamma}}.$
If $z<\widetilde{z}$, then

$\widetilde{u}(z,Y)=\frac{\left(\widehat{c}^{\,\alpha}L^{1-\alpha}\right)^{(1-\gamma)}}{\alpha(1-\gamma)}-(\widehat{c}+YL)z,$
(17)

where
$\widehat{c}=\left(zL^{(\alpha-1)(1-\gamma)}\right)^{\frac{1}{\alpha(1-\gamma)-1}}.$

###### Proof.

See [Barucci and Marazzina (2011), Proposition 3.1]. ∎

Let $\lambda>0$ be a Lagrange Multiplier. We consider a non-increasing process
$D(t)>0$ with $D(0)=1$ to take into account the no bankruptcy constraint, as
in [Choi et al. (2008), He and Pages (1993)], i.e., $D(t)$ is the integral of
the shadow price of the no bankruptcy constraint. Then, the following
proposition holds.

######  Proposition 3.2.

Let

$\displaystyle\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v})\\!\\!\\!\\!$
$\displaystyle:=$ $\displaystyle\\!\\!\\!\\!E\Big{[}\int_{0}^{\tau}e^{-\beta
t}\left(\widetilde{u}\left(\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v}),\,Y(t)\right)+\lambda Y(t)D(t)e^{\beta
t}H(t,\boldsymbol{v})\right)dt$ (18) $\displaystyle+$
$\displaystyle\\!\\!\\!\\!e^{-\beta\tau}\widetilde{U}\left(\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})\right)\Big{]},$

then

${\cal
J}(W_{0},Y_{0};c,l,\boldsymbol{\theta},\tau)\leq\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v})+\lambda
W_{0}.$ (19)

###### Proof.

Considering (13) and (14), it holds

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\cal
J}(W_{0},Y_{0};c,l,\boldsymbol{\theta},\tau)\\!=\\!E\Big{[}\\!\int_{0}^{\tau}\\!\\!\\!\\!e^{-\beta
t}\\!\left\\{u(c(t),l(t))\\!-\\!\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})(c(t)\\!+\\!Y(t)l(t))\right\\}dt$
$\displaystyle+e^{-\beta\tau}\left\\{U(W(\tau))-\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})W(\tau)\right\\}\Big{]}$
$\displaystyle+\lambda\,E\left[\int_{0}^{\tau}D(t)H(t,\boldsymbol{v})(c(t)+Y(t)l(t))dt+D(\tau)H(\tau,\boldsymbol{v})W(\tau)\right]$
$\displaystyle\leq E\left[\int_{0}^{\tau}e^{-\beta
t}\,\widetilde{u}\left(\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v}),\,Y(t)\right)dt+e^{-\beta\tau}\widetilde{U}\left(\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})\right)\right]$
$\displaystyle+\lambda\,E\left[\int_{0}^{\tau}D(t)H(t,\boldsymbol{v})(c(t)+Y(t)l(t))dt+D(\tau)H(\tau,\boldsymbol{v})W(\tau)\right].$
(20)

Exploiting the constraints (10), (11) and being $D(t)$ non increasing such
that $D(0)=1$ we obtain

$\displaystyle\\!\\!\\!\\!E\\!\left[\int_{0}^{\tau}\\!\\!D(t)H(t,\boldsymbol{v})(c(t)\\!+\\!Y(t)l(t))dt\\!+\\!D(\tau)H(\tau,\boldsymbol{v})W(\tau)\\!\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=E\\!\left[\int_{0}^{\tau}\\!\\!D(t)H(t,\boldsymbol{v})(c(t)\\!+\\!(l(t)\\!-\\!1)Y(t))dt\\!+\\!D(\tau)H(\tau,\boldsymbol{v})W(\tau)\\!+\\!\int_{0}^{\tau}\\!\\!Y(t)D(t)H(t,\boldsymbol{v})dt\\!\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=E\\!\left[\int_{0}^{\tau}\\!\\!Y(t)D(t)H(t,\boldsymbol{v})dt\\!\right]+E\\!\left[\\!H(\tau,\boldsymbol{v})W(\tau)\\!+\\!\int_{0}^{\tau}\\!\\!H(t,\boldsymbol{v})(c(t)\\!+\\!(l(t)\\!-\\!1)Y(t))dt\\!\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+E\\!\left[\int_{0}^{\tau}\\!\\!H(t,\boldsymbol{v})E_{t}\\!\left[\\!\frac{H(\tau,\boldsymbol{v})}{H(t,\boldsymbol{v})}W(\tau)\\!+\\!\int_{t}^{\tau}\\!\\!\frac{H(s,\boldsymbol{v})}{H(t,\boldsymbol{v})}(c(s)\\!+\\!(l(s)\\!-\\!1)Y(s))ds\right]dD(t)\right]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\leq
E\\!\left[\int_{0}^{\tau}\\!\\!Y(t)D(t)H(t,\boldsymbol{v})dt\\!\right]\\!+\\!W_{0}.$
(21)

Thus inequality (19) holds true. ∎

For a fixed stopping time $\tau\in{\cal S}$ and initial wealth $W_{0}$, we
denote by ${{\cal A}_{\tau}}$ the set of admissible processes
$(c,\,l,\,\boldsymbol{\theta})$. Thus we define

$V_{\tau}(W_{0},Y_{0}):=\sup_{(c,\,l,\,\boldsymbol{\theta})\in{\cal
A}_{\tau}}{\cal J}(W_{0},Y_{0};c,l,\boldsymbol{\theta},\tau).$ (22)

From Proposition 3.2, we obtain

$V_{\tau}(W_{0},Y_{0})\leq\inf_{\lambda>0,\,D(t)>0,\,\boldsymbol{v}\in
K(\sigma)}\left[\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v})+\lambda
W_{0}\right].$ (23)

and the following result holds true.

######  Proposition 3.3.

Equality in (23) holds if all the following conditions are satisfied:
\- for $0\leq t<\tau$

$\displaystyle\frac{\partial u}{\partial c}(c(t),l(t))=\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})\ \ \textrm{ and}\ \ \frac{\partial u}{\partial
l}(c(t),l(t))=\lambda Y(t)D(t)e^{\beta t}H(t,\boldsymbol{v})$ (24)

if $l(t)\leq L$, otherwise

$\displaystyle\frac{\partial u}{\partial c}(c(t),L)=\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})\ \ \textrm{ and}\ \ l(t)=L;$ (25)

\- it holds

$\displaystyle W(\tau)=I(\lambda D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})),$
$\displaystyle
E\left[\int_{0}^{\tau}H(s,\boldsymbol{v})(c(s)+(l(s)-1)Y(s))ds+H(\tau,\boldsymbol{v})W(\tau)\right]=W_{0};$

\- for any $t\in[0,\tau)$ such that $D(t)$ is not constant, i.e., $dD(t)\neq
0$,

$\displaystyle
E_{t}\left[\int_{t}^{\tau}\frac{H(s,\boldsymbol{v})}{H(t,\boldsymbol{v})}(c(s)+(l(s)-1)Y(s))ds+\frac{H(\tau,\boldsymbol{v})}{H(t,\boldsymbol{v})}W(\tau)\right]=0.$
(26)

###### Proof.

Equality in (23) holds if and only if equalities in (20) and (21) hold true.
Let us start with (20): this equality holds if

$\displaystyle u(c(t),l(t))-\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})(c(t)+Y(t)l(t))$ $\displaystyle=$
$\displaystyle\widetilde{u}\left(\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v}),\,Y(t)\right)$ $\displaystyle U(W(\tau))-\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})W(\tau)$ $\displaystyle=$
$\displaystyle\widetilde{U}\left(\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v})\right).$

By Proposition 3.1, this is equivalent to set for $0\leq t<\tau$

$\displaystyle\frac{\partial u}{\partial c}(c(t),l(t))=\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})\ \ \textrm{and}\ \ \frac{\partial u}{\partial
l}(c(t),l(t))=\lambda Y(t)D(t)e^{\beta t}H(t,\boldsymbol{v})$

if $l(t)\leq L$, otherwise

$\displaystyle\frac{\partial u}{\partial c}(c(t),L)=\lambda D(t)e^{\beta
t}H(t,\boldsymbol{v})\ \ \textrm{and}\ \ l(t)=L;$

and, reasoning as above for the utility function $U$, to impose

$\displaystyle U^{\prime}(W(\tau))=\lambda
D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v}),$

i.e., $W(\tau)=I(\lambda D(\tau)e^{\beta\tau}H(\tau,\boldsymbol{v}))$.
Now we deal with inequality (21): it becomes an equality if

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!E\\!\left[H(\tau,\boldsymbol{v})W(\tau)+\int_{0}^{\tau}H(t,\boldsymbol{v})(c(t)+(l(t)-1)Y(t))dt\right]=W_{0},$
(27)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!E\\!\left[\int_{0}^{\tau}\\!\\!\\!\\!H(t,\boldsymbol{v})E_{t}\\!\\!\left[\frac{H(\tau,\boldsymbol{v})}{H(t,\boldsymbol{v})}W(\tau)\\!+\\!\int_{t}^{\tau}\\!\\!\frac{H(s,\boldsymbol{v})}{H(t,\boldsymbol{v})}(c(s)+(l(s)-1)Y(s))ds\right]dD(t)\right]=0,$
(28)

and the last condition is equivalent to (26).
Notice that (27) and (28) imply (10) and (11), respectively. See [Choi et al.
(2008)] for further details. ∎

Optimality conditions in (24-25) provide a relation between the optimal
consumption, the leisure process and the triple ($\lambda,D(t),\,v$). As
described in [He and Pearson (1991), Section 3.2], we deal with the dual
Arrow-Debreu state price problem, which corresponds to minimize over
$\boldsymbol{v}\in K(\sigma)$ the maximum expected utility that the agent
could obtain in a complete market with state price density
$\Theta(\boldsymbol{v})$, defining the minimax local martingle measure. Thus
$\boldsymbol{v}^{*}$ is chosen in such a way to mimimize the agent wealth.
Once $\boldsymbol{v}^{*}$, the optimal multiplier $\lambda^{*}$ and the non-
increasing process $D^{*}(t)$ are computed, we can use these conditions to
obtain the optimal processes $c^{*}(t)$ and $l^{*}(t)$, as discussed in
[Barucci and Marazzina (2011), Appendix].

In order to compute the triple $(\lambda^{*},D^{*}(t),\boldsymbol{v}^{*})$, we
proceed as follows. Let

$\displaystyle\widetilde{V}(\lambda,Y_{0})$ $\displaystyle:=$
$\displaystyle\sup_{\tau\in[0,+\infty]}\inf_{\boldsymbol{v}\in
K(\sigma),\,D(t)>0}\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v})$
$\displaystyle=$
$\displaystyle\sup_{\tau\in[0,+\infty]}\inf_{D(t)>0}\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v}^{*}).$

Since (9) and (22) imply

${\cal V}(W_{0},Y_{0})=\sup_{\tau\in[0,+\infty]}V_{\tau}(W_{0},Y_{0}),$

then, under the conditions of Proposition 3.3, it holds

$\displaystyle{\cal V}(W_{0},Y_{0})$ $\displaystyle=$
$\displaystyle\sup_{\tau\in[0,+\infty]}\inf_{\lambda>0,\,D(t)>0}\left[\widetilde{V}(\lambda,D,\tau,Y_{0},\boldsymbol{v}^{*})+\lambda
W_{0}\right]$ $\displaystyle=$
$\displaystyle\inf_{\lambda>0}\left[\widetilde{V}(\lambda,Y_{0})+\lambda
W_{0}\right].$

Notice that this formulation is related to [He and Pearson (1991), Problem
(D’) page 277].

Let us define the process

$z(t,\boldsymbol{v})=\lambda D(t)e^{\beta t}H(t,\boldsymbol{v}),\ \
z(0,\boldsymbol{v})=\lambda$ (30)

and

$\displaystyle\phi(t,z,y):=\sup_{\tau>t}\inf_{\boldsymbol{v}\in
K(\sigma),\,D(t)>0}\\!\\!E\\!\Big{[}\\!\int_{t}^{\tau}\\!\\!\\!\\!e^{-\beta
s}\\!\left\\{\widetilde{u}(z(s,\boldsymbol{v}),Y(s))\\!+\\!Y(s)z(s,\boldsymbol{v})\right\\}ds\\!+\\!e^{-\beta\tau}\widetilde{U}(z(\tau,\boldsymbol{v}))$
$\displaystyle\Big{|}z(t,\boldsymbol{v})=z,Y(t)=y\Big{]}.$ (31)

Since $\phi(0,\lambda,Y_{0})=\widetilde{V}(\lambda,Y_{0})$, once $\phi$ is
computed, we can use (3) to obtain the value function ${\cal V}(W_{0},Y_{0})$
and thus the optimal couple $(\lambda^{*},D^{*}(t))$ and the optimal
strategies, as shown in [Barucci and Marazzina (2011), Appendix].

##  4 The Bellman equation and the free boundary problem

To address the optimization problem and to find out $\phi(\cdot)$ and the
function $\boldsymbol{v}^{*}$, we start from the dual function
$\widetilde{U}$. After retirement the optimal consumption-investment problem
is a classical Merton problem, considering the utility function (8), we have

$U(W(\tau))=\sup_{(c,\boldsymbol{\theta})}E\left[\int_{\tau}^{+\infty}e^{-\beta(t-\tau)}\frac{c^{\alpha(1-\gamma)}}{\alpha(1-\gamma)}dt\right],$

subject to the dynamic budget constraint (5) without labor income, i.e.,
$l(t)=1\,\forall t\geq\tau$.

It follows from standard results (see [Karatzas and Shreve (1998), Chapter 3])
that

$U(w)=\left(\frac{1}{\xi}\right)^{\Gamma}\frac{w^{1-\Gamma}}{1-\Gamma}$

where $\Gamma=1-\alpha(1-\gamma)$ and
$\xi=\frac{\Gamma-1}{\Gamma}\left(r+\frac{\boldsymbol{\Theta}^{2}(\boldsymbol{v}^{*})}{2\Gamma}\right)+\frac{\beta}{\Gamma}$.
In order to guarantee that the function $U$ is well defined, we assume
$\xi>0$. Thus, since $U$ is known analytically, we can also compute
$\widetilde{U}$ in (14):

$\widetilde{U}(z)=\frac{\Gamma}{\xi(1-\Gamma)}z^{\frac{\Gamma-1}{\Gamma}}.$

Since by Itô’s formula we have, for any $\boldsymbol{v}\in K(\sigma)$

$\displaystyle\frac{dz(t,\boldsymbol{v})}{z(t,\boldsymbol{v})}=\frac{dD(t)}{D(t)}+(\beta-r)dt-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})d\boldsymbol{B}(t),\
\ \ z(0)=\lambda,$

in order to derive the Bellman equation, we follow [He and Pages (1993)]
controlling the decreasing process $D$. Assume that $D$ is absolutely
continuous with respect to $t$, i.e., there exists $\psi\geq 0$ such that
$dD(t)=-D(t)\psi(t)dt$. If $\phi$ is twice differentiable with respect to $z$
and $y$, then the Bellman equation becomes

$\displaystyle\min_{\boldsymbol{v}\in K(\sigma),\,\psi\geq 0}\left\\{e^{-\beta
t}(\widetilde{u}\left(z,\,y\right)+zy)-\psi(t)z\frac{\partial\phi}{\partial
z}+\mathcal{L}_{t}^{0}\phi\right\\}=0,$ (32)

where

${\cal L}_{t}^{0}\phi=\frac{\partial\phi}{\partial
t}+(\beta-r)z\frac{\partial\phi}{\partial
z}+\mu_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}\boldsymbol{\Theta}^{2}(\boldsymbol{v})z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}+\frac{1}{2}\boldsymbol{\mu}_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{\mu}_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y}$

is the differential operator for the dual indirect utility function. Condition
(32) can be rewritten compactly as

$\displaystyle\min\left\\{e^{-\beta
t}(\widetilde{u}\left(z,\,y\right)+zy)+\mathcal{L}_{t}\phi,\
-\frac{\partial\phi}{\partial z}\right\\}=0,$ (33)

where

${\cal L}_{t}\phi=\frac{\partial\phi}{\partial
t}+\mu_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}\boldsymbol{\mu}_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}+(\beta-r)z\frac{\partial\phi}{\partial z}+\min_{\boldsymbol{v}\in
K(\sigma)}\left[\frac{1}{2}\boldsymbol{\Theta}^{2}(\boldsymbol{v})z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{\mu}_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y}\right]$

is the differential operator for the dual indirect utility function, i.e.,

${\cal L}_{t}\phi=\frac{\partial\phi}{\partial
t}+\mu_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}\boldsymbol{\mu}_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}+(\beta-r)z\frac{\partial\phi}{\partial
z}+\frac{1}{2}\boldsymbol{\Theta}^{2}(\boldsymbol{v}^{*})z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\boldsymbol{\Theta}^{\top}(\boldsymbol{v}^{*})\boldsymbol{\mu}_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y}$ (34)

and

$\boldsymbol{v}^{*}=arg\min_{\boldsymbol{v}\in
K(\sigma)}\frac{1}{2}\boldsymbol{\Theta}^{2}(\boldsymbol{v})z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{\mu}_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y}.$

Notice that the minimum over $\boldsymbol{v}\in K(\sigma)$ comes from the
minimax local martingale measure [He and Pearson (1991), Page 275]. Moreover,
since $\Theta(\boldsymbol{v})$ is a monotonic function with respect to
$\boldsymbol{v}$, it holds

${\cal L}_{t}\phi=\frac{\partial\phi}{\partial
t}+\mu_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}\boldsymbol{\mu}_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}+(\beta-r)z\frac{\partial\phi}{\partial
z}+\min_{\Theta(\boldsymbol{v}),\,\boldsymbol{v}\in
K(\sigma)}\left[\frac{1}{2}\boldsymbol{\Theta}^{2}(\boldsymbol{v})z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\boldsymbol{\Theta}^{\top}(\boldsymbol{v})\boldsymbol{\mu}_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y}\right]$

is the differential operator for the dual indirect utility function.
Furthermore, we recall that the space $K(\sigma)$ consists of vectors
$[0;x],\,x\in\mathbb{R}$, thus it holds

$[\boldsymbol{\Theta}]_{1}(\boldsymbol{v}^{*})=\frac{b-r}{\sigma}\ \ \ \
[\boldsymbol{\Theta}]_{2}(\boldsymbol{v}^{*})=\frac{[\boldsymbol{\mu}_{2}]_{2}y\frac{\partial^{2}\phi}{\partial
z\partial y}}{z\frac{\partial^{2}\phi}{\partial
z^{2}}}=[\boldsymbol{v}^{*}]_{2},$ (35)

where we denote with $[\cdot]_{i}$ the $i$-th component of a vector.

###  4.1 Correlated Wage

If the wage is given by

$dY(t)=M_{1}Y(t)dt+\rho M_{2}Y(t)dB_{1}(t)+M_{2}\sqrt{1-\rho^{2}}dB_{2}(t),$

i.e., $\boldsymbol{\mu_{2}}$ is equal to the vector $[\rho
M_{2},\,\sqrt{1-\rho^{2}}M_{2}]$, $\rho$ being the correlation between the
wage and the risky asset, then the above pde becomes

${\cal L}_{t}\phi=\frac{\partial\phi}{\partial
t}+M_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}M_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}+(\beta-r)z\frac{\partial\phi}{\partial
z}+\frac{1}{2}\Theta^{2}z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\rho\Theta M_{2}zy\frac{\partial^{2}\phi}{\partial z\partial
y}+\frac{1}{2}(\Theta^{*})^{2}z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\Theta^{*}\sqrt{1-\rho^{2}}M_{2}zy\frac{\partial^{2}\phi}{\partial
z\partial y},$

where

$\Theta=\sigma^{-1}(b-r),\ \ \textrm{and}\ \
\Theta^{*}=\frac{\sqrt{1-\rho^{2}}M_{2}y\frac{\partial^{2}\phi}{\partial
z\partial y}}{z\frac{\partial^{2}\phi}{\partial z^{2}}},$

i.e.,

${\cal L}_{t}\phi=\frac{\partial\phi}{\partial
t}+M_{1}y\frac{\partial\phi}{\partial
y}+\frac{1}{2}M_{2}^{2}y^{2}\frac{\partial^{2}\phi}{\partial
y^{2}}+(\beta-r)z\frac{\partial\phi}{\partial
z}+\frac{1}{2}\Theta^{2}z^{2}\frac{\partial^{2}\phi}{\partial
z^{2}}-\rho\Theta M_{2}zy\frac{\partial^{2}\phi}{\partial z\partial
y}-\frac{1}{2}(1-\rho^{2})M_{2}^{2}y^{2}\frac{\left(\frac{\partial^{2}\phi}{\partial
z\partial y}\right)^{2}}{\frac{\partial^{2}\phi}{\partial z^{2}}}.$

Notice that if $\rho=\pm 1$ we move back to the complete market case, and, if
$\rho=1$, we get the pde in [Barucci and Marazzina (2011), Section 5].

## References

  * [Barucci and Marazzina (2011)] Barucci, Emilio, and Marazzina, Daniele (2011) Optimal investment, stochastic labor income and retirement. Applied Mathematics and Computation 218-9, 5588-5604.
  * [Choi et al. (2008)] Choi, Kyoung Jin, Gyoocheol Shim and Yong Hyun Shin (2008) Optimal portfolio, consumption-leisure and retirement choice problem with CES utility. _Mathematical Finance_ 18, 445-472.
  * [Ding and Marazzina (2022)] Ding, Guodong, and Daniele Marazzina (2022) The impact of liquidity constraints and cashflows on the optimal retirement problem. _Finance Research Letters_ 49, 103159.
  * [Ding and Marazzina (2024)] Ding, Guodong, and Daniele Marazzina (2024) Effect of labour income on the optimal bankruptcy problem. _Annals of Operations Research_ 336, 773-795.
  * [Dybvig and Liu (2010)] Dybvig, Philip H., and Hong Liu (2010) Lifetime consumption and investment: retirement and constrained borrowing. _Journal of Economic Theory_ 145-3, 885-907.
  * [He and Pages (1993)] He, Hua, and Henri F. Pages (1993) Labor income, borrowing constraints, and equilibrium asset prices. _Economic Theory_ 3, 663-696.
  * [He and Pearson (1991)] He, Hua, and Neil D. Pearson (1991) Consumption and portfolio policies with incomplete market and short-sale constraints: the infinite dimensional case. _Journal of Economic Theory_ 54, 259-304.
  * [Karatzas and Shreve (1998)] Karatzas, Ioannis, and Steven E. Shreve (1998) _Methods of Mathematical Finance_. Series on Stochastic Modelling and Applied Probability, Vol. 39, Springer-Verlag.
  * [Tepla (2000)] Tepla, Lucie (2000) Optimal hedging and valuation of nontraded assets. European Finance Review 4, 231-251.


# Strain-gradient-induced switching of nanoscale domains in free-standing
ultrathin films

G. D. Belletti, S. D. Dalosto, Silvia Tinte  Instituto de Física del Litoral-
CONICET, Universidad Nacional del Litoral, Güemes 3450, (3000) Santa Fe,
Argentina

###### Abstract

We report first-principle atomistic simulations on the effect of local strain
gradients on the nanoscale domain morphology of free-standing PbTiO3 ultrathin
films. First, the ferroelectric properties of free films at the atomic level
are reviewed. For the explored thicknesses (10 to 23 unit cells), we find
flux-closure domain structures whose morphology is thickness dependent. A
critical value of 20 unit cells is observed: thinner films show structures
with 90∘ domain loops, whereas thicker ones develop, in addition, 180∘ domain
walls, giving rise to structures of the Landau-Lifshitz type. When a local and
compressive strain gradient at the top surface is imposed, the gradient is
able to switch the polarization of the downward domains, but not to the
opposite ones. The evolution of the domain pattern as a function of the strain
gradient strength consequently depends on the film thickness. Our simulations
indicate that in thinner films, first the $90^{o}$ domain loops migrate
towards the strain-gradient region, and then the polarization in that zone is
gradually switched. In thicker films, instead, the switching in the strain-
gradient region is progressive, not involving domain-wall motion, which is
attributed to less mobile $180^{o}$ domain walls. The ferroelectric switching
is understood based on the knowledge of the local atomic properties, and the
results confirm that mechanical flexoelectricity provides a means to control
the nanodomain pattern in ferroelectric systems.

###### pacs:

77.55.-g.,77.80.bn,77.80.Dj,77.90.+k

##  I Introduction

Flexoelectricity, the linear coupling between electric polarization $P$ and
inhomogeneous deformation, is always allowed by symmetry in all insulators
since strain gradients break inversion symmetry. This effect, found 50 years
ago kog64 , has traditionally received little attention owing to its
relatively weak effect in macroscopic samples. At the nanoscale, however,
strain gradients are typically much larger than at macroscopic scales, as
reduced dimensions imply larger gradients. The impressive progress in the
design and control of nanoscale structures has motivated a revival interest in
flexoelectricity, which has become an attractive topic in material science in
the past few years zub13 ; yud13 ; lee12 . In the context of perovskite thin
films, recent experimental breakthroughs leePRL11 ; gruv03 ; luSCI12  have
revealed that the effect of strain gradients can be enormous in thin films,
large enough to rotate or even switch ferroelectric domains. For instance, a
recent experimental work luSCI12 , which in part motivated ours, shows that an
inhomogeneous deformation caused by pressing the tip of an atomic force
microscope over the surface of BaTiO3 films is able to switch the electric
$P$.

To model in a more realistic manner the flexoelectric switching in ultrathin
ferroelectrics, it is important to include the domain morphology, which is
expected to depend strongly on the mechanical and electrical boundary
conditions yud13 . In the absence of electrodes or other sources of free
charges that compensate the depolarizing field, it is well known that
ferroelectrics respond by forming periodic arrays of $180^{o}$ domains, whose
character might be described using the traditional models of Kittel (open $P$
flow) or Landau-Lifshitz (closed $P$ flow). For instance, ordered $180^{o}$
stripe domains of the Kittel type have been observed in PbTiO3 (PTO) films as
thin as 3 unit cells deposited on SrTiO3 substrates, using x-ray photoelectron
diffraction fong04 .

From the theoretical side, to understand and predict properties of nanoscale
films at the atomistic level, there have been numerous studies using
$abinitio$ methods as well as simpler models constructed from them, such as
model Hamiltonians and atomistic core-shell models. In particular, in PTO
ultrathin films the ferroelectric and antiferrodistortive distortions in the
(001) surface have been described mey02 ; bung05 ; sepli0506 , indicating also
a decrease of tetragonality with stabilization of in-plane $P$ and
antiferrodistortive reconstruction in the region near the surface. Studies of
domain morphology reported vortex-like structures in films of a few unit cells
thick lai07 ; behe08 . The changes in $P$ and lattice deformation across the
$90^{o}$ and $180^{o}$ domain walls (DWs) hace also been investigated, with
calculations performed in the material bulk mey02 ; shim08  and in films of
3-unit-cell thickness shim10 ; qu13 . However, none of these works describe
the local properties of PTO ultrathin films and their correlation with the
domain morphology as the thickness varies. Concerning this kind of simulations
treating flexoelectricity, we can cite recent simulations involving a modified
model Hamiltonian technique for (BaxSr1-x)TiO3 ponoPRB12  which has been used
to calculate flexoelectric coefficients as a function of the film thickness
and temperature, but without making any reference to the domain presence.

Then, although the flexoelectricity effect on nanoscale is attracting much
attention, its exact role in the poling process of ultrathin films with a
nanodomain configuration has not been completely characterized from the
atomistic level. In this work, we address some particular questions, such as
How does a local strain gradient affect the domain pattern in free-standing
ultrathin films? and How does this depend on the film thickness and the strain
gradient strength? To answer these, using a first-principles-based atomistic
model we study the flexoelectric switching in ultrathin films of the
stereotypical ferroelectric PTO, when a longitudinal strain gradient is
applied locally over the center of one domain. Our simulations indicate that
for the explored films, thinner than 10 nm, the domain patterns in the
unstrained films are thickness dependent and responsible for the film being
locally deformed. In agreement with this, the film response to a local and
longitudinal strain gradient varies according to the nature of the domain
pattern and as a function of the gradient intensity. The paper is organized as
follows. Section II describes the computational details for modelling the film
and the conditions for mimicking a strain gradient. Section III presents the
results in two parts. First, the local properties of strain-free films are
reported for two thicknesses in order to understand later the gradient-induced
changes. Second, the domain pattern evolution for each of these films is
investigated as the strain gradient strength is enhanced. Finally, the paper
concludes with a summary in Sec.IV.

##  II Methodology: computational details

Prototype ferroelectric PTO is modeled using an atomistic core-shell model
fitted completely to first-principles calculations (i.e., no explicit
experimental data were used as input) as described in detail in ref. sepli04,
. The model contains fourth-order core-shell couplings, long-range Coulombic
interactions, and short-range interactions. It has shown to correctly
reproduce the cubic-tetragonal phase transition in bulk sepli04  and also was
able to describe the surface properties and interface effects of ultrathin
films sepli0506 .

Here we study free-standing ultrathin films of PTO under ideal open-circuit
electric boundary conditions performing classical molecular dynamics (MD)
simulations using the DL-POLY package dlpoly  at a constant temperature of 50
K, where PTO presents a tetragonal ferroelectric phase. Our films have two
PbO-terminated surfaces and are modeled using supercells of $N_{x}\times
N_{y}\times N_{z}$ unit cells. Periodic boundary conditions are applied along
the in-plane pseudocubic [100] and [010] directions (axes $x$ and $y$,
respectively) but not along the [001] one (axis $z$) defined as the growth
direction. The long-ranged electrostatic energy and forces are calculated by a
direct sum method wolf  which is a computationally efficient alternative to
the Ewald summation. In the $xy$ plane, the simulated supercell is kept
constant and constrained to form a simple square lattice ($N_{x}=N_{y}$) of
edge $N_{x}\times a$ with $a$=0.3866 nm, the in-plane lattice parameter in
bulk computed at 50 K with the present model. The film thicknesses range from
4.0 nm (10$c$) to 9.2 nm (23$c$), where $c$ is the lattice parameter in the
growth direction. Atoms are free to relax along all directions, except when
mentioned, as long as they remain inside the simulation box.

Local polarization is defined here as the dipole moment per unit volume of a
perovskite unit cell centered on the Ti ion and delimited by the eight Pb
nearest neighbors at the corners of the box. We consider contributions from
all atoms in the conventional cell and atomic positions with respect to that
center:

${\bf p}=\frac{1}{\it v}\sum_{\it i}{\frac{1}{\omega_{\it i}}{\it q_{i}}(\bf
r_{\it i}-\bf r_{\it Ti})}$ (1)

where $v$ is the volume of the cell, $q_{i}$ and $\bf r_{\it i}$ denote the
charge and the position of the $i$ particle, respectively, and $\omega_{i}$ is
a weight factor equal to the number of cells to which the particle belongs.
The reference position $\bf r_{\it Ti}$ corresponds to the core of the Ti
atom, and the sum runs over 29 particles, including cores and shells of the
surrounding ions (8 Pb and 6 O atoms) plus the shell of the Ti atom. Note that
$\bf p$ is independent of the origin for the $\bf r_{\it i}$ vectors and
vanishes when atoms are at their high-symmetry positions. The charge values
$q_{i}$ are given in Ref. sepli04, .

Figure 1: Optimal width $w$ of the simulation supercell for films of different
thickness $t$. For each film, $w$ and $t$ are informed by the pairs ($N_{x}a$,
$N_{z}c$) that label the supercell size along $x$ and $z$, respectively.
Inset: Schematic of the initial configuration of domains used to start the MD
simulations for all the thicknesses. Supercell boundaries are indicated by
dashed lines.

The initial configuration of the strain-free film for all the considered
thicknesses is set up by creating 180∘ stripe domains with polarizations
pointing upward and downward normal to the surface, alternated along the $x$
axis, as schematically shown in the inset of Fig. 1. The simulation supercell,
indicated by dashed lines, holds only two domains. Based on previous results
poy99 ; mey02  and also ours belle , which indicate that the 180∘ DWs centered
at the Pb-O lattice planes are more stable than those at the Ti-O2 ones, we
use the former in our starting configuration. Note that letting the structure
relax will allow this domain configuration to change and the optimized one
will depend on the film thickness as described below.

Figure 2: (Color online) Cell-by-cell polarization map projected on the $xz$
plane for strain-free films of 10$c$ (top panels) and 23$c$ (bottom panels)
thickness. Colors represent (left) the polarization $z$ component ($\mu$C/cm2)
and (right) the $c/a$ tetragonal relation. Domains are delimited by solid
lines: thick and thin for 90∘ and 180∘ domain walls, respectively. Large
arrows label the UP ($+P_{z}$) and DOWN ($-P_{z}$) domains. In the bottom-
right panel, the two types of 90∘ domain walls are indicated.

Strain gradient description. As we are interested in the effect of switching
the out-of-plane polarization, we apply strain gradients in the direction
normal to the surface, that is, a gradient of the out-of-plane strain
component $\eta_{zz}=c/c_{0}-1$ along the [001] direction denoted as
$\partial\eta_{zz}/\partial z$. We impose a linear gradient localized on an
in-plane square region of $4\times 4$ unit cells that penetrates until the
central plane, which compresses the upper plane by a given strain
$-\eta_{zz}^{max}$ and then decreases linearly in the intermediate planes
until vanishing at the central cells. One advantage of using the present
atomistic model to treat flexoelectric effects flex  is that the interatomic
potential does not need to be modified; it is sufficient to place the system
under appropriate mechanical-boundary conditions that mimic the required
gradient. To do that, one first has to establish the strain gradient in the
film, noting that its presence implies a stress gradient which produces net
forces acting on the atoms inside each unit cell. During MD simulations, the
atoms will try to relax in order to vanish these net forces. Thus, in order to
preserve the gradient a force pattern must be applied to the atoms in the unit
cell; i.e., certain atoms must be fixed along the full simulation.
Specifically, to set up the gradient, we start from the relaxed free film
(which already presents a domain structure) and displace rigidly the atomic
planes of the chosen unit cells until the average distance between consecutive
Pb-O lattice planes equals the required local strains $\eta_{zz}$ in order to
have $\partial\eta_{zz}/\partial z$. We choose to freeze out the $z$
components of the Pb-core positions, allowing all other internal degrees of
freedom (the Pb shell and other atoms, cores plus shells) to relax along all
directions. Other choices are possible, for instance, the force pattern could
be distributed between all atoms in the unit cell in different proportions
DHV_11 . Finally, note that the local strain gradient induces an internal
electric field which can compensate the depolarization field allowing that
part of the free-standing film is poled.

##  III Results

###  III.1 Free-standing ultrathin films of PbTiO3

We focus first on strain-free films under ideal open circuit boundary
conditions and revisit their local features. To begin with, the optimum
lateral size of the simulation supercell (i.e., the number of cells $N_{x}$)
is determined. For each considered thickness, we build up supercells of
different lateral sizes, relax their structures, compare their total energies,
and choose the lowest-energy one. The optimized structures are described
below; here we just present the optimum width $\omega$ of one domain period as
a function of the film thickness $t$. Figure 1 displays the linear relation
between $\omega$ and the root square of $t$, showing the fulfillment of the
Kittel law. We found that the optimum lateral sizes of the supercells are
26$a$, 30$a$, 34$a$, and 36$a$ for thickness of 10$c$, 16$c$, 20$c$, and
23$c$, respectively.

Hereafter we present results only for films of two thicknesses — 4.0 nm and
9.2 nm — because of their representative domain structures. Figure 2 shows a
side view of cell-by-cell polarization map for the two thicknesses, where
colors indicate the $z$-component values in the left panels and the $c/a$ cell
tetragonality in the right ones. Starting with the left panels, the 10$c$
ultrathin film (top panel) shows a flux-closure arrangement around central
cores or vortices. This vortex state, where the local $P$ orientation rotates
continuously around one geometrical core, can also be described as a
90∘-domain loop. Indeed, we can demarcate upward and downward $P$ domains
(hereafter UP and DOWN, respectively) with a rhombus shape surrounded by
triangular caps (where $P$ is almost parallel to the surface) that complete an
overall closed loop. Our simulations show that this domain morphology holds
for thicker films until thicknesses of $\sim$ 20 unit cells. Above that
critical value, the central vortices enlarge in the normal direction, giving
rise to 180∘ DWs that grow together with the film thickness. Consequently, the
shape of the UP and DOWN domains changes to hexagonal and a domain structure
as envisioned by Landau-Lifshitz is formed, as shown in the bottom panels for
a 23$c$-thickness film. In this structure, in particular, we see that the 180∘
DWs end up in vertices or confluences of three DWs as predicted by Srolovitz
and Scott scott86  and recently confirmed with high-resolution transmission
electron microscopy images of Pb(Zr1-xTix)O3 jia11 .

As regards the thicknesses of the DWs, they vary depending on the type and
location inside the film, the difference being more evident in the thicker
films. In the inner region, the orientation of the $P$ changes abruptly from
up to down over a distance of two lattice constants, creating very narrow 180∘
DWs, in agreement with previous theoretical results for defect-free walls
mey02 . The 90∘ DWs instead differ according to their position with respect to
the out-of-plane $P$ orientation, as indicated schematically in the bottom-
right panel in Fig. 2 for the UP domain as an example, but the same behavior
is repeated in the DOWN domain. In the top part, the $P$ rotates 90∘ in only
three unit cells, leading to narrow 90∘ DWs. In the bottom part, instead, the
$P$ rotates more continuously, requiring a few unit cells to change from the
in-plane $P$ at the surface until pointing perpendicular, making these 90∘ DWs
much broader. This behavior can be attributed to the electrostatic repulsion
between the head-to-head wu06 ; gur11  in-plane $P$s that confront each other
at the junction of two triangular domains.

The cell-by-cell tetragonal $c/a$ ratio is indicated by colors in the right
panels of Fig. 2, where $c$ and $a$ were obtained by taking the distances
between pairs of nearest-neighbor Pb atoms along $z$ and $x$, respectively.
The cells inside the triangular domains with $P$ parallel to the surfaces are
almost cubic, whereas those at the cores of the domains with $P$ perpendicular
are tetragonal and their values are thickness dependent. Further, and in
agreement with the high strain-$P$ interaction in PTO, these tetragonality
maps demonstrate the difference between the 90∘ DWs. The narrow upper ones are
well defined, presenting an abrupt change in the $c/a$ ratio as the tetragonal
axis rotates from perpendicular to parallel to the surface; however, the
broader ones at the bottom are characterized by a large tetragonality
(reaching, for instance, $\sim$8% in the 23$c$ film) which attenuates
gradually in a few unit cells towards both sides. Concerning the 180∘ DWs in
the thicker films, the tetragonality does not change notably, rather the wall
is defined by an abrupt relative displacement or offset of the Pb atoms along
$z$. Finally, we emphasize the significant lattice deformation observed at the
surfaces (shown by different colors in the right panels of Fig. 2) as well as
across the film thickness specially in the 90∘ DWs.

Figure 3: (Color online) Local lattice parameters $a$ and $c$ and modulus of
the polarization $z$-component across the thickness for the 10$c$ and 23$c$
films. Local parameters are averaged at the central cells of each domain —UP
[filled (blue) diamonds] and DOWN [filled (red) triangles]—and over the
complete $xy$ planes [open (green) circles].

To obtain further insights, we monitor those properties in detail through the
film thickness. Figure 3 shows the in-plane $a$ and out-of-plane $c$ lattice
parameters across the thickness for both films. Averaging these values over
complete planes parallel to the surface [open (green) circles], we find that
the average local $a$’s remain constant across the film due to the mechanical
boundary conditions imposed on the simulation supercell. The average $c$
values show the typical symmetric profile determined by the two free surfaces:
they are compressed at both surfaces and gradually increase towards the inner
film, consistent with previous results sepli0506 . However, examining closely
at the central region of each domain, the local $c$’s present an asymmetric
distribution: they are highly compressed in the two unit cells of the top
surface (with respect to the $P$ orientation) and slightly enlarged at the
bottom ones. The strain change along [001] is compensated in the plane only by
the local lattice parameters $a$, which behave oppositely to $c$, but not by
the $b$ parameters, which remain almost constants (not shown here),
corroborating that the unit cells at the surface neighborhood are highly
distorted. This particular asymmetry through the thickness was also predicted
in superlattices of PTO/SrTiO3 from first-principles calculations junq12  and
even in the limit case of PTO films of 1$c$ thickness whose unit cell is
distorted qu13 . Thus, this information together with the above $c/a$
description shows that the free surfaces induce spontaneous strain gradients
that break the inversion symmetry in the film, in agreement with previous
works shim10 ; qu13 . Furthermore, knowledge of the local strain maps is
essential to predict and understand the effect of an imposed strain gradient.
Moving to the right panels of Fig. 3, the profiles of the local $P_{z}$
modulus show that the averaged values per lattice plane are symmetric.
Remarkably, inside each domain, opposite to the highly asymmetric profiles of
local $c$, the $P_{z}$ profiles are almost symmetric (in thinner films,
actually they are slightly enhanced towards the top part). Then we note that
the modulus of $P$ and not $P_{z}$ follows the $c$ behavior across the film.

Table 1: Oxygen octahedra rotations (in degree) of the $AFD_{z}$ distortions for two thicknesses ($t$) in different atomic layers starting at the top surface, given at the center of the UP, DOWN and triangular domains, and averaging per full TiO2 atomic layer. $Abinitio$ results after Ref. bung05,  are reported in parenthesis. The last row of data corresponds to our calculations performed in the bulk material with two domains separated by 180∘ domain walls.  $t$ | Atomic layer | UP | DOWN | Triangular | Average
---|---|---|---|---|---
10$c$ | 1 | 14.2 | 13.3 | 13.0 | 13.0 (11.4)
| 2 | 1.0 | 2.9 | 3.5 | 2.8 (2.9)
| 3 | 6.0 | 6.2 | 6.4 | 6.5 (3.9)
| 4 | 4.9 | 5.1 | 5.2 | 4.9
| 5 | 4.8 | 5.0 | 5.8 | 5.2
23$c$ | 1 | 14.3 | 13.0 | 13.0 | 13.1 (11.4)
| 2 | 0.9 | 2.8 | 4.1 | 3.0 (2.9)
| 3 | 5.7 | 6.3 | 6.5 | 6.6 (3.9)
| 4 | 3.1 | 4.3 | 6.4 | 5.16
| 5 | 3.0 | 4.1 | 6.1 | 5.1
| 11 | 0.8 | 1.0 | 4.0 | 1.6
$bulk$ |  | 1.0 | 1.0 | 3.9 | 1.0 (3.3)

To end up this section, we analyze the local distribution of the
antiferrodistortive distortion which has been shown bung05 ; step01  to
coexist with the polar distortion at the surface, eventhough it essentially
vanishes in bulk. It involves the oxygen octahedron rotation around the
central Ti ion, where the rotation axis is the $z$ axis ($AFD_{z}$) with out-
of-phase movements between neighbor cells on the (001) plane but successive
in-phase movements along the [001]. Its strength is described by the rotation
angle of TiO4 in the normal layer at each unit cell. Table 1 lists the
$AFD_{z}$ rotation angle dependence of the atomic layers starting at the top
surface, for the two films of interest and in different regions of the film by
entering perpendicular to the surface: in the central cells of the UP and DOWN
domain, and crossing the center of one triangular domain (reaching the central
vortex of the 90∘ DW loop or the 180∘ DWs, depending on the film thickness).
The last column reports the values averaged for each complete TiO2 atomic
layer parallel to the surface. Moving along [001] in the center of one domain,
it is found that the rotation angle of $\sim 13^{o}$ in the superficial cells
decays towards the internal atomic layers: although it almost vanishes in the
second layer, recovers in the third one, and then keeps decreasing slowly. It
reaches the bulk value of $\sim$ 1o just at the eleventh layer of the 23$c$
film, while in the thinner film, the superficial effects predominate.
Furthermore, we note that cells in opposite surfaces differ in $\sim 1^{o}$,
which is also related to the different natures of the junction of the 90∘ DWs
as explained before. On the other hand, when entering the film through the
triangular domains, the $AFD_{z}$ distortions from the second layer towards
the inner ones are a bit larger than those in the center of the UP and DOWN
domains. Specifically, in the 180∘ DWs the $AFD_{z}$ distortion is about
$4^{o}$ (the octahedra rotate to allow the abrupt $P$ change, in addition to
the offset), versus the $\sim 1^{o}$ in the domain cores, where the cells are
tetragonal with ferroelectric distortions mostly along $z$. The averaged
values show a behavior similar to that in the cores of the out-of-plane $P$
domain but with values $\sim$ 1o larger than in the inner cells due to the
values in the DWs. Finally, note that our results for the first three atomic
layers can be compared, and indeed fully agree with previous calculations
reported using the present model sepli0506  as well as first-principles
methods bung05 . Although those calculations have provided evidence of
$AFD_{z}$ variation in the near-surface region, they were performed using
smaller supercells than here and without any reference to the effect of the
domain pattern. Hence, our inspection of the $AFD_{z}$ local distribution
gives insights into the dependence of the domain morphology, especially in the
DW neighborhood.

###  III.2 Film response to a local strain gradient

Figure 4: (Color online) Results for a 10$c$-thickness film. (a) Local lattice
parameter $c$ in the central region of the DOWN domain for a strain-free film
(triangles) and after the strain gradient is applied (circles). (b-d) Cell-by-
cell $P$ distributions for an unstrained film (b) and when a strain gradient
of -1.53 $\times~{}10^{7}$ m-1 (c) and -2.18$\times~{}10^{7}$ m-1 (d) imposed
in the zone delimited by a parallelepiped. For all cases, the $3D$ map and its
projection on the central $xy$ plane are shown. In (b) and (c) only one half
of the supercell projection on that plane is plotted, for clarity, and the
lines mark the position of the central vortices. In (d), the line in the $3D$
map also indicates the positions of the central vortices, and the projections
on three planes are presented: (d,i) on the central $xy$ plane, where dashed
lines indicate the positions of two lateral planes, $L_{x}$ (d,ii) and $L_{y}$
(d,iii). Colors denote the $z$ components of the polarization.

Now we investigate how the PTO film responds to a compressive strain gradient
applied locally as its intensity increases. We are particularly interested in
following the evolution of the domain configuration until eventually the
polarization switches in the strain-gradient zone. As the first step, we have
to choose where to apply the strain gradient. Based on the observed local
strain distribution in the central region of each domain (Fig. 3), we suppose
that a compressive strain gradient on the top surface of the film would only
accentuate the already existent inhomogeneous compression in the UP domain; on
the contrary, it would alter the almost-constant $c$ distribution in the DOWN
domain. Therefore, we apply to the latter a linear gradient on the top half of
the film along the $z$ axis as described in Sec. II and shown in Fig. 4(a)
where the local $c$ profile is plotted before and after the application for a
particular $\eta_{zz}^{max}$ of -4.25 %. The same profile is established on a
region of $4\times 4$ unit cells that extends from the top surface to the
central plane grad , as schematically indicated in the figure by a
parallelepiped. The strain gradient strength is increased progressively by
varying the compressive strain in the superficial cells $\eta_{zz}^{max}$ in
intervals of -0.25 %. For each strain gradient value, we carry out an MD
simulations until the system rearranges to that strain condition.

Thinner films. The results for the 10$c$ film are illustrated in Figs. 4(b)–
4(d) using three-dimensional (3D) maps of the local $P$ and projections of
them on different planes. The unstrained film is shown in Fig. 4(b), and how
it changes when the imposed local strain gradient $g$ takes two particular
values, namely, $g_{1}=-1.53\times 10^{7}$ m-1 and $g_{2}=-2.18\times 10^{7}$
m-1 in Figs. 4(c) and 4(d), respectively. Starting with the unstrained film,
which is repeated here so as to have a reference (same as the top panel of
Fig. 2), the 3D map displays a domain structure of 90∘ domain loops. Its
projection on the central $xy$ plane shows only two domains with the $P$
downward and upward perpendicular to the paper plane separated by the central
vortices, whose positions aligned along the $y$ axis are indicated by a line
in the figure. The square in the DOWN domain identifies the region where the
strain gradient will be applied. Thus, when the local strain gradient is
imposed and as its strength $g$ increases, we find that the response of the
film presents three different behaviors.
(I) For $g<g_{1}$, no appreciable changes are observed and the local $P$
distribution remains similar to that of the unstrained film.
(II) For $g_{1}\leq g<g_{2}$, when $g=g_{1}$, corresponding to an
$\eta_{zz}^{max}$ of -3.00%, the whole domain structure is displaced, almost
without changing either the size or the shape of the domains as shown in the
3D map of Fig. 4(c). The projection on the $xy$ plane clearly shows that the
central vortices have been shifted in $-x$ towards the gradient region and sit
just before its border. $P$ distributions similar to this one are observed as
the gradient is further increased.
(III) For $g=g_{2}$ or $\eta_{zz}^{max}=$ -4.25%, the dipoles rearrange and
the local $P$s in the gradient region, initially downward, are switched [3D
map of Fig. 4(d)] inver . This region has merged with the original UP domain,
extending its size and therefore changing its shape. The new domain
distribution can be appreciated more easily in the projection on the central
$xy$ plane [Fig. 4(d,i)]. Also, we examine two lateral planes that cross the
gradient region. The projection on a $yz$ plane shown in Fig. 4(d,ii) (where
all cells pointed downward before the gradient was applied) reveals that the
new domain pattern has developed a flux-closure arrangement around the central
vortices analogous to that found in the unstrained film. The projection on the
$L_{y}$ plane [Fig. 4(d,iii)] shows that most cells, including those in the
gradient region, are poled upward. The asymmetric distribution indicates that
the switching process has started in the vortices that were brought close to
the gradient region and now have disappeared, while the other vortices remain
on the right side.
Note that due to the supercell in-plane periodicity, it is equally likely that
in (II) the vortices had shifted in $+x$ towards the gradient region and
therefore the switching process would have started at the other side. Further,
we note that the $c$ profile in the gradient region, whose inferior part [open
circles in Fig. 4(a)] has been allowed to relax, has adopted a profile quite
similar to the one observed in the UP domains of the unstrained films, as was
assumed initially. Finally, we want to remark that the observed movements of
domains in our free-standing films are allowed because of the absence of a
substrate. Its presence, however, can hamper the movement of domains, as
shown, for instance, in piezoelectric force microscopy images of
Pb(Zr0.2Ti0.8)O3 patterns in epitaxial thin films naga02 .

Figure 5: (Color online) Surface of the hypertoroidal moment $z$ component
$h_{z}$ (in units of electronic charge $e$) as a function of the central $xy$
plane of the simulation cell, for a 10$c$-thickness film (a) unstrained and
(c) under a local strain gradient $g_{2}$. For clarity, we include, to the
right of each surface, the corresponding $3D$ polarization pattern. (b) Curves
on a transveral $xz$ plane (indicated by dashed lines) that crosses the strain
gradient region.

Hypertoroidal moment. In order to better describe the evolution of these
complex dipolar configurations, we use the recently defined order parameter
$\bf h$ named hypertoroidal moment prosan07 , which provides a measure of
subtle microscopic features. It involves the double cross product of the local
dipoles $\bf p_{\it i}$ with the vectors ${\bf r}_{\it i}$ locating their
positions,

${\bf h}=\frac{1}{4V}\sum_{i}{\bf r}_{\it i}\times({\bf r}_{\it i}\times{\bf
p}_{\it i_{t}})_{\it t}$ (2)

where $V$ denotes the supercell volume, the subindex $t$ labels the
transversal component, and the summation runs over all unit cells within the
supercell. The transversal component can be approximated as ${\bf p}_{\it
i_{t}}\approx{\bf p}_{\it i}-\langle{\bf p}\rangle$, $\langle{\bf p}\rangle$
being the average of the individual dipoles over all the sites. This
approximation is also valid for the expression in parenthesis. We refer the
reader to Ref. prosan07,  for a detailed discussion of this parameter. For 2D
systems (i.e., systems periodic in two Cartesian directions and finite in the
third one) such as the films considered here, once the simulation supercell is
chosen, the $\bf h$ value is independent of the origin used for $\bf r_{\it
i}$. However, the magnitude and even the sign of $\bf h$ do depend on the
choice made for the supercell, leading to the multivaluedness of the $\bf h$
parameter. We compute the hypertoroidal moment for different supercells whose
centers scan the complete central $xy$ plane of the original supercell.
Therefore we explore for the first time a surface of $\bf h$ and find that the
changes in value and/or sign reflect modifications in the full $P$
distribution.

Although $P$ switching is a dynamical and complex process, we may gain some
insight by examining the surface of only the $z$-component of the
hypertoroidal moment $h_{z}$ (the other two Cartesian components are almost
null) as the imposed strain gradient increases. Beginning with the unstrained
film shown in Fig. 5(a), we corroborate that the $h_{z}$ surface displays a
sinusoidal dependence along the $x$ direction with the same period as the
nanodomains and remains constant along the $y$ direction prosan07 .
Specifically, the maxima (minima) of $h_{z}$ extend along $y$, coinciding with
the positions of the central cores of the UP (DOWN) domain, while $h_{z}$
vanishes when rotation around the perpendicular $y$ axis exists, certainly
where the vortex centers sit. Then as the strain-gradient intensity increases
and following the $P$ distributions, the $h_{z}$ response also presents three
behaviors separated by the gradient values $g_{1}$ and $g_{2}$.
(I) For $g<g_{1}$, the sinusoidal shape remains almost unchanged.
(II) For $g_{1}\leq g<g_{2}$, when $g=g_{1}$, in agreement with the $P$
distribution, the $h_{z}$ surface shifts in $-x$ as shows by its projection on
a transverse cut that crosses the central gradient region [see the central
curve in Fig. 5(b)]. $h_{z}$ behaves similarly as the gradient increases.
(III) For $g=g_{2}$, needed to fully invert the $P$ in the affected region,
the $h_{z}$ surface shown in Fig. 5(c) has changed dramatically, indicating
that the 90∘-domain loop configuration has been modified. Two global extrema
are observed, one maximum and one minimum, revealing that the periodicity of
the domain structure and thus the extended cores of the each domain have
vanished. From the $\bf h$ definition, we interpret that the $h_{z}$ minimum
represents the part of the DOWN domain farthest away from the gradient region
and therefore the less affected one, whereas the $h_{z}$ maximum is located in
the region of the UP domain farthest away from the distribution of $-P_{z}$.
Fig. 5(b) shows that the $h_{z}$ profile in this particular direction has
slightly gained height in the positive branch, versus an important reduction
in the negative one, reflecting the size enhancement of the UP domain at
expenses of the DOWN domain.
Summarizing, the evolution of the hypertoroidal moment as a function of the
strain gradient strength provides information of the domain-structure
evolution localizing the domain cores and the vortices in a regular structure,
as well as their redistribution when the $P$s in the strain region are finally
switched.

Figure 6: (Color online) Two-dimensional projection of the polarization
patterns for an $xz$ plane that crosses the gradient region in a
23$c$-thickness film subject to a strain gradient of (a)
$g_{1}^{*}=-1.01\times 10^{7}$ m-1 and (b) $g_{2}^{*}=-1.57\times 10^{7}$ m-1,
just before and after the switching, respectively. The gradient region and the
immobile 180∘ DW are indicated.

Thicker films. Figure 6 displays the $P$-map projection on a transverse cut
crossing the gradient region of a $23c$-thickness film for two strain-gradient
values of $g_{1}^{*}=-1.01\times 10^{7}$ m-1 and $g_{2}^{*}=-1.57\times
10^{7}$ m-1, corresponding to $\eta_{zz}^{max}$ of -5.00 and -7.50 %,
respectively. As the gradient increases, starting at the surface the local
$P$s near and inside the gradient region rotate progressively until pointing
parallel to the surface, i.e., their $P_{z}$ components decrease gradually [as
shown in Fig. 6(a) for $g_{1}^{*}$] until they fully vanish. A lateral view in
a $yz$ plane, not shown here, indicates that the zone with $P$ parallel to the
surface starts to grow and expands towards the bottom of the film as the
gradient increases, then they are switched. For the $g_{2}^{*}$ value [Fig.
6(b)] the affected region is already poled upward and the asymmetric
distribution shows that one of the 180∘ DWs is still immobile. Remarkably, as
the gradient is increased further, the switched zone moves forward, generating
vortices at the bottom as it continues growing laterally. Thus, unlike the DW
movements observed in thinner films, we find that the $180^{o}$ DWs never move
forward to the gradient region. This fact can be attributed to the relatively
high energy barrier for the coherent motion of one entire $180^{o}$ DW
compared to the $90^{o}$ DW case.mey02

In summary, we investigate from an atomistic point of view flexoelectricity-
induced domain evolution. In free-standing PbTiO3 ultrathin films with
thickness-dependent domain patters, we study the ferroelectric switching
induced by a strain gradient acting locally in the direction normal to the
surface. We first review the properties of strain-free films and confirm that
the inhomogeneity of polarization, strain (strain gradient), and
antiferrodistortive rotational angles appears not only among the surface cells
but also in the neighborhood of the DWs. For the explored thicknesses (10 to
23 unit cells), we find flux-closure domain structures whose morphology is
thickness dependent. A critical value of 20 unit cells is found: thinner films
show structures with 90∘ domain loops, whereas thicker ones develop, in
addition, $180^{o}$ DWs giving rise to structures of the Landau-Lifshitz type.
When a compressive strain gradient normal to the surface acts locally at the
top surface of the film, we find that it is able to reverse the normal
polarization in DOWN domains but not in UP domains, which is understood based
on the knowledge of the local strain distribution before imposition of the
gradient. The evolution of the domain pattern as a function of the strain
gradient strength depends, accordingly, on the film thickness. In thinner
films, first one of the center of the $90^{o}$ DW loops moves towards the
gradient region and then the affected region is gradually switched. In thicker
films, instead, no $180^{o}$ DW movement is observed and the local
polarizations in the gradient region are progressively reversed as the
gradient strength is enhanced. This work contributes from an atomistic point
of view to the understanding of the mechanical flexoelectricity in tuning
domains in ultrathin films.

The authors acknowledge support from the Argentinian Agency ANPCyT through
Grants PICT-PRH 99 and PICT-PRH 195. G.D.B. thanks to CONICET-Argentina for
the Ph.D. fellowship.

## References

  * (1) S. M. Kogan, Sov. Phys. Solid State 5, 2069 (1964).
  * (2) P. Zubko, G. Catalan, and A.K. Tagantsev, Annu. Rev. Mater. Res. 43, 387 (2013).
  * (3) P. V Yudin and a K. Tagantsev, Nanotechnology 24, 432001 (2013).
  * (4) D. Lee, and T. W. Noh, Phil. Trans. R. Soc. A 370, 4944 (2012).
  * (5) D. Lee, A. Yoon, S. Y. Jang, J.-G. Yoon, J.-S. Chung, M. Kim, J. F. Scott, and T. W. Noh, Phys. Rev. Lett. 107, 057602 (2011).
  * (6) A. Gruverman, B.J. Rodriguez, a. I. Kingon, R.J. Nemanich, J.S. Cross, and M. Tsukada, Appl. Phys. Lett. 82, 3071 (2003).
  * (7) H. Lu, C.-W. Bark, D. Esque de los Ojos, J. Alcala, C. B. Eom, G. Catalan, A. Gruverman, Science 336, 59 (2012).
  * (8) Fong, G.B. Stephenson, S.K. Streiﬀer, J. A. Eastman, O. Auciello, P. H. Fuoss, and C. Thompson, Science 304, 1650 (2004).
  * (9) B. Meyer and D. Vanderbilt Phys. Rev. B 65, 104111 (2002)
  * (10) C. Bungaro and K. M. Rabe, Phys. Rev. B 71, 035420 (2005).
  * (11) M. Sepliarsky, M. G. Stachiotti, and R. L. Migoni, Phys. Rev. B 72, 014110 (2005), Phys. Rev. Lett. 96, 137603 (2006).
  * (12) B.-K. Lai, I. Ponomareva, I. Kornev, L. Bellaiche, and G. J. Salamo, Appl. Phys. Lett. 91, 152909 (2007).
  * (13) R.K. Behera, B.B. Hinojosa, S.B. Sinnott, A. Asthagiri, and S.R. Phillpot, J. Phys. Condens. Matter 20, 395004 (2008).
  * (14) T. Shimada, K. Wakahara, Y. Umeno, and T. Kitamura, J. Phys. Condens. Matter 20, 325225 (2008).
  * (15) T. Shimada, S. Tomoda, and T. Kitamura, Phys. Rev. B 81, 144116 (2010).
  * (16) B. Yin and S. Qu, J. Appl. Phys. 114, 063703 (2013).
  * (17) I. Ponomareva, A. K. Tagantsev, L. Bellaiche, Phys. Rev. B 85, 104101 (2012).
  * (18) M. Sepliarsky, Z. Wu, A. Asthagiri, and R. E. Cohen, Ferroelectrics 301, 55 (2004).
  * (19) W. Smith and T. R. Forester, computer code DL-POLY, Daresbury and Rutherford Appleton Laboratory, Daresbury, England.
  * (20) D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999).
  * (21) S. Pöykkö and D.J. Chadi, Appl. Phys. Lett. 75, 2830 (1999); J. Phys. and Chem. of Solids 61, 291 (2000).
  * (22) G. Belletti and S. Tinte, unpublished (2014).
  * (23) We have checked that the present atomistic model for PTO is able to reproduce the longitudinal flexoelectric coefficient in the range of $nC/m$. belle
  * (24) J. Hong and D. Vanderbilt, Phys. Rev. B 84, 180101 (2011).
  * (25) D. J. Srolovitz, and J. F. Scott, Phys. Rev. B 34, 1815 (1986).
  * (26) C.-L. Jia, K. W. Urban, M. Alexe, D. Hesse and I. Vrejoiu. Science 331, 1420 (2011).
  * (27) X. Wu and D. Vanderbilt, Phys. Rev. B 73, 020103 (2006).
  * (28) M. Y. Gureev, A. K. Tagantsev, and N. Setter, Phys. Rev. B 83, 184104 (2011).
  * (29) Pablo Aguado-Puente and Javier Junquera, Phys. Rev. B 85, 184105 (2012).
  * (30) A. Munkholm, S.K. Streiffer, M.V. Ramana Murty, J.A. Eastman, C. Thompson, O. Auciello, L. Thompson, J.F. Moore, and G.B. Stephenson, Phys. Rev. Lett. 88, 016101 (2001).
  * (31) Two comments on the arbitrary size of the affected region: (i) that its deep depends on the film thickness indicates that, for the same compressive strain on the superficial unit cell, the gradient is weaker in thicker films; (ii) that the lateral size is kept constant means that the affected region is farther away from the 180∘ DW as the domain size increases together with the film thickness.
  * (32) For a given value of strain gradient, the switching process ocurrs as a function of time during a molecular dynamics simulation, noting that the switching kinetics is beyond the scope of this work.
  * (33) V. Nagarajan, A. Roytburd, A. Stanishevsky, S. Prasertchoung, T. Zhao, et al. Nat. Mater. 2, 43–4 (2002).
  * (34) S. Prosandeev, L. Bellaiche, J. Mater. Sci. 44, 5235 (2009).