Is there a simpler way to get the Vectors of a local axis?
The two essential information that you need are the world matrix and the axis / vector in local space that you're interested in. In your case the latter is pointing in the direction of the local y-axis which is equivalent to the vector: 

$$axis_{local\ y} = \begin{pmatrix}0.0 \\ 1.0 \\ 0.0 \end{pmatrix}$$

The world matrix contains the information how the object is rotated in the global coordinate system. Therefore we need to decompose the world matrix to get the rotation matrix*, which can then be multiplied with the vector to get the direction of the local y-axis in the global coordinate system.

(*) Internally this is a quaternion

For example, if the object isn't rotated at all, the rotation matrix is:

$$rot = \begin{pmatrix}1.0 &amp; 0.0 &amp; 0.0 \\ 0.0 &amp; 1.0 &amp; 0.0 \\ 0.0 &amp; 0.0 &amp; 1.0 \end{pmatrix}$$

$$\begin{split}axis_{global\ y} &amp;= rot \cdot axis_{local\ y}\\ &amp;= \begin{pmatrix}1.0 \times 0.0 + 0.0 \times 1.0 + 0.0 \times 0.0 \\ 0.0 \times 0.0 + 1.0 \times 1.0 + 0.0 \times 0.0 \\ 0.0 \times 0.0 + 0.0 \times 1.0 + 0.0 \times 0.0 \end{pmatrix}\\ &amp;= \begin{pmatrix}0.0 \\ 1.0 \\ 0.0 \end{pmatrix}\end{split}$$

The axis doesn't change, which is what we would expect since the local coordinate system is aligned with the global coordinate system. When the object is rotated by 90° around the global x-axis, the local y-axis now points upwards in the global coordinate system:

$$rot = \begin{pmatrix}1.0 &amp; 0.0 &amp; 0.0 \\ 0.0 &amp; 0.0 &amp; -1.0 \\ 0.0 &amp; 1.0 &amp; 0.0 \end{pmatrix}$$

$$\begin{split}axis_{global\ y} &amp;= rot \cdot axis_{local\ y}\\ &amp;= \begin{pmatrix}1.0 \times 0.0 + 0.0 \times 1.0 + 0.0 \times 0.0 \\ 0.0 \times 0.0 + 0.0 \times 1.0 - 1.0 \times 0.0 \\ 0.0 \times 0.0 + 1.0 \times 1.0 + 0.0 \times 0.0 \end{pmatrix}\\ &amp;= \begin{pmatrix}0.0 \\ 0.0 \\ 1.0 \end{pmatrix}\end{split}$$

import bpy
import mathutils

obj = bpy.data.objects["Cube"] # Replace with your own object selection

(translation, rotation, scale) = obj.matrix_world.decompose()
local_y_axis = mathutils.Vector((0.0, 1.0, 0.0))
local_y_axis_global_coords = rotation @ local_y_axis 
print(local_y_axis_global_coords)


The approach works with any vector in local space, not just the axis.



If you're only interested in the $x$-, $y$- or $z$-axis, then you don't need to multiply with the rotation matrix as Robin Betts correctly noted. When the world matrix is normalized or the object scale is 1.0, then the local $x$-, $y$- and $z$-axis in global coordinates are equivalent to the first three column vectors.

import bpy
import mathutils

obj = bpy.data.objects["Cube"] # Replace with your own object selection

normalize_matrix_world = obj.matrix_world.normalized()

local_x_axis_global_coords = [row[0] for row in normalize_matrix_world[:3]]
local_y_axis_global_coords = [row[1] for row in normalize_matrix_world[:3]]
local_z_axis_global_coords = [row[2] for row in normalize_matrix_world[:3]]
print(local_x_axis_global_coords)
print(local_y_axis_global_coords)
print(local_z_axis_global_coords)


The world matrix is the result of multiplication of the scaling, rotation and translation matrices.

Scaling with matrix $S$
$$\begin{pmatrix}x'\\y'\\z'\\w'\end{pmatrix} = \begin{pmatrix}\alpha &amp; 0 &amp; 0 &amp; 0 \\ 0 &amp; \beta &amp; 0 &amp; 0 \\ 0 &amp; 0 &amp; \gamma &amp; 0 \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$

Translation with matrix $T$
$$\begin{pmatrix}x'\\y'\\z'\\w'\end{pmatrix} = \begin{pmatrix} 1 &amp; 0 &amp; 0 &amp; t_x \\ 0 &amp; 1 &amp; 0 &amp; t_y \\ 0 &amp; 0 &amp; 1 &amp; t_z \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$

Rotation around x-axis with $R_x$:
$$\begin{pmatrix}x'\\y'\\z'\\w'\end{pmatrix} = \begin{pmatrix} 1 &amp; 0 &amp; 0 &amp; 0 \\ 0 &amp; cos(\theta) &amp; -sin(\theta) &amp; 0 \\ 0 &amp; sin(\theta) &amp; cos(\theta) &amp; 0 \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$

Rotation around y-axis with $R_y$:
$$\begin{pmatrix}x'\\y'\\z'\\w'\end{pmatrix} = \begin{pmatrix} cos(\theta) &amp; 0 &amp; sin(\theta) &amp; 0 \\ 0 &amp; 1 &amp; 0 &amp; 0 \\ -sin(\theta) &amp; 0 &amp; cos(\theta) &amp; 0 \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$

Rotation around z-axis with $R_z$:
$$\begin{pmatrix}x'\\y'\\z'\\w'\end{pmatrix} = \begin{pmatrix} cos(\theta) &amp; -sin(\theta) &amp; 0 &amp; 0 \\ sin(\theta) &amp; cos(\theta) &amp; 0 &amp; 0 \\ 0 &amp; 0 &amp; 1 &amp; 0 \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$

The rotation matrix is:
$$\begin{split}R &amp;= R_x \cdot R_y \cdot R_z \\ &amp;= \begin{pmatrix}cos(\theta_y) \cdot cos(\theta_z) &amp; -cos(\theta_y) \cdot sin(\theta_z) &amp; sin(\theta_y) &amp; 0 \\ sin(\theta_x) \cdot \sin(\theta_y) \cdot cos(\theta_z) + cos(\theta_x) \cdot sin(\theta_z) &amp; -sin(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + cos(\theta_x) \cdot cos (\theta_z) &amp; -sin(\theta_x) \cdot cos(\theta_y) &amp; 0\\ -cos(\theta_x) \cdot sin(\theta_y) \cdot cos(\theta_z) + sin(\theta_x) \cdot sin(\theta_z) &amp; cos(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + sin(\theta_x) + cos(\theta_z) &amp; cos(\theta_x) \cdot cos(\theta_y) &amp; 0 \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix}\end{split}$$

The world matrix is:

$$ world = T \cdot R \cdot S  $$

If the scale factors $\alpha$, $\beta$ and $\gamma$ are $1$ or we normalize first, then $T \cdot R = T \cdot R \cdot S$ because the $S$ is the identity matrix:

$$ T \cdot R = \begin{pmatrix}cos(\theta_y) \cdot cos(\theta_z) &amp; -cos(\theta_y) \cdot sin(\theta_z) &amp; sin(\theta_y) &amp; t_x \\ sin(\theta_x) \cdot \sin(\theta_y) \cdot cos(\theta_z) + cos(\theta_x) \cdot sin(\theta_z) &amp; -sin(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + cos(\theta_x) \cdot cos (\theta_z) &amp; -sin(\theta_x) \cdot cos(\theta_y) &amp; t_y\\ -cos(\theta_x) \cdot sin(\theta_y) \cdot cos(\theta_z) + sin(\theta_x) \cdot sin(\theta_z) &amp; cos(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + sin(\theta_x) + cos(\theta_z) &amp; cos(\theta_x) \cdot cos(\theta_y) &amp; t_z \\ 0 &amp; 0 &amp; 0 &amp; 1\end{pmatrix}$$

The first three entries in the first three columns give us the local axis in global coordinates for the $x$-, $y$- and $z$- axis respectively. This is because if we multiply the local axis vector with the matrix, the result is the same as the column in the world matrix:

$$ axis_{local x} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0\end{pmatrix}$$

$$ axis_{local y} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0\end{pmatrix}$$

$$ axis_{local z} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0\end{pmatrix}$$

$$ world \cdot axis_{local\ x} = \begin{pmatrix}cos(\theta_y) \cdot cos(\theta_z) \\ sin(\theta_x) \cdot \sin(\theta_y) \cdot cos(\theta_z) + cos(\theta_x) \cdot sin(\theta_z) \\ -cos(\theta_x) \cdot sin(\theta_y) \cdot cos(\theta_z) + sin(\theta_x) \cdot sin(\theta_z)\\ 0\end{pmatrix}$$

$$ world \cdot axis_{local\ y} = \begin{pmatrix} -cos(\theta_y) \cdot sin(\theta_z) \\ -sin(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + cos(\theta_x) \cdot cos (\theta_z) \\ cos(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + sin(\theta_x) + cos(\theta_z) \\ 0 \end{pmatrix}$$

$$ world \cdot axis_{local\ z} = \begin{pmatrix} sin(\theta_y) \\ -sin(\theta_x) \cdot cos(\theta_y) \\ cos(\theta_x) \cdot cos(\theta_y) \\ 0 \end{pmatrix}$$