Shearing Deformations of β-Cristobalite-Like Boron Arsenate

Abstract

Boron arsenate, BAsO₄, is crystalline material ($I\overline{4}$ group) that was recently shown to be auxetic in its (001) plane for loading in any direction in this plane, and, which exhibits negative linear compressibility at elevated pressured in its [001] direction. This work presents and discusses the results of extensive density functional theory (DFT) based simulations aimed at studying deformations that such crystals undergo when subjected to shear loading in an attempt to obtain a better insight into the manner in which this material responds to mechanical loads. The deformations for shearing in the (001) plane are described in terms of the ‘rotating squares’ model, which was used to explain the auxeticity in the same plane where it was shown that shear loading results primarily in deformations which make the ‘squares’ become ‘parallelogram-like’ rather than rotate. This lack of rigidity in projected ‘squares’ was discussed by looking at changes in bond lengths and bond angles.

1. Introduction

Boron arsenate, BAsO₄, see Figure 1, is a spectacular material having both a history and a prehistory, particularly in terms of its synthesis and crystal properties. Prior to being characterised and noted as a material with its own importance by Schulze in the early 1930s [1,2,3], it was first mentioned by Berger in a note published a century ago which presented a number of novel reactions [4]. In one of Berger’s new reactions, arsenic was obtained by evaporation of the resultant compound formed of the mixture of arsenic acids and boric acids. Here, boron arsenate was not characterized and no particulars were given for the synthesis of the compound, as it was not the scope of their study, in which it was used as a reducible compound to form arsenic vapours. This crystal was then first characterized using X-Ray diffraction by Schulze in 1933 [1], when it was described as consisting of connected tetrahedra of BO₄ and AsO₄, which is the structural model being used in this work.

Schulze’s method for the synthesis of the crystal, was described in a later publication the following year as a solid state reaction between arsenic oxide (As₂O₅) and boron trioxide (B₂O₃) at high temperature [2]. This method produced fine single crystals with edges of approximately 0.1 mm in length. Within the same publication, Schulze reported that the crystal was tetragonal, a result which was inferred from the tetragonal growth of the crystals. This 1934 publication also made the first comparison of BAsO₄ to the structure of cristobalite, a report which led to a number of important advances throughout the years. According to Gruner [5], another article published at around this time was by Levi and Ghiron [6], where it was determined that the boron arsenate that grows from such crystallizations is in one crystalline form, i.e., it is isomorphic.

Another two methods for preparing boron arsenate were described within the same year (1934), one by Schumb and Hartford [7] and the other by Gruner [5]. In the former, 100 g of arsenic acid (H₃AsO₄) were dissolved in the minimum amount of water before any solid suspended matter was filtered off. The filtered solution was then heated to its boiling point and 20 g of boric acid (B(OH)₃) were added. The solution was kept at boiling point until a white precipitate formed. This white precipitate contained small crystals which were examined and characterized in a conglomerate form. The specifics on the size of the crystals and the characterization carried out were actually reported by Schulze in 1935 [3], citing a private communication he had with Schumb and Hartford. The method for the preparation of BAsO₄ described by Schumb and Hartford appears to be the same method which was only briefly described by Berger back in 1920 [4], and it gives important details which were not previously disclosed, thus making this synthesis reproducible.

The second method for the preparation of boron arsenate was reported by Gruner, and it involved the reaction between arsenic pentoxide and boronic acid [5]. The resulting arsenic salt product was placed in water and studied further. It was found that the boron arsenate produced, once placed in water, forms two different hydrate complexes depending on the temperature: the trihydrate (BAsO₄·3 H₂O) and the hexahydrate (BAsO₄·6 H₂O) forms. The solid forms of these two complexes have different X-ray diffraction patterns (obtained in Debye–Scherrer geometry). Both also differ from the X-ray diffraction pattern of the anhydrous form of boron arsenate [2]. Comparison with the diffraction data for corresponding boron phosphates showed a very high similarity between the powder patterns (both hydrates and anhydrates) [5]. An analysis of the literature of time, also reveals that in some of the early publications, BAsO₄ was referred to as having a low-cristobalite-like structure by Schulze and later by Nieuwenkamp and also by Shafer, which is now known to not be the case; however, they correctly noted the similarity between this structure and the structure of cristobalite in general [2,8,9].

Much more recently, a microwave synthesis version of the method by Schulze was carried out successfully by Baykal et al. (2006), which decreased both the energy and the time required to carry out such a synthesis [10].

This paper will look at a different aspect in BAsO₄ namely, from the perspective of its mechanical properties and will focus on the atomic level deformations that occur when this material, as a single crystal, is subjected to mechanical loads. Such a perspective is important in view of the anomalous negative Poisson’s ratio (auxetic) [11,12] and negative linear compressibility characteristics [12,13], properties it shares with a number of other anomalous materials, models and structures. In this respect it is important to highlight that both fields of auxetics and that of negative linear compressibility has grown substantially in the past few years, as discussed elsewhere [14,15,16]. For example, in the case of auxetics, apart from the very early and fundamental modelling of cellular solids [17,18], studies include production and/or experimental characterisation of polymeric auxetics [19,20] foams or other cellular systems [21,22,23], biomaterials [24], modelling and/or experimental characterisation of crystalline materials such as hydrophthalates [25] silicates [26] and other oxides [11,12,27], the developments include deign and optimisation of molecular-level auxetics [28,29,30,31,32,33], studying of various auxetic model structures and macromodels [34,35,36,37,38,39,40,41,42,43,44,45,46] including smart tuneable auxetics [47,48]. Reported applications of auxetics include sport [49], filtration [50,51], textiles [44] and medicine [52,53]. Similar advances have also been made in the field of negative compressibility [16,32,54,55,56,57,58,59,60], including a fundamental study that looks into stability aspects of such negative systems [61]. Crystal symmetry plays an important role on the manifestation or otherwise of auxetic behaviour and a number of key studies have focused on looking at crystalline materials with specific symmetries [62]. Worth highlighting are the studies on crystals with cubic symmetry [63,64,65,66] and tetragonal symmetry [67,68], to where BAsO₄ belongs. The theoretical framework required to study crystals in this manner is well established and explained in detail elsewhere [69].

2. Rationale

BAsO₄, see Figure 1, has a particularly interesting crystalline structure with the ‘rotating squares’ motif, well known for its auxetic characteristics [37], being projected in its (001) plane (see Figure 1b). It is hence not surprising that a number of recent studies have looked at this material to study its mechanical properties. In particular, one should mention a seminal study by Haines et al. [13] that has looked at this material as this is subjected to hydrostatic pressure. This work was able to identify, for the first time, that BAsO₄ exhibits negative linear compressibility (NLC) in its [001] direction at elevated pressures. This anomalous NLC property has been explained this effect in terms of deformation of the tetrahedra which appear to be the main cause of this effect [12] (the ‘demi wine-rack mechanism’), rather than the more conspicuous tetrahedral tilting about a local axis parallel to the [001] direction. A more recent study, which looked at uniaxial loading with a focus on the [001] direction, has further confirmed this amenable mode of deformation thus re-enforcing this explanation [27].

In addition to the experimentally measured NLC, BAsO₄ has also been shown to exhibit a negative Poisson’s ratio in the (001) plane [11,12], see Figure 1c, a property which was explained in terms of tetrahedral rotations about a local axis parallel to the orthogonal [001] direction. More specifically, as shown in Figure 1, the tetrahedra project in the auxetic (001) plane with a ‘rotating squares’ motif [37], a geometry which is amenable to auxeticity. More importantly, a DFT based study was able to confirm that uniaxial loading in the [100] or [010] direction is predicted to result in molecular-level deformations where these 2D projected squares appear to rotate relative to each other, i.e., deform with an auxeticity inducing deformation mechanism. Here, it should be noted that as discussed by Lisovenko et al. [67,68], in general, a tetragonal crystal such as BAsO₄ having a seven-constant tetragonal anisotropy (i.e., s₁₆ = −s₂₆ ≠ 0, see

Table 1. The elements of the compliance matrix of BAsO₄ as reported by Grima-Cornish et al. (2020) [12], which correspond to shear modulus G_xy = s₆₆⁻¹ = 40.5GPa. Note that these elements fulfil symmetry requirements for a tetragonal crystal [69]. After the application of standard axis transformations [69], the elements of this matrix may be used to compute the mechanical properties of the crystal in any direction/plane.

sij (x10−12 Pa−1)	j = 1	2	3	4	5	6
i = 1	10.4	3.4	−5.1	0.0	0.0	0.3
2	3.4	10.4	−5.1	0.0	0.0	−0.3
3	−5.1	−5.1	12.6	0.0	0.0	0.0
4	0.0	0.0	0.0	18.9	0.0	0.0
5	0.0	0.0	0.0	0.0	18.9	0.0
6	0.3	−0.3	0.0	0.0	0.0	24.7

) is expected to manifest some asymmetric aspects in the dependence of the Poisson’s ratio with the direction of loading compared to its six-constant tetragonal anisotropy counterpart (i.e., s₁₆ = −s₂₆ = 0). This is hardly visible in Figure 1c due to the fact that the magnitude of the coefficient s₁₆ = −s₂₆ is too small relative to the other coefficients to induce appreciable asymmetry in the xy-plane. Obviously, asymmetry is not precluded in the other planes.

However, despite these important discoveries, there are various aspects in the mechanical behaviour BAsO₄ that deserve further studies. For example, the idealised ‘rotating squares’ model, where the squares are perfectly rigid and simply rotate relative to each other [37], assumes that the system cannot shear in plane. This is the equivalent of having the compliance term s₆₆ = 0, which would correspond to an infinite shear modulus G_xy. (Prima facie, one might ask if zero value of s₆₆ is not in conflict with the stability of the system, which is typically granted by positive definiteness of the (Gibbs free) energy. The latter, inter alia, implies s₆₆ > 0. The rotating square model, however, is so internally constrained that it is only a ‘one degree of freedom’ model. Thus, since shear cannot be realised in the idealised ‘rotating squares’ model, the positivity of s₆₆ is not implementable—nor necessary.) However, as evident from Table 1, which lists the simulated compliances of BAsO₄, this is not the case with s₆₆ being predicted as 24.7 × 10⁻¹² Pa⁻¹, corresponding to a finite G_xy = s₆₆⁻¹ = 40.5GPa. Furthermore, although the idealised ‘rotating squares’ model predicts a constant in-plane Poisson’s ratio of −1 irrespective of the direction of loading, the negative Poisson’s ratio in the (001) plane of BAsO₄ is anisotropic and much lower in magnitude. In fact, although the DFT simulations suggest that single crystalline BAsO₄ is predicted to exhibit negative Poisson’s ratio in its (100) plane for loading in any direction (see Figure 1c), the Poisson’s ratio reached a maximum negative value of just c. −0.33 on-axis with auxeticity being at a minimum at 45° off-axis.

To explain such profile of mechanical properties, it is important to not only look at the behaviour of the material when it is subjected to uniaxial loading, but also when it is subjected to shear, something which so far has never been considered in detail in the case of BAsO₄. In view of this, the present paper, with the help of results from simulation, will analyse the shear properties of boron arsenate with a special focus on the nanoscale molecular-level deformations in an attempt to explain better the trends in the anomalous Poisson’s ratio of this material.

3. Simulations and Analysis

The procedure used for the simulation is based on density functional theory (DFT) approach performed using the CASTEP code [70] as implemented within the Materials Studio V6.1 modelling environment (Accelrys, now Biova, 2019). Details of the simulations are provided elsewhere [12,27]. In particular, since this study will focus on shear deformations, meaning that the unit cell will not retain its 90° unit cell angles, due importance is to be given to the manner how the crystal is aligned in the 3D global space. For this work, the single crystal of boron arsenate with its standard $I\overline{4}$ symmetry reduced to P1 was aligned within the global coordinate system according to the convention adopted by the Institute of Radio Engineers (IRE) [71], i.e., in a manner that the [001] crystal direction is parallel to the global z-axis and the [010] crystal direction is aligned in the global yz-plane with no constraints being placed on the [100] direction. For this type of orientation, the unit cell has its shape defined by the triangular matrix H such that:

$$H=\left(\begin{array}{c}a\\ b\\ c\end{array}\right)=\left(\begin{array}{ccc}{a}_x& a_y& a_z\\ 0& b_y& b_z\\ 0& 0& c_z\end{array}\right)$$

where $a,b,c$ are the unit cell vectors which, for this type of orientation, are of the form [72]:

$$a=\left[a_x,a_y,a_z\right]=\left[a\sin \left(\beta \right)\sin \left(\gamma \right),-a\sin \left(\beta \right)\cos \left(\gamma \right),a\cos \left(\beta \right)\right]$$

$$b=\left[0,b_y,b_z\right]=\left[0,b\sin \left(\alpha \right),b\cos \left(\alpha \right)\right]$$

$$c=\left[0,0,c_z\right]=\left[0,0,c\right]$$

where $a,b,c$ are the unit cell lengths whilst $\alpha ,\beta ,\gamma$ are the unit cell angles. Here, it should be noted that the reduction of symmetry from the standard $I\overline{4}$ symmetry (which can be described in terms of just three independent atoms, as reported by Schulze in 1934 [2]) reduced to P1 results in a system with twelve independent atoms: two boron atoms, labelled B1 and B2; two arsenic atoms, labelled A1 and A2; and eight oxygen atoms, labelled O1, O2,…, O8) as detailed in Figure 2 and

Table 2. (a) Boron arsenate as reported by Schulze [2]. (b) Relationship between the original fractional coordinates with symmetry applied and the fractional coordinates of the P1 system, together with the unique labelling system used to identify each atom within the unit cell.

(a) The crystal structure of boron arsenate as reported by Schulze.
Unit Cell Properties
Space Group				$I\overline{4}$
a = b				4.4580
c				6.7960
$\alpha =\beta =\gamma$				90°
Fractional Atomic Coordinates
B				0.00000			0.50000			0.25000
As				0.00000			0.00000			0.00000
O				0.16000			0.26000			0.14000
(b) The atom labelling system used in this work, see Figure 2.
Element	$I\overline{4}$ Fract. Coordinates			Transformation Applied			P1 Fract. Coordinates			Label
As	0	0	0	x	y	z	0	0	0	A1
				−x	−y	z	0	0	0
				y	−x	−z	0	0	0
				−y	x	−z	0	0	0
				x + ½	y + ½	z + ½	0.5	0.5	0.5	A2
				−x + ½	−y + ½	−z + ½	0.5	0.5	0.5
				y + ½	−x + ½	−z + ½	0.5	0.5	0.5
				−y + ½	x + ½	−z + ½	0.5	0.5	0.5
B	0	0.5	0.25	x	y	z	0	0.5	0.25	B1
				−x	−y	z	0	0.5	0.25
				y	−x	−z	0.5	0	0.75	B2
				−y	x	−z	0.5	0	0.75
				x + ½	y + ½	z + ½	0.5	0	0.75
				−x + ½	−y + ½	−z + ½	0.5	0	0.75
				y + ½	−x + ½	−z + ½	0.5	0	0.75
				−y + ½	x + ½	−z + ½	0.5	0	0.75
O	0.16	0.26	0.14	x	y	z	0.16	0.26	0.14	O1
				−x	−y	z	0.84	0.74	0.14	O3
				y	-x	−z	0.26	0.84	0.86	O2
				−y	x	−z	0.74	0.16	0.86	O4
				x + ½	y + ½	z + ½	0.66	0.76	0.64	O5
				−x + ½	−y + ½	−z + ½	0.34	0.24	0.64	O7
				y + ½	−x + ½	−z + ½	0.76	0.34	0.36	O6
				−y + ½	x + ½	−z + ½	0.24	0.66	0.36	O8

.

All density functional theory (DFT) simulations were carried out with full periodic boundary conditions being applied using the Generalised Gradient Approximation with the Perdew–Burke–Ernzerhof exchange-correlation functional (GGA–PBE) [73,74] utilising a Monkhorst–Pack Grid of 6 × 6 × 4 (k-point separation was set to approximately 0.04 Å⁻¹). The overall charge and initial spin were set to zero. An energy cut-off value of 1200 eV was used. A geometry optimization was first carried out using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) minimiser with no external pressure p being applied where the convergence criteria utilised included an energy cut off per atom of 5 × 10⁻⁶ eV atom⁻¹ and a maximum stress on the crystal of 0.01 GPa. During the geometry optimisation process, the unit cell was being periodically realigned to the IRE default alignment [71] with no constraints on the atoms or crystal lattice being imposed, apart from those enforced through the periodic boundary conditions. Additional simulations were then performed with the system being subjected to shear stresses in the range of −10 GPa to +10 GPa using the same settings. These simulations were performed with the load being applied in an incremental manner in steps of 2GPa. These simulations compliment other simulations meant to study the behaviour of the system when subjected to uniaxial on-axis loading and a hydrostatic pressure p reported elsewhere [12,27]. Note that throughout the simulations, particularly at larger strains, care was taken to monitor the symmetry of the system, so as to ascertain, for example, that the system did not undergo a phase transition.

For all systems subjected to an external stress $\sigma =\left({\sigma }_{ij}\right)$, the Lagrangian strains were calculated from the unit cell matrix $H_{\sigma }$ at a stress $\sigma$ relative to $H_0$, the unit cell matrix system prior to any application of external stresses, using the procedure outlined by Schlenker et al. [72]. More specifically, the Lagrangian strain tensor $\epsilon =\left({\epsilon }_{ij}\right)$ for an applied stress $\sigma$ was obtained through:

$$\epsilon =\left({\epsilon }_{ij}\right)=\frac{1}{2}\left[{\left(H_0^{-1}{H}_{\sigma }\right)}^T+\left(H_0^{-1}{H}_{\sigma }\right)\right]-I_3$$

where I₃ is a 3 × 3 identity matrix and $\epsilon =\left({\epsilon }_{ij}\right)$ is the strain tensor. The on-axis uniaxial strains ${\epsilon }_i$ in the x, y and z directions and the shear strains ${\gamma }_{ij}$ in the ij plane were then calculated as:

$${\epsilon }_x={\epsilon }_{11}\hspace{1em}{\epsilon }_y={\epsilon }_{22}\hspace{1em}{\epsilon }_z={\epsilon }_{33}$$

$${\gamma }_{yz}={\epsilon }_{23}+{\epsilon }_{32}\hspace{1em}{\gamma }_{xz}={\epsilon }_{13}+{\epsilon }_{31}\hspace{1em}{\gamma }_{xy}={\epsilon }_{12}+{\epsilon }_{21}$$

Furthermore, in order to understand the way the crystal structure changes upon the application of a mechanical loads, a number of lengths and angles between atoms in both the three-dimensional space (i.e., bond lengths, length of sides of tetrahedra, bond angles, tetrahedral angles, angles between tetrahedra, torsion angles) and the two-dimensional projection in the xy-plane (i.e., lengths of ‘projected squares’, inter- and intra-square angles) were determined using vector algebra (see Figure 3).

More specifically, the length l between arbitrary atoms A and B, and the angle θ in 3D between arbitrary atoms A, B, and C, having Cartesian coordinates $\left(A_x,A_y,A_z\right)$, $\left(B_x,B_y,B_z\right)$ and $\left(C_x,C_y,C_z\right)$ respectively, are given by:

$$l=\left|\overrightarrow{AB}\right|=\left|\overrightarrow{OB}-\overrightarrow{OA}\right|$$

$$\theta ={\cos }^{-1}\left(\frac{\overrightarrow{BA}\cdot \overrightarrow{BC}}{\left|\overrightarrow{BA}\right|\left|\overrightarrow{BC}\right|}\right)={\cos }^{-1}\left(\frac{\left(\overrightarrow{OA}-\overrightarrow{OB}\right)\cdot \left(\overrightarrow{OC}-\overrightarrow{OB}\right)}{\left|\overrightarrow{OA}-\overrightarrow{OB}\right|\left|\overrightarrow{OC}-\overrightarrow{OB}\right|}\right)$$

where:

$$\overrightarrow{OA}=\left(\begin{array}{c}{A}_x\\ A_y\\ A_z\end{array}\right)\hspace{1em}\overrightarrow{OB}=\left(\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right)\hspace{1em}\overrightarrow{OC}=\left(\begin{array}{c}{C}_x\\ C_y\\ C_z\end{array}\right)$$

Thus since, for example:

$$\overrightarrow{AB}=\overrightarrow{OB}-\overrightarrow{OA}=\left(\begin{array}{c}{B}_x-A_x\\ B_y-A_y\\ B_z-A_z\end{array}\right),$$

the distance between A and B is given by:

$$l=\sqrt{{\left(B_x-A_x\right)}^2+{\left(B_y-A_y\right)}^2+{\left(B_z-A_z\right)}^2}.$$

whilst the angle between AB and BC around the common atom B is given by:

$$\theta ={\cos }^{-1}\left(\frac{\left(A_x-B_x\right)\left(C_x-B_x\right)+\left(A_y-B_y\right)\left(C_y-B_y\right)+\left(A_z-B_z\right)\left(C_z-B_z\right)}{\sqrt{{\left(A_x-B_x\right)}^2+{\left(A_y-B_y\right)}^2+{\left(A_z-B_z\right)}^2}\sqrt{{\left(C_x-B_x\right)}^2+{\left(C_y-B_y\right)}^2+{\left(C_z-B_z\right)}^2}}\right).$$

Similarly, the projected length l^* between A and B in the xy-plane is given by:

$$l^{*}=\sqrt{{\left(B_x-A_x\right)}^2+{\left(B_y-A_y\right)}^2}$$

whilst the projected angle ${\theta }^{*}$ between AB and AC, in the xy-plane are given by:

$${\theta }^{*}={\cos }^{-1}\left(\frac{\left(A_x-B_x\right)\left(C_x-B_x\right)+\left(A_y-B_y\right)\left(C_y-B_y\right)}{\sqrt{{\left(A_x-B_x\right)}^2+{\left(A_y-B_y\right)}^2}\sqrt{{\left(C_x-B_x\right)}^2+{\left(C_y-B_y\right)}^2}}\right)$$

Furthermore, various torsion angles φ related to connected atoms A–B–C–D, were also measured. The magnitudes of these angles were taken as the magnitude of the angles between the normal of the planes A–B–C and B–C–D, i.e., referring to Figure 3c:

$$\left|\varphi \right|={\cos }^{-1}\left(\frac{\left(\overrightarrow{BA}\times \overrightarrow{BC}\right)\cdot \left(\overrightarrow{CB}\times \overrightarrow{CD}\right)}{\left|\overrightarrow{BA}\times \overrightarrow{BC}\right|\left|\overrightarrow{CB}\times \overrightarrow{CD}\right|}\right)$$

The sign of these torsion angles was taken as the sign of the vector triple product $\left(\overrightarrow{BA}\times \overrightarrow{BC}\right)\times \left(\overrightarrow{CB}\times \overrightarrow{CD}\right)\cdot k$ where k is the unit vector in the z-direction so that a clockwise angle is assigned a positive sign and an anticlockwise angle is assigned a negative sign.

4. Results and Discussion

Images of the system when subjected to a shear load of +10GPa in the xy and xz planes are shown in Figure 4. To facilitate the discussion, these images of the sheared system, see Figure 4b,c, are shown alongside images of the equivalent unloaded system, see Figure 4a, and images of the system when subjected to uniaxial loads in the x and z direction, see Figure 4d,e. A quantification of the deformation is provided in the various plots shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, where Figure 5 reports the cell parameters; Figure 6 reports the Lagrangian strains; Figure 7 reports deformations in terms of changes in bond lengths, bond angles and torsion angles; Figure 8 reports deformations in terms of changes in to the tetrahedral shapes; and Figure 9 reports deformations in terms of changes in terms of the 2D ‘connected squares’ motif as projected by the 3D tetrahedra in the xy-plane. In all cases, it should be noted that, due to symmetry, the behaviour in the xz and yz planes are equivalent to each other (only the behaviour in the xz plane is reported here). Similarly, uniaxial loading in the x has the same effect loading in the y direction and thus only the behaviour for uniaxial loading in the x-direction is presented here. From these results, it is evident that the systems are not resistant to shear and, the molecular level deformations, which will be discussed in more detail below, are comparable in magnitude to those arising from uniaxial loading.

Before looking at the effect of the different loads, it is important to recall that BAsO₄, like other crystalline systems, should not be considered as a purely mechanical system. What is remarkable, nonetheless, is that although the simulation protocol looks at the system through a rather complex quantum mechanical formulation, the interpretation of the results can still, to quite a reasonable approximation, be performed through a rather simple mechanical interpretation that looks at BAsO₄ simply as a ‘mechanical structure’ obeying rather simple rules. General rules that seem to apply are that the deformations are predominantly due to changes in the angles, i.e., changes in angles between the tetrahedra (the B–O–As angles and in the O-B-O-As and O-As-O-B torsion angles) as well as changes in the O–M–O (M = B, As) within he tetrahedra themselves. Bond lengths also change, but to a lower extent where, in general, B-O bond lengths tend to change more than the As-O bonds. In fact, a general trend is that the smaller BO₄ tetrahedra have more pronounced bond length changes whilst the larger AsO₄ tetrahedra have more pronounced O–M–O angle changes. This difference in the manner of deformation of the BO₄ and AsO₄ tetrahedra can, at least in part, be due to actual differences in the size: from a purely mechanical perspective, the longer O–As bond is expected to be a more effective lever than the shorter O–B bond, thus facilitating changes in O–As–O angles compared to the O–B–O angles. From a chemistry perspective, it is also known that arsenic can form slightly more covalent and less polarised bonds with oxygen when compared to boron. Such a difference in the covalent character is due to the difference of the electronegativities of boron and arsenic (arsenic is slightly more electronegative than boron). In addition, the two atoms have a different oxidation state in this crystal structure (boron: 3+, arsenic: 5+) which contributes to make the BO₄ tetrahedra being smaller than the AsO₄ tetrahedra. All this results in slightly more “deformation resistant” O–As bonds when compared to O–B bonds, and a more deformable O–As–O bond angle as opposed to O–B–O.

Looking in more detail at the shear deformations in the x–y plane and attempting to interpret these both in terms of the 2D projections in the (001) x–y plane, one may note that shear loading in the plane of the squares does not result any appreciable relative rotation of the squares. Instead, shear loading in the x-y plane results in extensive deformations of the squares themselves which tend to adopt a parallelogram-like shape of the Type IIα form [75]. Such a shape may be considered as a sheared variation of the Type II connected rectangles [76], a motif which is known to be present in the (100) and (010) planes of α-cristobalite [77,78]. All this is very clearly evident from images of the deformed system in Figure 4 and from Figure 9i which quantifies the deformation in terms of the 2D model, from a ‘rotating squares’ perspective and from the quantification in terms of percentages changes in length shown in Figure 10i. These deformations may be even more visible if one had to look at the deformed structures where the BO₄ and AsO₄ are depicted as actual tetrahedra (see Figure 11). From a quantitative perspective, we note that the change in the side lengths of the projected ‘squares’ is so excessive that these change by as much as 25% from τ_xy = −10 GPa to τ_xy = +10 GPa. The ‘squares’ that deform the most are those which relate to the AsO₄ tetrahedra where the deformations and are caused primarily by changes in the O-M-O bond angles (which change by as much as c. 20° for a change from τ_xy = −10 GPa to τ_xy = +10 GPa) rather than the bond lengths, which change by less than 1%. Interestingly, the bond lengths do not follow a linear relation with applied shear stress and tend to increase in length for both a positive and a negative shear. In terms of angles, one may note that the 90° square angles also change in a non-insignificant manner when subjected to a shear, and much more than the angle between the ‘squares’ themselves, as evident to comparing Figure 9b-i and Figure 9d-i with Figure 9e-i. All this suggests that, from the perspective of the 2D projected squares, shear loading results primarily (and almost exclusively) in deformations of the ‘squares’ (which tend to become more like the ‘Type IIα parallelograms’ [75]) rather than their relative rotation. Thus, a comprehensive look at the deformations in the x-y plane suggests that while uniaxial loading in the [100] x-direction (and, by symmetry, uniaxial loading in the [010] y-direction) favours a ‘rotating squares’ mode of deformation leading to auxetic behaviour, shear loading in the xy-plane favours deformation of the squares themselves, a non-auxetic inducing mechanism. This finding is rather important as it explains both the anisotropy in the Poisson’s ratio and its deviation from the idealised value of −1 for rotating squares. Although the structural requirements are present (0 GPa system projects in the (001) plane as almost perfect squares, see Figure 1b, the mechanistic requirement to have perfectly rigid units which simply rotate relative to each other is not there. Such deviations from the ideal behaviour are to be expected, since the ‘squares’ in this system are mere 2D projections of 3D molecular tetrahedra which cannot be expected to behave as perfectly rigid units. In fact, if one had to look at the tetrahedra and how they deform when a shear load in the x-y plane is applied, see Figure 8i, one would notice that the four tetrahedra in the system are becoming quite irregular upon shearing, hence explaining the observation that their projection in the x-y plane no longer results in squares.

The magnitude of compliance term s₆₆ (=G_xy⁻¹) is rather high compared to the other s_ii signifying that the crystal is quite weak in shear in this plane. From an auxeticity perspective, this is rather unfortunate as the low-shear modulus indicates deviation from ideal behaviour for the maximisation of auxeticity. It is more than likely that, had the material been more resistant to shear deformation, the extent of auxeticity would have been more pronounced. This weakness in shear is also noticeable in the other (010) and (100) planes where, for example, we may note that the angle changes that occur on shearing in the x-z plane are some of the more pronounced angle changes for any form of loading. The same can be said for the bond length changes. In this case we note that, once again, there are non-negligible deformations of the tetrahedra, but this time, the most pronounced changes occur in the angles between the tetrahedra through changes in the B–O–As bond angles and the torsion angles. An interesting feature to note is that, although the shear load is being applied in the x–z plane, there are still some observable deformations in the x–y plane where, once again, the deformation projects primarily as deformations of ‘squares’ rather than their relative rotation. This once again emphasises the fact that the tetrahedra, or their 2D projections as ‘squares’, cannot be considered as rigid units. It is also interesting to note the rather large inter-tetrahedral bond angle changes do not project in the x-y plane as relative rotations of squares, which further emphasises the fact that this material should not be treated as a simple ‘unimode’ structure.

Before concluding it is important to note some of the strengths and limitations of this work. First and foremost, a strength of this work is that one would not have a comprehensive picture of how a material behaves under applied mechanical loads without looking at the behaviour of materials when these are sheared. The availably of such information and its analysis as presented here helps to complete this discussion on the properties of BAsO₄. In this respect, it is rather unfortunate that shear behaviour is very often ignored in such studies, particularly given the wealth of information that can be derived from it. In this particular case, this work has identified an important mechanistic feature related to the ‘rotating squares’ description of BAsO₄ in the auxetic (001) plane that would come in useful for formulation of analytical models of such materials. In fact, it has been shown that for a proper formulation of such a model, one would need to ensure that the ‘squares’ need to have the capability to become ‘Type IIα parallelograms’ [75] upon shearing. On the same note, this work provides further evidence that it is essential that any model for BAsO₄ would need to incorporate the possibility of the tetrahedral/square units to deform rather than behave as simple ‘rotating rigid units’. A limitation of this work is that it was based entirely on results of simulations. In this case, this was mitigated through the use of a DFT methodology, rather than a force-field based approach as used in the past, which methods are known to be less prone to transferability errors.

5. Conclusions

This work has presented and discussed the behaviour of BAsO₄ when subjected to shear loading thus enabling a better understanding of this material, which is now known be both auxetic and exhibiting NLC. From this work, further evidence was obtained that deformation of the BO₄ and AsO₄ tetrahedra, which project as ‘squares’ in the auxetic (001) plane, cannot be assumed to be negligible. In fact, it was now shown that such deformations become the most prominent mode of deformation upon shearing. It was further discussed that this can explain the particular values of Poisson’s ratio for this material, which are less than the −1 predicted for idealised models. All this provides highly valuable information should there be attempts to formulate analytical models for these types of materials.