Hadi Daneshmand
email: seyed.daneshmand@inria.fr

Amir Joudaki
email: amir.joudaki@inf.ethz.ch

Francis Bach
email: francis.bach@inria.fr

Batch Normalization Orthogonalizes Representations in Deep Random Networks

This paper underlines a subtle property of batch-normalization (BN): Successive batch normalizations with random linear transformations make hidden representations increasingly orthogonal across layers of a deep neural network. We establish a non-asymptotic characterization of the interplay between depth, width, and the orthogonality of deep representations. More precisely, under a mild assumption, we prove that the deviation of the representations from orthogonality rapidly decays with depth up to a term inversely proportional to the network width. This result has two main implications: 1) Theoretically, as the depth grows, the distribution of the representation -after the linear layers-contracts to a Wasserstein-2 ball around an isotropic Gaussian distribution. Furthermore, the radius of this Wasserstein ball shrinks with the width of the network. 2) Practically, the orthogonality of the representations directly influences the performance of stochastic gradient descent (SGD). When representations are initially aligned, we observe SGD wastes many iterations to orthogonalize representations before the classification. Nevertheless, we experimentally show that starting optimization from orthogonal representations is sufficient to accelerate SGD, with no need for BN.

Introduction

Batch Normalization (BN) [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] enhances training across a wide range of deep network architectures and experimental setups [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] Huang | Densely connected convolutional networks[END_REF][START_REF] Silver | Mastering the game of go without human knowledge[END_REF]. The practical success of BN has inspired research into the underlying mechanism of BN [START_REF] Santurkar | How does batch normalization help optimization?(no, it is not about internal covariate shift)[END_REF][START_REF] Karakida | The normalization method for alleviating pathological sharpness in wide neural networks[END_REF][START_REF] Arora | Theoretical analysis of auto rate-tuning by batch normalization[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF]. BN influences first-order optimization methods by avoiding the rank collapse in deep representation [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF], direction-length decoupling of optimization [START_REF] Kohler | Exponential convergence rates for batch normalization: The power of length-direction decoupling in non-convex optimization[END_REF], influencing the convergence of the steepest descent [START_REF] Arora | Theoretical analysis of auto rate-tuning by batch normalization[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF], and smoothing the optimization objective function [START_REF] Santurkar | How does batch normalization help optimization?(no, it is not about internal covariate shift)[END_REF][START_REF] Karakida | The normalization method for alleviating pathological sharpness in wide neural networks[END_REF]. However, the benefits of BN go beyond its critical role in optimization. For example, [START_REF] Frankle | Training batchnorm and only batchnorm: On the expressive power of random features in cnns[END_REF] shows that BN networks with random weights also achieve surprisingly high performance after only minor adjustments of their weights. This striking result motivates us to study the representational power of random networks with BN.

We study hidden representations across layers of a laboratory random BN with linear activations. Consider a batch of samples passing through consecutive BN and linear layers with Gaussian weights. The representations of these samples are perturbed by each random linear transformation, followed by a non-linear BN. At first glance, the deep representations appear unpredictable 35th Conference on Neural Information Processing Systems (NeurIPS 2021), Sydney, Australia.

after many stochastic and non-linear transformations. Yet, we show that these transformations orthogonalize the representations. To prove this statement, we introduce the notion of "orthogonality gap", defined in Section 3, to quantify the deviation of representations from a perfectly orthogonal representation. Then, we prove that the orthogonality gap decays exponentially with the network depth and stabilizes around a term inversely related to the network width. More precisely, we prove

E orthogonality gap = O (1 -α)
depth + batchsize α √ width holds for α > 0 that is an absolute constant under a mild assumption. In probability theoretic terms, we prove stochastic stability of the Markov chain of hidden representations [START_REF] Kushner | Stochastic stability and control[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF][START_REF] Khasminskii | Stochastic stability of differential equations[END_REF]. The orthogonality of deep representations allows us to prove that the distribution of the representations after linear layers contracts to a Wasserstein-2 ball around isotropic Gaussian distribution as the network depth grows. Moreover, the radius of the ball is inversely proportional to the network width. Omitting details, we prove the following bound holds:

Wasserstein 2 (representations, Gaussian) 2 = O (1 -α) depth (batchsize) + (batchsize) 2 α √ width .
The above equation shows how depth, width, and batch size, interact with the Gaussian approximation of the representations. Since the established rate is exponential with depth, the distribution of the representations stays in a Wasserstein ball around isotropic Gaussian distribution after a few layers. Thus, BN not only stabilizes the distribution of the representations, which is its main promise [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], but also enforces Gaussian isotropic distribution in deep layers.

There is growing interest in bridging the gap between neural networks, as the most successful parametric methods for learning, and Gaussian processes and kernel methods, as well-understood classical models for learning [START_REF] Jacot | Neural tangent kernel: Convergence and generalization in neural networks[END_REF][START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF][START_REF] Lee | Deep neural networks as gaussian processes[END_REF][START_REF] Bietti | On the inductive bias of neural tangent kernels[END_REF][START_REF] Huang | Kernel methods match deep neural networks on timit[END_REF]. This link is inspired by studying random neural networks in the asymptotic regime of infinite width. The seminal work by [START_REF] Neal | Bayesian learning for neural networks[END_REF] sparks that a single-layer network resembles a Gaussian process as its width goes to infinity. However, increasing the depth may significantly shift the distribution of the representations away from Gaussian [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. This distributional shift breaks the link between Gaussian processes and deep neural networks. To ensure Gaussian representations, [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF] suggests increasing the width of the network proportional to the network depth. Here, we show that BN ensures Gaussian representations even for deep networks with finite width. This result bridges the link between deep neural networks and Gaussian processes in the regime of finite width. Many studies rely on deep Gaussian representations in an infinite width setting [START_REF] Yang | A mean field theory of batch normalization[END_REF][START_REF] Schoenholz | Deep information propagation[END_REF][START_REF] Pennington | The emergence of spectral universality in deep networks[END_REF][START_REF] Klambauer | Self-normalizing neural networks[END_REF][START_REF] De Palma | Random deep neural networks are biased towards simple functions[END_REF]. Our non-asymptotic Gaussian approximation can be incorporated into their analysis to extend these results to the regime of finite width.

Since training starts from random networks, representations in these networks directly influence training. Hence, recent theoretical studies has investigated the interplay between initial hidden representations and training [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF][START_REF] Frankle | Training batchnorm and only batchnorm: On the expressive power of random features in cnns[END_REF][START_REF] Schoenholz | Deep information propagation[END_REF][START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF][START_REF] Bahri | Statistical mechanics of deep learning[END_REF]. In particular, it is known that hidden representations in random networks without BN become correlated as the network grows in depth, thereby drastically slowing training [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF][START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF][START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF]. On the contrary, we prove that deep representations in networks with BN are almost orthogonal. We experimentally validate that initial orthogonal representations can save training time that would otherwise be needed to orthogonalize them. By proposing a novel initialization scheme, we ensure the orthogonality of hidden representations. Such an initialization effectively avoids the training slowdown with depth for vanilla networks, with no need for BN. This observation further motivates studying the inner workings of BN to replace or improve it for deep learning. Theoretically, we made the following contributions:

1. For MLPs with batch normalization, linear activation, and Gaussian weights, we prove that representations across layers become increasingly orthogonal up to a constant inversely proportional to the network width. 2. Leveraging the orthogonality, we prove that the distribution of the representations contracts to a Wasserstein ball around a Gaussian distribution as the depth grows. Up to the best of our knowledge, this is the first non-asymptotic Gaussian approximation for deep neural networks with finite width.

Experimentally, we made the following contribution1 :

3. Inspired by our theoretical understanding, we propose a novel weight initialization for standard neural networks that ensure orthogonal representations without BN. Experimentally, we show that this initialization effectively avoids training slowdown with depth in the absence of BN.

Preliminaries

Notations

Akin to [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF], we focus on a Multi-Layer Perceptron (MLP) with batch normalization and linear activation. Theoretical studies of linear networks is a growing research area [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF][START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF][START_REF] Bartlett | Gradient descent with identity initialization efficiently learns positive-definite linear transformations by deep residual networks[END_REF]Arora et al., 2019a). When weights are initialized randomly, linear and non-linear networks share similar properties such as the rank collapse issue studied in [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF]. For ease of analysis, we assume activations are linear. Yet, we will argue that our findings extend to non-linear in Appendix G.

We use n to denote batch size, and d to denote the width across all layers, which we further assume is larger than the batch size d ≥ n. Let H ∈ R d×n denote representations for n samples in layer , with H 0 ∈ R d×n corresponding to n input samples in the batch with d features. Successive representations are connected by Gaussian weight matrices W ∼ N (0, I d /d) as

H +1 = 1 √ d BN (W H ), BN (M ) = diag(M M ) -1 /2 M, (1)
where diag(.) zeros out off-diagonal elements of its input matrix, and the scaling factor 1/ √ d ensures that all matrices {H k } k have unit Frobenius norm (see Appendix 2). The BN function in Eq. ( 1) differs slightly from the commonly used definition for BN as the mean correction is omitted. However, [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] shows this difference does not change the network properties qualitatively. Appendix G further discusses about the role of the mean correction term in our analysis. Readers can find all the notations in Appendix Table 2. [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] observe that if inputs are linearly independent, then their hidden representations remain linearly independent in all layers as long as d = Ω(n 2 ). Under technical assumptions, [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] establishes a lower-bound on the average of the rank of hidden representations over infinite layers. Based on this study, we assume that the linear independence holds and build our analysis upon that. This avoids restating the technical assumptions of [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] and also further technical refinements of their theorems. The next assumption presents the formal statement of the linear independence property. Assumption A 1 (α, ). There exists an absolute positive constant α such that the minimum singular value of H k is greater than (or equal to) α for all k = 1, . . . , .

The linear independence of hidden representations

The linear independence of the representations is a shared property across all layers. However, the representations constantly change when passing through random layers. In this paper, we mathematically characterize the dynamics of the representations across layers.

Orthogonality of deep representations 3.1 A warm-up observation

To illustrate the difference between BN and vanilla networks, we compare hidden representations of two input samples across the layers of these networks. Figure 1 plots the absolute value of cosine similarity of these samples across layers. This plot shows a stark contrast between vanilla and BN networks: While representations become increasingly orthogonal across layers of a BN network, they become increasingly aligned in a vanilla network. More specifically, we observe that BN is able to orthognalized almost aligned representations; while, the vanilla network provides almost align representation of two orthogonal samples in deep layers. While this behaviour has been theoretically studied for vanilla networks [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF][START_REF] Bjorck | Understanding batch normalization[END_REF][START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF], up to the best of our knowledge, it is not theoretically analyzed for BN networks for networks with finite width regime. In the following section, we formalize and prove this orthogonalizing property for BN networks with finite widths.

Theoretical analysis

The notion of orthogonality gap plays a central role in our analysis. Given the hidden representation H ∈ R d×n , matrix H H constitutes inner products between representations of different samples.

Note that H H ∈ R n×n is different form the covariance matrix HH /n ∈ R d×d . The orthogonality gap of H is defined as the deviations of H H from identity matrix, after proper scaling. More precisely, define V : R d×n \ 0 → R + as

V (H) := 1 H 2 F H H - 1 I n 2 F I n F .
(2)

The following theorem establishes a bound on the orthogonality of representation in layer . Theorem 1. Under Assumption A 1 (α, ), the following holds:

E [V (H +1 )] ≤ 2 1 - 2 3 α + 3n α √ d . (3)
Assumption A 1 is studied by [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] who note that A 1 (α, ∞) holds as long as d = Ω(n 2 ). If A 1 does not hold, one can still prove that there is a function of representations that decays with depth up to a constant (see Appendix C).

The above result implies that BN is an approximation algorithm for orthogonalizing the hidden representations. If we replace diag(M ) -1 by (M ) -1 in BN formula, in Eq. ( 1), then all the hidden representation will become exactly orthogonal. However, computing the inverse of a non-diagonal d × d matrix is computationally expensive, which must repeat for all layers throughout training, and differentiation must propagate back through this inversion. The diagonal approximation in BN significantly reduces the computational complexity of the matrix inversion. Since the orthogonality gap decays at an exponential rate with depth, the approximate orthogonality is met after only a few layers. Interestingly, this yields a desirable cost-accuracy trade-off: For a larger width, the orthogonality is more accurate, due to term 1/ √ d in the orthogonality gap, and also the computational gain is more significant.

From a different angle, Theorem 1 proves the stochastic stability of the Markov chain of hidden representations. In expectation, stochastic processes may obey an inherent Lyapunov-type of stability [START_REF] Kushner | Stochastic stability and control[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF][START_REF] Khasminskii | Stochastic stability of differential equations[END_REF]. One can analyze the mixing and hitting times of Markov chains based on the stochastic stability of Markov chains, [START_REF] Kemeny | Markov chains[END_REF][START_REF] Eberle | Markov processes[END_REF]. In our analysis, the orthogonality gap is a Lyapunov function characterizing the stability of the chain of hidden representations. This stability opens the door to more theoretical analysis of this chain, such as studying mixing and hitting times. For, example the stability can be used for mixing analysis of the chain {H 1 , . . . , H , . . . }. Let π denote the stationary distribution of this chain. Under a particular stochastic stability condition called geometric drift condition,

|E [ϕ(H )] -E H∼π [ϕ(H)]| ≤ α |E [ϕ(H 0 )] -E H∼π [ϕ(H)]|
holds for α ∈ (0, 1) and measurable function ϕ : R d×n → R (see Theorem 3.6 of [START_REF] Hairer | Convergence of markov processes[END_REF]). The drift condition holds if there exists a Lyapunov function L : R d×n → [0, 1] and constant K ≥ 0 such that the following holds:

E [L(H +1 )|H ] ≤ γL(H ) + K.
We prove that the above condition holds in Appendix Eq (12). However, such theoretical analyses may not be of interest to the machine learning community, hence we focus on the implications of these results for understanding the underpinnings of neural networks.

It is helpful to compare the orthogonality gap in BN networks to studies on vanilla networks [START_REF] Bjorck | Understanding batch normalization[END_REF][START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF][START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF]. The function implemented by a vanilla linear network is a linear transformation as the product of random weight matrices. The spectral properties of the product of i.i.d. random matrices are the subject of extensive studies in probability theory [START_REF] Bougerol | Products of random matrices with applications to Schrödinger operators[END_REF]. Invoking these results, one can check that the orthogonality gap of the hidden representations in vanilla networks rapidly increases since the rank of representations converges to a one matrix as the depth grows [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF]. We refer readers to Appendix D for more details about vanilla networks.

Remarkably, [START_REF] Agarwal | A deep conditioning treatment of neural networks[END_REF] proves if activation functions obey a self-normalization property, then a specific kernel of hidden representations becomes well-condition as the network depth grows. However, it is not clear how to impose the self-normalization property. Here, we establish the whitening for an explicitly defined normalization used in practice.

Experimental validations

Our experiments presented in Fig. 2a validate the exponential decay rate of V with depth. In this plot, we see that log(V ) linearly decreases for = 1, . . . , 20, then it wiggles around a small constant.

Our experiments in Fig. 2b suggest that the O(1/ √ d) dependency on width is almost tight. Since V (H ) rapidly converges to a ball, the average of V (H ) over layers estimates the radius of this ball. This plot shows how the average of V (H ), over 500 layers, changes with the network width, validating the O(1/ √ d) dependency implied by Theorem 1. Consider an input matrix H 0 ∈ R d×n containing n samples in R d . Assuming that elements of this matrix are i.i.d. zero-mean and unit variance random variables, the minimum singular value of H 0 is greater than

√ d - √ n.
In practical applications, the batchsize used for normalization is relatively smaller than the network width, hence A 1 (Θ(1), 0) holds. For the intermediate representations, we also observed that A 1 (α, ) holds -for an α independent from -as long as d n. Let α 0 denote the minimum singular value of H 0 . For different values of d and n, we check whether A 1 (0.1α 0 , 1000) holds. Fig. 3 illustrates that this assumption hold for d = Ω(n 2 ), which is also confirmed by [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF]. 4 Gaussian approximation

Orthogonality yields Gaussian approximation

The established result on the orthogonality gap allows us to show that the after linear layers are approximately Gaussian. If H is orthogonal, the random matrix W H is equal in distribution to a standard Gaussian matrix due to the invariance of standard Gaussian distribution under linear orthogonal transformations. We can formalize this notion of Gaussianity by bounding the Wasserstein-2 distance between the distribution of W H and standard Gaussian distribution. Let W 2 (R 1 , R 2 ) denote the Wasserstein-2 distance between probability distributions of two random variables R 1 , and R 2 . The next lemma formally establishes the link between the orthogonality and the distribution of the representations. Lemma 2. Given G ∈ R d×n with i.i.d. zero-mean 1/d-variance Gaussian elements, the following Gaussian approximation holds:

W 2 W H , G/ √ n 2 ≤ 2nE [V (H )] . (4)
Combining the above result with Theorem 1 yields the result presented in the next corollary. Corollary 3. For G ∈ R d×n with i.i.d. zero-mean 1/d-variance Gaussian elements,

W 2 W H , G/ √ n 2 ≤ 4n 1 - 2 3 α + 6n 2 α √ d (5)
holds under Assumption A 1 (α, ).

In other words, the distribution of the representations contracts to a Wasserstein 2 ball around an isotropic Gaussian distribution as the depth grows. The radius of the Wasserstein 2 ball is at most O(1/ √ width). As noted in the last section, A 1 is extensively studied by [START_REF] Daneshmand | Batch normalization provably avoids rank collapse for randomly initialised deep networks[END_REF] where it is shown that A 1 (α > 0, ∞) holds as long as d = Ω(n 2 ).

Deep neural networks as Gaussian processes

Leveraging BN, Corollary 3 establishes the first non-asymptotic Gaussian approximation for deep random neural networks. For vanilla networks, the Gaussianity is guaranteed only in the asymptotic regime of infinite width [START_REF] Garriga-Alonso | Deep convolutional networks as shallow gaussian processes[END_REF][START_REF] Neal | Bayesian learning for neural networks[END_REF][START_REF] Lee | Deep neural networks as gaussian processes[END_REF][START_REF] Neal | Bayesian learning for neural networks[END_REF][START_REF] Hazan | Steps toward deep kernel methods from infinite neural networks[END_REF]. Particularly, [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF] links vanilla networks to Gaussian processes when their width is infinite and grows in successive layers, while our Gaussian approximation holds for networks with finite width across layers. The link between Gaussian processes and infinite-width neural networks has inspired several studies to rely on Gaussian representations deep random networks [START_REF] Klambauer | Self-normalizing neural networks[END_REF][START_REF] De Palma | Random deep neural networks are biased towards simple functions[END_REF][START_REF] Pennington | The emergence of spectral universality in deep networks[END_REF][START_REF] Schoenholz | Deep information propagation[END_REF][START_REF] Yang | A mean field theory of batch normalization[END_REF]. Assuming the representations are Gaussian, [START_REF] Klambauer | Self-normalizing neural networks[END_REF] designed novel activation functions that improve the optimization performance, De Palma et al. ( 2019) studies the sensitivity of random networks, [START_REF] Pennington | The emergence of spectral universality in deep networks[END_REF] highlights the spectral universality in random networks, [START_REF] Schoenholz | Deep information propagation[END_REF] studies information propagating through the network layers, and [START_REF] Yang | A mean field theory of batch normalization[END_REF] studies gradients propagation through the depth. Indeed, our analysis implies that including BN imposes the Gaussian representations required for these analyses. Although our result is established for linear activations, we conjecture a similar result holds for non-linear MLPs (see Appendix G).

The orthogonality and optimization

In the preceding sections, we elaborated on the theoretical properties of BN networks in controlled settings. In this section, we demonstrate the practical applications of our findings. In the first part, we focus on the relationship between depth and orthogonality. Increasing depth drastically slows the training of neural networks with BN. Furthermore, we observe that as depth grows, the training slowdown highly correlates with the orthogonality gap. This observation suggests that SGD needs to orthogonalize deep representations in order to start classification. This intuition leads us to the following question: If orthogonalization is a prerequisite for training, can we save optimization time by starting from orthogonal representations? To test this experimentally, we devised a weight initialization that guarantees orthogonality of representations. Surprisingly, even in a network without BN, our experiments showed that this initialization avoids the training slow down, affirmatively answering the question.

Throughout the experiments, we use vanilla MLP (without BN) with a width of 800 across all hidden layers, ReLU activation, and used Xavier's method for weights intialization [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. We use SGD with stepsize 0.01 and batch size 500 and for training. The learning task is classification with cross entropy loss for CIFAR10 dataset (Krizhevsky et al., 2009, MIT license). We use PyTorch (Paszke et al., 2019, BSD license) and Google Colaboratory platform with a single Tesla-P100 GPU with 16GB memory in all the experiments. The reported orthogonality gap is the average of the orthogonality gap of representation in the last layer.

Orthogonality correlates with optimization performance

In the first experiment, we show that the orthogonality of representations at the initialization correlates with optimization speed. For networks with 15, 30, 45, 60, and 75 widths, we register training loss after 30 epochs and compare it with the initial orthogonality gap. Figure 4a shows the training loss (blue) and the initial orthogonality gap (red) as a function of depth. We observe that representations are more entangled, i.e., orthogonal, when we increase depth, coinciding with the training slowdown. Intuitively, the slowdown is due to the additional time SGD must spend to orthogonalize the representations before classification. In the second experiment, we validate this intuitive argument by tracking the orthogonality gap during training. Figure 4b plots the orthogonality gap of output and training loss for a network with 20 layers. We observe that SGD updates are iteratively orthogonalizing representations, marked by the reduction in the orthogonality gap.

Learning with initial orthogonal representations

We have seen that the slowdown in SGD for deeper networks correlates with the orthogonality gap before training. Here we show that by preemptively orthogonalizing representations, we avoid the slowdown with depth. While in MPL with linear activations, hidden representations remain orthogonal simply by taking orthogonal matrices as weights [START_REF] Pennington | The emergence of spectral universality in deep networks[END_REF][START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF], the same does not hold for networks with non-linear activations, such as ReLU. To enforce the orthogonality in the absence of BN, we introduce a dependency between weights of successive layers that ensures deep representations remain orthogonal. More specifically, we incorporate the SVD decomposition of the hidden representation of each layer into the initialization of the subsequent layer. To emphasis this dependency between layers and to distinguish it from purely orthogonal weight initialization, we refer to this as iterative orthogonalization.

We take a large batch of samples n ≥ d, as the input batch for initialization. Let us assume that weights are initialized up to layer W -1 . To initialize W , we compute SVD decomposition of the representations H = U Σ V where matrices U ∈ R d×d and V ∈ R n×d are orthogonal. Given this decomposition, we initialize W by

W = 1 Σ 1 /2 F V Σ -1 /2 U , (6)
where V ∈ R d×d is an orthogonal matrix obtained by slicing V ∈ R n×d . Notably, the inverse in the above formula exists when n is sufficiently larger than d2 . It is easy to check that V (W H ) < V (H ) holds for the above initialization (see Appendix F), similar to BN. By enforcing the orthogonality, this initialization significantly alleviates the slow down of training with depth (see Fig. 5), with no need for BN. This initialization is not limited to MLPs. In Appendix H, we compare iterative orthogonalization with a BN-replacement method recently proposed by [START_REF] Brock | High-performance large-scale image recognition without normalization[END_REF]. In Appendix I, we propose a similar SVD-based initialization for convolutional networks that effectively accelerates the training of deep convolutional networks.

Discussion

To recap, we proved the recurrence of random linear transformations and BN orthogonalizes samples. Our experiments underline practical applications of this theoretical finding: starting from orthogonal representations effectively avoids the training slowdown with depth for MLPs. Based on our experimental observations in Appendix I, we believe that this result extends to standard convolution networks used in practice. In other words, a proper initialization ensuring the orthogonality Although our theoretical bounds hold for MLPs with linear activations, our experiments confirm similar orthogonal stability for various neural architectures. Appendix Figure 6 demonstrates this stability for MLPs with ReLU activations and a convolutional networks. This plot compares the evolution of the orthogonality gap, through layers, for BN networks with vanilla networks. For vanilla networks, we observe that the gap increases with depth. On the contrary, the gap decreases by adding BN layers and stabilizes at a term that is constant with regard to depth. Based on our observations, we conjecture that hidden representations of modern BN networks obey similar stability.

A more formal statement of the conjecture is presented in Appendix G. [START_REF] Lubana | Beyond batchnorm: Towards a general understanding of normalization in deep learning[END_REF] experimentally compares properties of different normalization techniques including layer normalization (LN) [START_REF] Ba | Layer normalization[END_REF]. According to this study, LN does not necessarily orthogonalize the outputs of deep neural networks. Hence, more theoretical studies are required to understand the essence of different normalization techniques in deep learning.

Figure 1 :

 1 Figure 1: Orthogonality: BN vs. vanilla networks. The horizontal axis shows the number of layers, and the vertical axis shows the absolute value of cosine similarity between two samples across the layers (d = 32). Mean and 95% confidence intervals of 20 independent runs.

Figure 2 :

 2 Figure 2: Orthogonality gap vs. depth and width. Left: log(V (H )) vertically versus horizontally. Right: log( 1 500

Figure 3 :

 3 Figure 3: Validations for A1 Pixel (n, d) marks whether A1 holds in all 10 independent runs: The black color indicates A1(0.1α0, 1000) failed in at least once (where α0 is the minimum singular value of H0). The blue curve marks d = (n -28) 1.8 highlighting A1 holds for d = Ω(n 2 ).

Figure 4 :

 4 Figure 4: Orthogonality and Optimization Left: the orthogonality gap at initialization (red, left axis) and the training loss after 30 epochs (blue, right axis) with depth. Right: the orthogonality gap (red, left axis) and the training loss in each epoch (blue, right axis). Mean and 95% confidence interval of 4 independent runs.

Figure 5 :

 5 Figure 5: Iterative orthogonalization. Horizontal axis: depth. Vertical axis: the training loss after 30 epochs for Xavier's initialization (blue), our initialization (red). Mean and 95% confidence interval of 4 independent runs.

Table 1 :

 1 Table 1 briefly compares Gaussian approximation for vanilla and BN networks. BN MLP (Cor. 3) (in)finite (in)finite In a O(width -1 /4 )-W 2 ball around Gaussian Distribution of representations in random vanilla and BN networks. For the convolutional network, the width refers to the number of channels. Results for Vanilla MLPs and Vanilla convolution networks are establish by Matthews et al. (2018), and Garriga-Alonso et al. (2019), respectively. Remarkably, Corollary 3 holds for MLP with linear activations.

	Network	Width	Depth	Distribution of Outputs
	Vanilla MLP	infinite	finite	Converges to Gaussian as width → ∞
	Vanilla Convnet	infinite	finite	Converges to Gaussian as width → ∞

Implementations are available at https://github.com/hadidaneshmand/batchnorm21.git

We may inductively assume that H is almost orthogonal by the choice of W1, . . . , W -1 . Thus, Σ is invertible.

Acknowledgements

We thank Gideon Dresdner, Vincent Fortuin, and Ragnar Groot Koerkamp for their helpful comments and discussions. This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001(PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grant SEQUOIA 724063).


Anna N Morozovska

Eugene A Eliseev

Sergei V Kalinin
email: sergei2@ornl.gov

Riccardo Hertel
email: riccardo.hertel@ipcms.unistra.fr

Flexo-Sensitive Polarization Vortices in Thin Ferroelectric Films

We consider theoretically the influence of the flexoelectric coupling on the spatial distribution and temperature behavior of spontaneous polarization for several types of stable domain structure in thin ferroelectric films, such as stripe domains and vortices. Finite element modelling (FEM) for BaTiO3 films and analytical calculations within the Landau-Ginzburg-Devonshire approach reveals that an outof-plane polarization component can be very sensitive to the flexoelectric coupling for periodic quasi-2D stripe domains and 3D vortex-antivortex structures. However, the influence is rather different for these structures. The flexoelectric coupling increases significantly the amplitude of a small out-of-plane polarization component in the stripe domains, and the "up" or "down" direction of the component is defined by the sign of the flexoelectric coefficients. Concerning the vortex-antivortex pairs, their antivortices with in-plane anti-сirculation have smooth wide dipolar cores through the entire film, whose shape and other features are almost insensitive to the coupling. The vortices with in-plane vorticity have spike-like cores with an out-of-plane quadrupolar moment induced by the flexoelectric coupling. The cores are located near the film-dead layer interfaces. FEM results corroborated by analytical calculations prove that a change of the flexoelectric coefficient sign leads to a reorientation of the core axial polarization, making the flexo-sensitive 3D vortices similar to the recently introduced "flexons" in cylindrical nanoparticles. The relatively wide temperature range (from 200 to 400 K) of the flexosensitive vortices' existence gives us the hope that they can be observed experimentally in thin ferroelectric films by scanning probe and nonlinear optical microscopy methods.

I. Introduction

Since its appearance and until now, nanoscale ferroics (ferromagnets, ferroelectrics, ferroelastics) have been the main object of fundamental research on the physical nature of long-range order of polar, magnetic, and structural properties [START_REF] Glinchuk | Nanoferroics[END_REF][START_REF] Hlinka | Skyrmions in ferroelectric materials[END_REF][START_REF] Gruverman | Piezoresponse force microscopy and nanoferroic phenomena[END_REF]. The leading role is played by the emergence of a domain structure of long-range order parameters, such as electric polarization, magnetization, and (anti)ferrodistortion, and its interaction with the surface of a nanoferroic [START_REF] Kalinin | Surface-screening mechanisms in ferroelectric thin films and their effect on polarization dynamics and domain structures[END_REF][START_REF] Wang | Understanding, Predicting, and Designing Ferroelectric Domain Structures and Switching Guided by the Phase-Field Method[END_REF].

In the case of ferroelectrics, which we discuss here, the long-range order of the domains is governed by both electrostatic fields and the distribution of strain fields and their gradients. A strong and spatially extended gradient of elastic fields generated by the surface, domain walls, and/or other inhomogeneities creates a sufficiently large flexoelectric polarization, which is proportional to this gradient (this phenomenon is a direct flexoelectric effect) [START_REF] Tagantsev | Piezoelectricity and flexoelectricity in crystalline dielectrics[END_REF]. On the other hand, the electrical polarization gradient gives rise to an inhomogeneous strain (inverse flexoelectric effect) [7]. The thermodynamic description of the flexoelectric effect is given by Lifshitz invariants [START_REF] Tolédano | Landau Theory of Phase Transitions[END_REF][START_REF] Wang | Flexoelectricity in solids: Progress, challenges, and perspectives[END_REF].

The flexoelectric effect exists in all ferroics [START_REF] Tagantsev | Piezoelectricity and flexoelectricity in crystalline dielectrics[END_REF][7][START_REF] Tolédano | Landau Theory of Phase Transitions[END_REF][START_REF] Wang | Flexoelectricity in solids: Progress, challenges, and perspectives[END_REF], in contrast to the flexomagnetic [START_REF] Lukashev | Flexomagnetic effect in frustrated triangular magnetic structures[END_REF] and flexomagnetoelectric [START_REF] Pyatakov | Flexomagnetoelectric interaction in multiferroics[END_REF][START_REF] Tanygin | Symmetry theory of the flexomagnetoelectric effect in the magnetic domain walls[END_REF] effects, the existence of which is critically sensitive to the presence of time inversion and its connection with other operations of point symmetry of a particular material [START_REF] Eliseev | Linear magnetoelectric coupling and ferroelectricity induced by the flexomagnetic effect in ferroics[END_REF][START_REF] Eliseev | Effective flexoelectric and flexomagnetic response of ferroics[END_REF].

The influence of the flexoelectric effect on the macroscopic properties of ferroics (as well as any other materials) is relatively small [START_REF] Tagantsev | Piezoelectricity and flexoelectricity in crystalline dielectrics[END_REF][7][START_REF] Tolédano | Landau Theory of Phase Transitions[END_REF][START_REF] Wang | Flexoelectricity in solids: Progress, challenges, and perspectives[END_REF]. However, in nanosized and nanostructured ferroics, the influence of the flexoelectric effect on their polar, magnetic, electronic properties and phase transitions can be significant and may lead to fundamental changes in their properties [START_REF] Majdoub | Title[END_REF][START_REF] Eliseev | Spontaneous flexoelectric/flexomagnetic effect in nanoferroics[END_REF][START_REF] Ahmadpoor | Flexoelectricity in two-dimensional crystalline and biological membranes[END_REF]. This size-dependence of the effect arises from the gradients of physical quantities playing a leading role in such nanosystems [START_REF] Catalan | The effect of flexoelectricity on the dielectric properties of inhomogeneously strained ferroelectric thin films[END_REF][START_REF] Catalan | Strain gradients in epitaxial ferroelectrics[END_REF].

Recently, it was predicted theoretically that a decrease in the correlation-gradient polarization energy leads to a significant increase in the polarization gradient, which in turn leads to spontaneous bending of otherwise uncharged domain walls in thin multiferroic films [START_REF] Eliseev | Intrinsic structural instabilities of domain walls driven by gradient couplings: meandering anferrodistortive-ferroelectric domain walls in BiFeO3[END_REF] and ferroelectric nanoparticles [START_REF] Eliseev | Labyrinthine domains in ferroelectric nanoparticles: Manifestation of a gradient-induced morphological phase transition[END_REF]. Such domain walls can form meandering [START_REF] Eliseev | Intrinsic structural instabilities of domain walls driven by gradient couplings: meandering anferrodistortive-ferroelectric domain walls in BiFeO3[END_REF] and/or labyrinthine [START_REF] Eliseev | Labyrinthine domains in ferroelectric nanoparticles: Manifestation of a gradient-induced morphological phase transition[END_REF] structures.

Later, similar structures were discovered experimentally in thin BiFeO3 and Pb(Zr0.4Ti0.6)O3 films by HR-STEM [START_REF] Han | Mapping gradient-driven morphological phase transition at the conductive domain walls of strained multiferroic films[END_REF] and PFM [START_REF] Nahas | Inverse transition of labyrinthine domain patterns in ferroelectric thin films[END_REF] methods, respectively, and corroborated by ab initio calculations [START_REF] Nahas | Inverse transition of labyrinthine domain patterns in ferroelectric thin films[END_REF].

However, the influence of flexo-effects on the morphology of domain structures in thin films has not yet been studied systematically.

Most published experimental studies of flexoelectric phenomena in proper [START_REF] Ma | Flexoelectric effect in ceramic lead zirconate titanate[END_REF][START_REF] Ma | Flexoelectricity of barium titanate[END_REF][START_REF] Lee | Giant Flexoelectric Effect in Ferroelectric Epitaxial Thin Films[END_REF] and incipient [START_REF] Zubko | Strain-gradient-induced polarization in SrTiO3 single crystals[END_REF][START_REF] Zubko | Erratum: Strain-Gradient-Induced Polarization in SrTiO3 Single Crystals[END_REF] ferroelectrics, atomistic quantum-mechanical [START_REF] Bersuker | Pseudo Jahn-Teller effect in the origin of enhanced flexoelectricity[END_REF][START_REF] Maranganti | Atomistic determination of flexoelectric properties of crystalline dielectrics[END_REF] and first-principles calculations [START_REF] Ponomareva | Finite-temperature flexoelectricity in ferroelectric thin films from first principles[END_REF][START_REF] Hong | First-principles theory and calculation of flexoelectricity[END_REF][START_REF] Stengel | Flexoelectricity from density-functional perturbation theory[END_REF][START_REF] Dreyer | Current-density implementation for calculating flexoelectric coefficients[END_REF] are overwhelmingly aimed at determining the magnitude of the flexoelectric coefficients [START_REF] Zubko | Erratum: Strain-Gradient-Induced Polarization in SrTiO3 Single Crystals[END_REF][START_REF] Bersuker | Pseudo Jahn-Teller effect in the origin of enhanced flexoelectricity[END_REF][START_REF] Maranganti | Atomistic determination of flexoelectric properties of crystalline dielectrics[END_REF][START_REF] Ponomareva | Finite-temperature flexoelectricity in ferroelectric thin films from first principles[END_REF][START_REF] Hong | First-principles theory and calculation of flexoelectricity[END_REF][START_REF] Stengel | Flexoelectricity from density-functional perturbation theory[END_REF][START_REF] Dreyer | Current-density implementation for calculating flexoelectric coefficients[END_REF] and the structure of the flexoelectric tensor [START_REF] Shu | Symmetry of flexoelectric coefficients in crystalline medium[END_REF][START_REF] Quang | The number and types of all possible rotational symmetries for flexoelectric tensors[END_REF]. The majority of the phenomenological papers on this topic are devoted to the influence of the flexo-effects on the macroscopic properties and phase transitions of the ferroic film as a whole (see e.g. Refs. [START_REF] Morozovska | Flexocoupling impact on the size effects of piezo-response and conductance in mixed-type ferroelectrics-semiconductors under applied pressure[END_REF][START_REF] Vorotiahin | Tuning the Polar States of Ferroelectric Films via Surface Charges and Flexoelectricity[END_REF][START_REF] Eliseev | Defect driven flexo-chemical coupling in thin ferroelectric films[END_REF] and refs therein), and only a few of them address the actual influence of the flexoelectric effect on the structure of domain walls [START_REF] Morozovska | Interfacial Polarization and Pyroelectricity in Antiferrodistortive Structures Induced by a Flexoelectric Effect and Rotostriction[END_REF][START_REF] Morozovska | Anisotropic conductivity of uncharged domain walls in BiFeO3[END_REF][START_REF] Yudin | Bichiral structure of ferroelectric domain walls driven by flexoelectricity[END_REF][START_REF] Eliseev | Structural phase transitions and electronic phenomena at 180-degree domain walls in rhombohedral BaTiO3[END_REF][START_REF] Gu | Non-Ising Character of a Ferroelectric Wall Arises from a Flexoelectric Effect[END_REF]. The main result of the papers [40 -44] is the prediction that the flexoelectric effect induces the appearance of a small but sufficiently strong polarization component perpendicular to the plane of the nominally uncharged domain wall of the ferroic, resulting in an effective "charging" of the wall. PFM and cAFM experiments registering the conductivity of nominally uncharged domain walls in ferroelectrics and multiferroics confirm the theoretical predictions [START_REF] Maksymovych | Tunable metallic conductance in ferroelectric nanodomains[END_REF][START_REF] Balke | Enhanced electric conductivity at ferroelectric vortex cores in BiFeO3[END_REF][START_REF] Vasudevan | Domain wall geometry controls conduction in ferroelectrics[END_REF]. Non-Ising and chiral ferroelectric domain walls have been revealed by nonlinear optical microscopy [START_REF] Cherifi-Hertel | Non-ising and chiral ferroelectric domain walls revealed by nonlinear optical microscopy[END_REF].

Note that only a few works analyze the influence of the flexoelectric coupling on the curled vortex-like domain structures in ferroelectric thin films [START_REF] Li | Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling[END_REF] and nanoparticles [START_REF] Eliseev | Flexo-Elastic Control Factors of Domain Morphology in Core-Shell Ferroelectric Nanoparticles: Soft and Rigid Shells[END_REF][START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF]. Namely, using machine learning and phase-field modeling, Li et at. [START_REF] Li | Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling[END_REF] analyze the role of flexoelectricity on polar 2D vortices in a PbTiO3/SrTiO3 superlattice. The axis of these vortices is parallel to the film surface, and the flexoelectric coupling influences their details in a quantitative way, but not sharply or critically. The influence of flexoelectricity on polarization vortices in spherical [START_REF] Eliseev | Flexo-Elastic Control Factors of Domain Morphology in Core-Shell Ferroelectric Nanoparticles: Soft and Rigid Shells[END_REF] and cylindrical [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF] core-shell ferroelectric nanoparticles can be significant because it can induce a small axial polarization of the vortex core. Analytical calculations and simulations [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF] have proven that a change of the flexoelectric coefficient sign leads to a reorientation of the axial polarization of the vortex and that an anisotropy of the flexoelectric coupling critically influences the vortex core formation and its related domain morphology. Here we consider the influence of the flexoelectric coupling on the spatial distribution and temperature behavior of the spontaneous polarization for several types of stable domain structures in thin BaTiO3 films, such as periodic quasi-2D stripe domains and arrays of 3D vortex-antivortex pairs, whose axes are perpendicular to the film surface.

II. Problem Statement

Using a finite element modelling (FEM) and Landau-Ginzburg-Devonshire (LGD) phenomenological approach combined with electrostatic equations and elasticity theory, we model the polarization, internal electric field, elastic stresses, and strains in a thin BaTiO3 film. We consider a BaTiO3 film [001] sandwiched between two ultra-thin paraelectric dead layers with a high relative dielectric permittivity. A top electrode and conducting substrate are in contact with the layers, and the voltage is applied between the conductors (see Fig. 1a). The dead layers are required for the thermodynamic stability of the domain structure. Only a single-domain distribution is stable in the case of perfect electric contact between the ferroelectric film surfaces and ideal conducting electrodes. The film thickness varies from 4 to 20 nm, and the thickness of each dead layer varies from 0 to 0.8nm. The relative dielectric permittivity is 300.

For most cases, we impose a slight tensile strain 0 < 𝑢 𝑚 <0.25% at the interface BaTiO3 filmsubstrate. The strain leads to the disappearance of the metastable orthorhombic phase and shifts the temperatures of structural and polar transitions [START_REF] Pertsev | Effect of mechanical boundary conditions on phase diagrams of epitaxial ferroelectric thin films[END_REF]. Also, the strain supports a rhombohedral ferroelectric phase in a single-domain film at temperatures below 360 K. At the same conditions, the formation of in-plane domains is favorable in the film due to depolarization effects. The transition from the rhombohedral to the tetragonal phase occurs at a temperature above 360 K, while the transition from the tetragonal ferroelectric to the paraelectric cubic phase occurs above 420 K (see Fig. 1b).

An LGD free energy functional 𝐺 of the BaTiO3 film includes a Landau energyan expansion on 2-4-6 powers of the polarization components 𝑃 𝑖 , 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 ; a polarization gradient energy, 𝐺 𝑔𝑟𝑎𝑑 ; the electrostatic energy, 𝐺 𝑒𝑙 ; an elastic, electrostriction contribution 𝐺 𝑒𝑠 , a flexoelectric contribution, 𝐺 𝑓𝑙𝑒𝑥𝑜 ; and a surface energy term, 𝐺 𝑆 [START_REF] Morozovska | Flexocoupling impact on the size effects of piezo-response and conductance in mixed-type ferroelectrics-semiconductors under applied pressure[END_REF][START_REF] Eliseev | Defect driven flexo-chemical coupling in thin ferroelectric films[END_REF]. The free energy functional, the Euler-Lagrange equations obtained from its variation, a mathematical formulation of the electrostatic and elastic sub-problem, and FEM simulation details are given in Appendix A of Supplementary Materials. The ferroelectric, dielectric, and elastic properties of the BaTiO3 [001] are listed in Table AI therein.

For the initial distribution of the polarization in the film, we used either a periodic 90-degree zigzag pattern of in-plane domain stripes or a regular arrangement of crossed 90-degree domain walls in the XZ-plane (see Fig. A1a-c in Appendix A). When we used a purely random noise as the initial distribution of polarization at room temperature, it relaxed to in-plane a-domain stripes without any fluxclosure at the film surfaces. This circumstance can indicate that a-domain stripes are the energetic ground state configuration for a film in the rhombohedral phase. In fact, the energy density of the periodic stripe domains is -1.560 MPa, while the energy density of the vortex-antivortex pairs is -1.373 MPa (for the same polarity of antivortex cores) or -1.368 MPa (for the opposite polarity of antivortex cores) for a 10nm BaTiO3 film at room temperature.

FIGURE 1. (a)

The ferroelectric film placed between two paraelectric dead layers with a very high dielectric permittivity. Dead layers thickness is 0.8 nm and relative dielectric permittivity is 300. The mismatch strain exists at the interface "dead layer 1conducting substrate". The coordinate frame we are using is shown in the bottom left corner. (b) Phase diagram of a thick single-domain BaTiO3 film on a rigid substrate. PE and FE denote the paraelectric and ferroelectric phases, respectively. Indices "c", "r" and "aa" corresponding to the tetragonal cphase with an out-of-plane polarization, rhombohedral phase all three components of polarization, and a-phase with two in-plane components of polarization, respectively.

III. Flexoelectric Coupling Influence on the Polarization Distribution

Our FEM studies revealed that the out-of-plane polarization component 𝑃 𝑦 can be very sensitive to the flexoelectric coupling (shortly: "flexocoupling") for both stripe domains and vortices. However, the influence is rather different for each of them.

The influence of the flexocoupling on the in-plane periodic stripe domains in a 5-nm BaTiO3 film is shown in Fig. 2. The dependence of the polarization magnitude 𝑃 = |𝑃 ⃗ | on 𝐹 𝑖𝑗 is negligible, and thus Fig. 2a is the same for positive, zero, and negative 𝐹 𝑖𝑗 . This happens because the value of 𝑃 is governed by the in-plane polarization components, 𝑃 𝑥 and 𝑃 𝑧 , and the amplitude of 𝑃 𝑥 (about 20 C/cm 2 ) is almost flexo-insensitive (see Fig. 2e). Although the flexocoupling modulates the amplitude of 𝑃 𝑧 , the modulation amplitude is very small (less than 0.6 C/cm 2 ) in comparison with its almost constant value ~19 C/cm 2 (see Fig. 2g). The flexocoupling significantly increases the amplitude of a small out-of-plane polarization component 𝑃 𝑦 (from 0.1 C/cm 2 to 0.4 C/cm 2 ), and the "up" or "down" direction of 𝑃 𝑦 is defined by the sign of flexocoupling coefficients 𝐹 𝑖𝑗 . The amplitude of 𝑃 𝑦 is proportional to the strength |𝐹 𝑖𝑗 |. This effect is illustrated qualitatively in Fig. 2b, 2c and 2d plotted for negative, zero and positive 𝐹 𝑖𝑗 , respectively. The quantitative influence of 𝐹 𝑖𝑗 on 𝑃 𝑦 is clearly seen from z-profiles of 𝑃 𝑦 in Fig. 2f are calculated at the film top surface for negative (blue curves), zero (black curves), and positive (red curves) flexoelectric coefficients 𝐹 𝑖𝑗 . Note that the small value of 𝑃 𝑦 at 𝐹 𝑖𝑗 = 0 results from the electrostrictive coupling, and its phase depends on the initial conditions in the case.

Let us underline that the flexocoupling influence on 𝑃 𝑖 revealed in this work and shown in Fig. 2, is in a qualitative agreement with the earlier predictions for the flexocoupling-induced Néel component of polarization at a "nominally uncharged" single domain wall in a bulk rhombohedral BaTiO3 [START_REF] Eliseev | Structural phase transitions and electronic phenomena at 180-degree domain walls in rhombohedral BaTiO3[END_REF] and BiFeO3 [START_REF] Morozovska | Anisotropic conductivity of uncharged domain walls in BiFeO3[END_REF]. However, the flexocoupling influence on the periodic stripe a-domains in thin films is quantitatively different from that found in the case of single walls in a bulk rhombohedral material [START_REF] Morozovska | Anisotropic conductivity of uncharged domain walls in BiFeO3[END_REF][START_REF] Yudin | Bichiral structure of ferroelectric domain walls driven by flexoelectricity[END_REF][START_REF] Eliseev | Structural phase transitions and electronic phenomena at 180-degree domain walls in rhombohedral BaTiO3[END_REF]. Using the results of Ref. [START_REF] Morozovska | Anisotropic conductivity of uncharged domain walls in BiFeO3[END_REF], an approximate analytical expression for the Néel-type component -c).

𝑃 𝑦 ≈ 𝜀 0 𝜀 𝑏 𝑓 𝑄 1+2𝛽𝜀 0 𝜀 𝑏 𝜕𝑃 𝑥 2 𝜕𝑧 + 𝜀 0 𝜀 𝑏 𝑞𝑃 𝑥 1+2𝛽𝜀 0 𝜀 𝑏 (𝑃 𝑠 2 -𝑃 𝑥 2 ) (1)
Besides topological aspects, vortices and antivortices also differ concerning their energy.

Although the correlation energy is similar in both cases, vortex structures are virtually divergence-free, such that they tend to form spontaneously in inhomogeneous structures to reduce the dipolar energy. By contrast, antivortices display strong divergences, leading to a distinct distribution of bound charges and a corresponding increase in dipolar energy. This behavior is known from magnetic structures, and the situation in ferroelectrics is analogous: here the divergence 𝑑𝑖𝑣 𝑃 ⃗ is relatively small in the vortex core case and quite significant in the antivortex one (see Fig. A1d in Appendix A). Because of these differences in energy, individual antivortices tend to dissolve in magnetic systems, and particular geometric shapes are required to stabilize them [START_REF] Shigeto | Magnetic force microscopy observation of antivortex core with perpendicular magnetization in patterned thin film of permalloy[END_REF]. Antivortex structures typically occur for topological reasons as they form a link between two neighboring vortices with equal circulation. In magnetic structures, such situations are well-known in the form of "cross-tie" domain walls [START_REF] Middelhoek | Domain walls in thin Ni-Fe films[END_REF], which are essentially linear chains of alternating vortex-antivortex structures. However, we are not aware of a generalization of cross-tie domain wall structures to two-dimensional checkerboard-type vortexantivortex arrays.

It appears that the dependence of the polarization magnitude 𝑃 on 𝐹 𝑖𝑗 is very weak, and so the images in Fig. 3d and4c, as well as 𝑃-profiles shown in Figs. 5a and5b, are almost 𝐹 𝑖𝑗 -independent. In contrast to the magnitude, the distribution of the out-of-plane polarization component 𝑃 𝑦 does depend on 𝐹 𝑖𝑗 , i.e., it is "flexo-sensitive", but the sensibility is only significant for a vortex region (as will be shown below).

The antivortex has a wide and smooth dipolar core, whose rounded rectangular-like shape and other features are weakly sensitive to the flexocoupling (see Figs. 3a-c and Fig. B2a-c). The direction of the component 𝑃 𝑦 in the antivortex core is defined by initial conditions (compare Fig. 4a and4b), but not by the 𝐹 𝑖𝑗 structure. The magnitude of 𝑃 𝑦 in the antivortex core is relatively high, reaching (10 -15) C/cm 2 for thin films (see Fig. 5c,e,g).

The vortex has a very prolate asymmetrical spike-like quadrupolar core, whose polarity is defined by the sign of 𝑃 𝑦 , is controlled by the flexocoupling (compare blue and reddish vortex cores in Figs. 3a-c and the quadrupolar areas inside dotted ellipses in Figs. 4d-f). The sign of 𝐹 𝑖𝑗 determines the sign of the polarization and its value in the vortex core, as it is seen from the color images in Figs. 4d-f, and especially from the 𝑦-profiles in Fig. 5d, 5f and 5h calculated for 𝐹 𝑖𝑗 < 0, 𝐹 𝑖𝑗 = 0 and 𝐹 𝑖𝑗 > 0, respectively. The vortex core polarization becomes significantly smaller for 𝐹 𝑖𝑗 = 0 (see Fig. 5f). At the nominal values of 𝐹 𝑖𝑗 (listed in Table AI) the maximal polarization in the vortex core reaches (1 -1.5) C/cm 2 , depending on the film surface (see Fig. 5d, 5f and5h). An axial asymmetry of the quadrupolar core is due to different elastic boundary conditions at the mechanically free top surface, where the normal stress is absent, and at the bottom surface clamped to a rigid substrate with a mismatch strain. The asymmetry is clearly seen by comparing the small and high-contrast reddish part with the smooth and more extended blue part of the quadrupolar core in Figs. 4d and4f, as well as from the comparison of solid and dashed curves in Figs. 5d and5h. Notably that the transformation 𝐹 𝑖𝑗 → -𝐹 𝑖𝑗 inverts the vortex core (𝑃 𝑦 → -𝑃 𝑦 ) with respect to the y-axis. The behavior of the flexo-sensitive polarization vortex is qualitatively the same as the flexon-type polarization in cylindrical nanoparticles [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF]. We argue that the flexo-sensitive vortex cores can be interpreted as a manifestation of "flexons" [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF] in thin films. To ascertain this conjecture, a topological index of these structures is calculated and analyzed below. .

(

) 2
Similarly to the flexon discussed in Ref. [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF], 𝑛(𝑦) is a normalized profile of the out-of-plane polarization component, and the coordinate origin is placed in the vortex axis. Note that the expression ( 2) is an approximation since the distribution of 𝑃 𝑦 (𝑥, 𝑦, 𝑧) is not fully axially symmetric with respect to the vortex axis. The axial asymmetry appears in the form of a square mesh between the vortices (see e.g., Fig. 3).

Although 𝑃 𝑦 is nonzero for 𝐹 𝑖𝑗 = 0 due to the electrostriction coupling, two Bloch points (𝑃 = 0) located symmetrically under the film surfaces exist for that case (see Fig. 6d). The depth profile of 𝑃 𝑦 becomes asymmetric with respect to the film surfaces for 𝐹 𝑖𝑗 ≠ 0, and only one Bloch point exists in this case (see Fig. 6e). The transformation 𝐹 𝑖𝑗 → -𝐹 𝑖𝑗 leads to the transformation 𝑃 𝑦 → -𝑃 𝑦 (compare red and blue curves in Fig. 6e). Hence, the sign of 𝑃 𝑦 , and so the sign of 𝑛(𝑦), is defined by the sign of 𝐹 𝑖𝑗 , since the flexo-induced contribution to the polarization component 𝑃 𝑦 significantly dominates over the electrostriction contribution. The result corroborates the "flexon" character [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF] of the vortex polarization. The dependence 𝑛(𝑦) is shown in Fig. 6e for zero, positive and negative 𝐹 𝑖𝑗 , respectively.

Since the value 𝑃(0, 𝑦, 0) coincides with |𝑃 𝑦 (0, 𝑦, 0)| at the vortex axis, and 𝑃 𝑦 (0, 𝑦, 0) = 𝑃(0, 𝑦, 0) = 0 in the Bloch point, the topological index jumps in that point from -½ to +½ (red curve in Fig. 6e), or from +½ to -½ (blue curve in Fig. 6e), depending on the 𝐹 𝑖𝑗 sign. In any case, 𝑛(𝑦) = ± 1 2 ⁄ at the film surfaces. The result demonstrates the surface localization of the flexo-sensitive vortices, similarly to the "edge" localization of flexons in cylindrical nanoparticles [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF]. The topological index, which can be interpreted as the degree to which a structure is chiral, changes sign from one surface to the other, and To the best of our understanding, the experiment [START_REF] Balke | Enhanced electric conductivity at ferroelectric vortex cores in BiFeO3[END_REF] indicates the possible existence of flexosensitive vortices. However, most phase-field simulations (see e.g., [START_REF] Winchester | Electroelastic fields in artificially created vortex cores in epitaxial BiFeO3 thin films[END_REF]) either did not include the flexoelectric effect, or they reveal other trends [START_REF] Li | Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling[END_REF]. Indeed, the flexo-sensitive vortices are "hidden" by flexo-insensitive antivortices, whose out-of-plane polarization is much higher. Moreover, when we increase the transverse size of the computational cell, the core of the antivortex becomes a complex structure resembling an r-phase domain. Still, the polarization remains strictly vertical in the core region, as for the c-phase domain. To define the temperature interval, where the flexo-sensitive stripes or vortices are stable or metastable, we performed FEM studies of the polarization distribution in the film in the temperature range from 200 K to 440 K. A stripe-domain configuration with a flexo-sensitive out-of-plane polarization component 𝑃 𝑦 is stable at temperatures lower than 415 K; at higher temperatures it becomes faint as the film transforms to the paraelectric phase.

A vortex-like configuration with a flexo-sensitive out-of-plane polarization component 𝑃 𝑦 is metastable at temperatures lower than 410 K (we use the word "metastable" here to underline that the domain stripes have lower energy). The flexo-insensitive antivortex core disappears above 320 K. The dipolar core of the vortex gradually disappears above 415 K, when the vortex-antivortex structure transforms into a circular vortex-antivortex structure without both cores, which exists up to the paraelectric transition at about 420 K. The relatively wide temperature range (200 < 𝑇 < 410 K) of the flexo-sensitive vortices meta-stability gives us the hope that the domain morphology can be observed experimentally in thin ferroelectric films by PFM, cAFM [START_REF] Maksymovych | Tunable metallic conductance in ferroelectric nanodomains[END_REF][START_REF] Balke | Enhanced electric conductivity at ferroelectric vortex cores in BiFeO3[END_REF][START_REF] Vasudevan | Domain wall geometry controls conduction in ferroelectrics[END_REF] and nonlinear optical microscopy [START_REF] Cherifi-Hertel | Non-ising and chiral ferroelectric domain walls revealed by nonlinear optical microscopy[END_REF] methods.

The disappearance of antivortex cores when the temperature is decreased below (315 -320) K is related to the transition from the monoclinic r-phase (with all three polarization components) to the orthorhombic a-phase (with two in-plane components of polarization). The monoclinic-orthorhombic transition occurs at about 340 K for thick films; and its temperature decrease to 320 K in thin films with dead layers is caused by the depolarization field arising from the layers.

IV. Conclusion

In this theoretical work, we have considered the influence of the flexocoupling on the spatial distribution and the temperature behavior of the ferroelectric polarization for different types of stable domain structures in thin ferroelectric films, such as periodic stripe domains and arrays of vortexantivortex pairs. Our FEM simulations and analytical calculations reveal that an out-of-plane polarization component can be very sensitive to the flexocoupling for both stripe domains and vortices. However, the influence is rather different. Namely, the flexocoupling significantly increases the amplitude of a small out-of-plane polarization component in the stripe domains. The "up" or "down" direction of this out-ofplane component depends on the sign of the flexocoupling coefficients. The vortex has a spike-like quadrupolar core controlled by the flexocoupling, whereas the antivortex has a wide smooth dipolar core, whose shape and other features are weakly insensitive to the coupling. The origin of the flexo-sensitive vortex is due to the system's tendency to minimize its elastic and electrostatic energy, because a vortex domain wall (in fact, a flux-closure domain) creates a much weaker depolarization field compared to the field of a charged domain wall.

We interpret the flexo-sensitivity of the vortex cores as a variant of the effects that stabilize the recently revealed "flexons" [START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF] in thin films. To the best of our knowledge, such flexo-sensitive vortices have not been reported previously, while several experiments [START_REF] Balke | Enhanced electric conductivity at ferroelectric vortex cores in BiFeO3[END_REF][START_REF] Li | Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling[END_REF] suggest their existence.

The relatively wide temperature range (200 < 𝑇 < 410 K) in which the flexo-sensitive vortices are meta-stable gives us hope that these structures can be observed experimentally in thin ferroelectric films by PFM, cAFM, and nonlinear optical microscopy methods. However, in thin films, the PFM observation of the flexo-sensitive vortices is complicated by the presence of much higher contrast antivortices, which are flexo-insensitive.

Acknowledgements. This material is based upon work (S.V.K.) supported by the U.S. The boundary condition for polarization at the film-dead layer interface accounts for the flexoelectric effect:

𝑎 𝑖𝑗 (𝑆) 𝑃 𝑗 + (𝑔 𝑖𝑗𝑘𝑙 𝜕𝑃 𝑘 𝜕𝑥 𝑙 -𝐹 𝑘𝑙𝑖𝑗 𝜎 𝑘𝑙 ) 𝑛 𝑗 | 𝑥 3 =0,ℎ = 0 (A.3)
where n is the outer normal to the surface, i = 1, 2, 3. In our FEM studies, we use the so-called "natural" boundary conditions corresponding to 𝑎 𝑖𝑗 (𝑆) = 0. To model the film infinity in the lateral directions "1"

and "2", The periodic boundary conditions are valid at the sidewalls of the computation region.

The electric field components 𝐸 𝑖 in Eqs.(A.2) are derived from the electric potential 𝜑 in a conventional way, 𝐸 𝑖 = -𝜕𝜑 𝜕𝑥 𝑖 ⁄ . The potential 𝜑 𝑓 satisfies the Poisson equation in the ferroelectric film (subscript "f"):

𝜀 0 𝜀 𝑏 ( 𝜕 2 𝜕𝑥 1 2 + 𝜕 2 𝜕𝑥 2 2 + 𝜕 2 𝜕𝑥 3 2 ) 𝜑 𝑓 = 𝜕𝑃 𝑖 𝜕𝑥 𝑖 , 0 ≤ 𝑥 3 ≤ ℎ. (A.4a)
The electric potential 𝜑 𝑑 in the dead layer satisfies the Laplace equation (subscript "d"):

𝜀 0 𝜀 𝑒 ( 𝜕 2 𝜕𝑥 1 2 + 𝜕 2 𝜕𝑥 2 2 + 𝜕 2 𝜕𝑥 3
2 ) 𝜑 𝑑 = 0, -𝑑 < 𝑥 3 < 0 ∪ ℎ < 𝑥 3 < ℎ + 𝑑. (A.4b) Equations (A.4) are supplemented with the continuity conditions for electric potential and the normal components of the electric displacements at the film surfaces:

(𝜑 𝑑 -𝜑 𝑓 )| 𝑥 3 =0,ℎ = 0, 𝒏(𝑫 𝑑 -𝑫 𝑓 )| 𝑥 3 =0,ℎ = 0. (A.4a)
to the unit cell size, 0.4 nm, the maximal size is 1.6 nm. The dependence on the mesh size is verified by increasing the minimal size to 0.8 nm. We verified that this only results in minor changes in the electric polarization, electric field, and elastic stress and strain, such that the spatial distribution of each of these quantities becomes less smooth (i.e., they contain numerical errors in the form of a small random noise).

However, when using these larger cell sizes, all significant details remain visible, and more importantly, the system energy remains essentially the same with an accuracy of about 0.1%.

LGD coefficients and other material parameters of BaTiO3 are listed in Table AI. Polarization components in the initial distribution, whose relaxation leads to the lateral grid of vortices and antivortices, are shown in Fig. A1. Note that we superimposed a slight random noise to prevent the formation of unstable equilibrium states arising from idealized initial conditions.

Table AI.

LGD coefficients and other material parameters of BaTiO3

FIGURE 2 .

 2 FIGURE 2. Regular stripe domains inside a thin strained BaTiO3 film. Polarization magnitude 𝑃 (a) and out-ofplane component 𝑃 𝑦 calculated for negative (b), zero (c) and positive (d) flexoelectric coefficients 𝐹 𝑖𝑗 . (e)-(g) Polarization z-profiles calculated at the film top surface for negative (blue curves), zero (black curves) and positive (red curves) flexoelectric coefficients 𝐹 𝑖𝑗 . The film thickness is 5 nm, mismatch strain 𝑢 𝑚 =0.2%, and the temperature 𝑇 is 300 K. The dependence of 𝑃 on 𝐹 𝑖𝑗 is negligible and not shown in the figure. Arrows in the plot (a) show the direction of polarization vector.

2 )

 2 Polarization profiles: blue -𝐹 𝑖𝑗 < 0, black -𝐹 𝑖𝑗 = 0, red -𝐹 𝑖𝑗 > 0.The influence of the flexocoupling on a periodic lattice of polarization vortices and antivortices is shown in Figs.3-5and Figs.B1-3(in the Supplement) for a 9-nm BaTiO3 film stretched with a misfit strain 0.1%. The computation cell, which translates periodically, is a rectangular parallelepiped consisting of two vortices separated by two antivortices (see Fig.3and Fig.B1). The vortex and antivortex can be topologically distinguished by the winding number 𝑊 = ∮ 𝛼(𝜑) 2𝜋 𝑑𝑆, that is the normalized line integral on a closed loop S over the 𝛼(𝜑) that the in-plane component of the polarization 𝑃 ⃗ encloses with the x-axis[START_REF] Hertel | Exchange Explosions: Magnetization Dynamics during Vortex-Antivortex Annihilation[END_REF]. The winding number counts the number of vortex (𝑊 = +1) and antivortex (𝑊 = -1) structures within the loop 𝑆 (see Figs.6a

FIGURE 3 .and Fig. 4 )

 34 FIGURE 3. The XZ-sections of the polarization out-of-plane component 𝑃 𝑦 at the film top surface calculated for negative (a), zero (b) and positive (c) flexoelectric coefficients 𝐹 𝑖𝑗 . The XZ-sections of the polarization magnitude 𝑃 (d) and hydrostatic stress 𝜎 (e, f), which are almost flexo-independent. Arrows show polarization direction in the vortex-antivortex region. The thickness of BaTiO3 film is 9 nm, mismatch strain 𝑢 𝑚 =0.1%, and 𝑇 = 300 K.

FIGURE 4 .

 4 FIGURE 4. The 𝜉𝑌-sections of the out-of-plane polarization component 𝑃 𝑦 for different seedings (a-b) and the magnitude 𝑃 (c), which are almost flexo-independent. Flexo-sensitive 𝜁𝑌-sections of the 𝑃 𝑦 calculated for negative (d), zero (e) and positive (f) flexoelectric coefficients 𝐹 𝑖𝑗 . Coordinates 𝜉 = (𝑥 + 𝑧)/√2 and 𝜁 = (𝑥 -𝑧)/√2. Arrows in the plots (a)-(c) show polarization direction. The thickness of BaTiO3 film is 9 nm, mismatch strain 𝑢 𝑚 =0.1%, and 𝑇 = 300 K.

μC/cm 2 (FIGURE 5 .

 25 FIGURE 5. Profiles of the flexo-independent antivortices (left column) and flexo-sensitive vortices (right column) at the top (solid curves) and bottom (dashed curves) surfaces of a thin strained BaTiO3 film. (a)-(b) Flexoindependent 𝜉and 𝜁-profiles of the polarization magnitude 𝑃. (c)-(h) Flexo-sensitive 𝜉and 𝜁-profiles of the out-

  polarization orientation and the integration is performed over the vortex cross-section 𝑆 = {𝑥, 𝑧}, quantifies the chirality of the polarization structure. Using the method explained in Appendix D of Ref.[START_REF] Morozovska | Introducing the Flexon -a New Chiral Polarization State in Ferroelectrics[END_REF], the approximate y-dependence of the topological index is 𝑛(𝑦) ≈ -𝑃 𝑦 (𝑥=0,𝑦,𝑧=0) 2𝑃(𝑥=0,𝑦,𝑧=0)

  changes sign upon reversal of the sign of 𝐹 𝑖𝑗 . Since |𝑃 𝑦 |~|𝐹 12 |, the index 𝑛(𝑦) increases in magnitude with increasing absolute value of |𝐹 𝑖𝑗 |. These properties make obvious the clear correlation between the flexocoupling and the formation of chiral polar vortices in thin ferroelectric films.

FIGURE 6 .

 6 FIGURE 6. (a) A ferroelectric vortex (𝑊 = +1) and (b) antivortex (𝑊 = -1), whose topology can be distinguished by the winding number 𝑊 (c). Part (c) is redrawn from Ref. [54]. (d-e) Depth profiles of the out-ofplane polarization component 𝑃 𝑦 (solid curves) and magnitude 𝑃 (dotted curves) calculated at the vortex axis for zero (black curves, d), negative (blue curves, e) and positive (red curves, e) flexoelectric coefficients 𝐹 𝑖𝑗 . Circles are Bloch points, where 𝑃 = 0. (f) The y-profile of the polarization's topological index 𝑛(𝑦) for zero (black dotted lines), positive (solid red lines), and negative (solid blue lines) 𝐹 𝑖𝑗 . The film thickness is 9 nm, mismatch strain 𝑢 𝑚 =0.1%, and 𝑇 = 300 K.

2 FijF11 = 2 . 4 ,

 224 b = 7 (core background) e = 10 (surrounding) ai (in mJ/C 2 ) a1 = 3.34(T-381)10 5 , T = 3.3410 5 (a1 = -2.9410 7 at 298 K) aij (in m 5 J/C 4 ) a11 = 4.69(T-393)10 6 -2.0210 8 , a12 = 3.23010 8 , (a11 = -6.7110 8 at 298 K) aijk (in m 9 J/C 6 ) a111 = -5.52(T-393)10 7 +2.7610 9 , a112 = 4.4710 9 , a123 = 4.9110 9 (at 298 K a111 = 82.810 8 , a112 = 44.710 8 , a123 = 49.110 8 ) Qij (m 4 /C 2 ) Q11=0.11, Q12= -0.043, Q44=0.059 sij (in 10 -12 Pa -1 ) s11=8.3, s12= -2.7, s44=9.24 gij (in 10 -10 m 3 J/C 2 ) g11=5.0, g12= -0.2, g44= 0.(in 10 -11 m 3 /C) fij (in V) F12 = 0.5, F44 = 0.6 (the first two values are recalculated from [a] values f11 = 5.1, f12 = 3.3, f44 = 0.065 V. The equality 𝐹 44 = 𝐹 11 -𝐹 12 is valid in the isotropic case. [a] I. Ponomareva, A.K. Tagantsev, L. Bellaiche, Finite-temperature flexoelectricity in ferroelectric thin films from first principles, Phys. Rev. B 85 (2012) 104101.Analytical calculations are performed and visualized in the Mathematica 12.2 software (https://www.wolfram.com/mathematica).

Figure A1 .

 A1 Figure A1. (a)-(c) Polarization components in the initial distribution, whose relaxation leads to the lateral grid of vortices and antivortices. Two red spots in the plot (c) are the seeding of antivortex cores. Arrows in the plot (c) show the polarization direction. The small random component superimposed on the regular structure is not shown. (d) Polarization divergence is the relaxed structure consisting of two vortex-antivortex pairs (final state).

Figure B1 .Figure B2 .Figure B3 .

 B1B2B3 Figure B1. XZ-section of the out-of-plane polarization component 𝑃 𝑦 at the film top (a-c) and bottom (d-f) surfaces calculated for negative (a,d), zero (b,e) and positive (c,f) flexoelectric coefficients 𝐹 𝑖𝑗 . Arrows in the plots show the direction of polarization in the vortex-antivortex region. The thickness of BaTiO3 film is 9 nm, mismatch strain 𝑢 𝑚 =0.1%, and the temperature is 300 K.

  can be used as an estimate in the vicinity of the domain walls. Here the first term has flexoelectric nature, and the second term originates from electrostriction coupling. Constants 𝑓 𝑄 ≈ 𝐹 12 𝑄 11 𝑠 12 -𝑠 11 𝑄 12

			𝑠 11 2 -𝑠 12 2	, 𝛽 ≈
	𝑎 1 + (𝑎 12 + 𝑄 44 2 2𝑠 44 ⁄	)𝑃 𝑠 2 , and 𝑞 ≈ -	𝑄 11 𝑠 12 -𝑠 11 𝑄 12 𝑠 11 2 -𝑠 12

2

𝑄

[START_REF] Gu | Non-Ising Character of a Ferroelectric Wall Arises from a Flexoelectric Effect[END_REF]

, where 𝐹 12 is the flexoelectric coefficient, 𝑄 𝑖𝑗 are electrostriction coefficients, 𝑠 12 are elastic compliances, 𝑎 1 and 𝑎 12 are LGD expansion coefficients, 𝑃 𝑠 is the spontaneous polarization, 𝜀 𝑏 is a background permittivity

[START_REF] Tagantsev | Interface-induced phenomena in polarization response of ferroelectric thin films[END_REF]

, and 𝜀 0 is the vacuum permittivity (see Table

AI

).

Mathematical Formulation of the Problem and FEM Details A. Mathematical Formulation of the Problem The

  Department of Energy, Office of Science, Office of Basic Energy Sciences, and performed at the Center for Nanophase Materials Sciences, a US Department of Energy Office of Science User Facility. A portion of FEM was conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility (CNMS Proposal ID: 257). A.N.M. work is supported by the National Academy of Sciences of Ukraine (the Target Program of Basic Research of the National Academy of Sciences of Ukraine LGD free energy functional G additively includes a Landau expansion on powers of 2-4-6 of the polarization (𝑃 𝑖 ), 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 ; a polarization gradient energy contribution, 𝐺 𝑔𝑟𝑎𝑑 ; an electrostatic contribution, 𝐺 𝑒𝑙 ; the elastic, electrostriction, flexoelectric contributions, 𝐺 𝑒𝑠+𝑓𝑙𝑒𝑥𝑜 ; and a surface energy, Here 𝑉 𝑓 and 𝑉 𝑑 are the film and dead layer volume, respectively. The coefficient 𝑎 𝑖 linearly depends on temperature T, 𝑎 𝑖 (𝑇) = 𝛼 𝑇 [𝑇 -𝑇 𝐶 ], where 𝛼 𝑇 is the inverse Curie-Weiss constant and 𝑇 𝐶 is the ferroelectric Curie temperature renormalized by electrostriction and surface tension. Tensor components 𝑎 𝑖𝑗 are regarded as temperature-independent. The tensor 𝑎 𝑖𝑗 is negatively defined for the BaTiO3undergoing the first order transition to the paraelectric phase. The higher nonlinear tensor 𝑎 𝑖𝑗𝑘 and the gradient coefficients tensor 𝑔 𝑖𝑗𝑘𝑙 are positively defined and regarded as temperature-independent. The following designations are used in Eq.(A.1e): 𝜎 𝑖𝑗 is the stress tensor, 𝑠 𝑖𝑗𝑘𝑙 is the elastic compliances tensor, 𝑄 𝑖𝑗𝑘𝑙 is the electrostriction tensor, and 𝐹 𝑖𝑗𝑘𝑙 is the flexoelectric tensor. 𝑎 1 -𝑄[START_REF] Tanygin | Symmetry theory of the flexomagnetoelectric effect in the magnetic domain walls[END_REF] (𝜎 11 + 𝜎 33 ) -𝑄 11 𝜎 22 ) -𝑄 44 (𝜎 12 𝑃 1 + 𝜎 23 𝑃 3 ) + 4𝑎 11 𝑃 2 The temperature-dependent Khalatnikov coefficient 𝛤 [60] determines the relaxation time of the polarization 𝜏 𝐾 = 𝛤 |𝛼| ⁄ . Consequently, 𝜏 𝐾 typically varies in the range (10 -9 -10 -6 ) seconds for temperatures far from TC. As argued by Hlinka et al. [61], we assumed that 𝑔 44 ′ = -𝑔 12 in Eqs.(A.2).

	𝛤	𝜕𝑃 2 𝜕𝑡	Supplementary Materials + 2𝑃 2 (3 + 2𝑎 12 𝑃 2 (𝑃 1 2 + 𝑃 3 2 ) +
	6𝑎 111 𝑃 2 5 + 2𝑎 112 𝑃 2 (𝑃 1 4 + 2𝑃 2 2 𝑃 1 2 + 𝑃 3 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 2 𝑃 1 2 𝑃 3 2 -𝑔 11 -𝐹 11 𝜕𝜎 22 𝜕𝑥 2 -𝐹 12 ( 𝜕𝜎 11 𝜕𝑥 2 + 𝜕𝜎 33 𝜕𝑥 2 𝜕𝑃 3 𝜕𝑡 + 2𝑃 3 (𝑎 1 -𝑄 12 (𝜎 11 + 𝜎 22 ) -𝑄 11 𝜎 33 ) -𝑄 44 (𝜎 13 𝑃 1 + 𝜎 23 𝑃 2 ) + 4𝑎 11 𝑃 3 𝜕 2 𝑃 2 𝜕𝑥 2 2 -𝑔 44 ( 𝜕 2 𝑃 2 𝜕𝑥 1 2 + ) -𝐹 44 ( 𝜕𝜎 12 𝜕𝑥 1 + 𝜕𝜎 23 𝜕 2 𝑃 2 𝜕𝑥 3 2 ) = 𝜕𝑥 3 ) + 𝐸 2 (A.2b) 3 + 2𝑎 12 𝑃 3 (𝑃 1 2 + 𝑃 2 2 ) + APPENDIX A. 𝐺 𝑆 . It has the form: 𝛤 𝐺 = 𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 + 𝐺 𝑔𝑟𝑎𝑑 + 𝐺 𝑒𝑙 + 𝐺 𝑒𝑠+𝑓𝑙𝑒𝑥𝑜 + 𝐺 𝑆 , (A.1a) 6𝑎 111 𝑃 3 5 + 2𝑎 112 𝑃 3 (𝑃 1 4 + 2𝑃 3 2 𝑃 1 2 + 𝑃 2 4 + 2𝑃 2 2 𝑃 3 2 ) + 2𝑎 123 𝑃 3 𝑃 1 2 𝑃 2 2 -𝑔 11 𝜕 2 𝑃 3 𝜕𝑥 3 2 -𝑔 44 ( 𝜕 2 𝑃 3 𝜕𝑥 1 2 + 𝜕 2 𝑃 3 𝜕𝑥 2 2 ) =
			𝐺 𝐿𝑎𝑛𝑑𝑎𝑢 = ∫ 𝑑 3 𝑟 𝑉 𝑓	[𝑎 𝑖 𝑃 𝑖 2 + 𝑎 𝑖𝑗 𝑃 𝑖 2 𝑃 𝑗 2 + 𝑎 𝑖𝑗𝑘 𝑃 𝑖 2 𝑃 𝑗 2 𝑃 𝑘 2 ], -𝐹 11 𝜕𝜎 33 𝜕х 3 -𝐹 12 ( 𝜕𝜎 11 𝜕х 3 + 𝜕𝜎 22 𝜕х 3 ) -𝐹 44 ( 𝜕𝜎 13 𝜕𝑥 1	+	𝜕𝜎 23 𝜕𝑥 2	(A.1b) ) + 𝐸 3
					𝐺 𝑔𝑟𝑎𝑑 = ∫ 𝑑 3 𝑟 𝑉 𝑓	𝑔 𝑖𝑗𝑘𝑙 2	𝜕𝑃 𝑖 𝜕𝑥 𝑗	𝜕𝑃 𝑘 𝜕𝑥 𝑙	,	(A.1c) (A.2c)
			𝐺 𝑒𝑙 = -∫ 𝑑 3 𝑟 𝑉 𝑓	(𝑃 𝑖 𝐸 𝑖 +	𝜀 0 𝜀 𝑏 2	𝐸 𝑖 𝐸 𝑖 ) -	𝜀 0 2 ∫ 𝜀 𝑖𝑗 𝑑 𝐸 𝑖 𝐸 𝑗 𝑑 3 𝑟 𝑉𝑑	,	(A.1d)
			𝐺 𝑒𝑠+𝑓𝑙𝑒𝑥𝑜 = -∫ 𝑑 3 𝑟 𝑉 𝑓	[	𝑠 𝑖𝑗𝑘𝑙 2	𝜎 𝑖𝑗 𝜎 𝑘𝑙 + 𝑄 𝑖𝑗𝑘𝑙 𝜎 𝑖𝑗 𝑃 𝑘 𝑃 𝑙 + 𝐹 𝑖𝑗𝑘𝑙 (𝜎 𝑖𝑗	𝜕𝑃 𝑘 𝜕𝑥 𝑙	-𝑃 𝑘	𝜕𝜎 𝑖𝑗 𝜕𝑥 𝑙	)]	(A.1e)
				𝐺 𝑆 =	1 2 ∫ 𝑑 2 𝑟 𝑆	𝑎 𝑖𝑗 (𝑆) 𝑃 𝑖 𝑃 𝑗 .	(A.1f)
			Allowing for the Khalatnikov mechanism of polarization relaxation [ 59 ], minimization of the free
	energy (A.1) with respect to polarization leads to three coupled time-dependent Euler-Lagrange
	equations for polarization components inside the film,	𝛿𝐺 𝛿𝑃 𝑖	= -𝛤	𝜕𝑃 𝑖 𝜕𝑡	, where 𝑖 = 1, 2, 3. The explicit form
	of the equations for a ferroelectric crystal with m3m parent symmetry is:
	𝛤	𝜕𝑃 1 𝜕𝑡	+ 2𝑃 1 (𝑎 1 -𝑄 12 (𝜎 22 + 𝜎 33 ) -𝑄 11 𝜎 11 ) -𝑄 44 (𝜎 12 𝑃 2 + 𝜎 13 𝑃 3 ) + 4𝑎 11 𝑃 1 3 + 2𝑎 12 𝑃 1 (𝑃 2 2 + 𝑃 3 2 ) +
		6𝑎 111 𝑃 1 5 + 2𝑎 112 𝑃 1 (𝑃 2 4 + 2𝑃 1 2 𝑃 2 2 + 𝑃 3 4 + 2𝑃 1 2 𝑃 3 2 ) + 2𝑎 123 𝑃 1 𝑃 2 2 𝑃 3 2 -𝑔 11	𝜕 2 𝑃 1 𝜕𝑥 1 2 -𝑔 44 ( 𝜕 2 𝑃 1 𝜕𝑥 2 2 +	𝜕 2 𝑃 1 𝜕𝑥 3 2 ) =
								-𝐹 11	𝜕𝜎 11 𝜕𝑥 1	-𝐹 12 ( 𝜕𝜎 22 𝜕𝑥 1	+	𝜕𝜎 33 𝜕𝑥 1	) -𝐹 44 ( 𝜕𝜎 12 𝜕𝑥 2	+	𝜕𝜎 13 𝜕𝑥 3	) + 𝐸 1
										(A.2a)

"Prospective basic research and innovative development of nanomaterials and nanotechnologies for 2020 -2024", Project № 1/20-H, state registration number: 0120U102306). R.H. acknowledges funding from the French National Research Agency through contract ANR-18-CE92-0052 "TOPELEC".

J. Hlinka, Mobility of Ferroelectric Domain Walls in Barium Titanate, Ferroelectrics 349, 49 (2007).

The voltage is fixed at the electrodes: 𝜑 𝑑 | 𝑥 3 =-𝑑 = 0, 𝜑 𝑑 | 𝑥 3 =ℎ+𝑑 = 𝑈.

(A.4c) and periodic boundary conditions for electric potential are valid at the sidewalls of the computation region.

Elastic equations of state follow from the variation of the energy (A.

where 𝑢 𝑖𝑗 is the strain tensor components related to the displacement components 𝑈 𝑖 in the following way: 𝑢 𝑖𝑗 = (𝜕𝑈 𝑖 /𝜕𝑥 𝑗 + 𝜕𝑈 𝑗 /𝜕𝑥 𝑖 )/2.

Equations (A.5) should be considered along with equations of mechanical equilibrium

compatibility equations, 𝑒 𝑖𝑘𝑙 𝑒 𝑗𝑚𝑛 𝜕 2 𝑢 𝑙𝑛 (𝒙) 𝜕𝑥 𝑘 𝜕𝑥 𝑚 ⁄ = 0, which are equivalent to the mechanical displacement vector 𝑈 𝑖 continuity [START_REF] Timoshenko | Theory of Elasticity[END_REF]. The boundary conditions for elastic stresses and displacements at the film surfaces are conventional continuity of elastic displacement vector and normal stress components:

The boundary conditions for elastic stresses and displacements at the dead layer surfaces account for the free top surface and mismatch strain at the substrate electrode (𝑢 11 = 𝑢 22 = 𝑢 𝑚 ):

B. Finite Element Modelling and Analytical Calculations Details

FEM simulations are performed in COMSOL@MultiPhysics software, using electrostatics, solid mechanics, and general math (PDE toolbox) modules. The size of the computational region is not less than 4040160 nm 3 , and is commensurate with the cubic unit cell constant (about 0.4 nm) of BaTiO3 at room temperature. The minimal size of a tetrahedral element in a mesh with fine discretization is equal