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 & 0.0 & 0.0 \\ 0.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 1.0 \end{pmatrix}$$
$$\begin{split}axis_{global\ y} &= rot \cdot axis_{local\ y}\\ &= \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}\\ &= \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 & 0.0 & 0.0 \\ 0.0 & 0.0 & -1.0 \\ 0.0 & 1.0 & 0.0 \end{pmatrix}$$
$$\begin{split}axis_{global\ y} &= rot \cdot axis_{local\ y}\\ &= \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}\\ &= \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 & 0 & 0 & 0 \\ 0 & \beta & 0 & 0 \\ 0 & 0 & \gamma & 0 \\ 0 & 0 & 0 & 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 & 0 & 0 & t_x \\ 0 & 1 & 0 & t_y \\ 0 & 0 & 1 & t_z \\ 0 & 0 & 0 & 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 & 0 & 0 & 0 \\ 0 & cos(\theta) & -sin(\theta) & 0 \\ 0 & sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 0 & 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) & 0 & sin(\theta) & 0 \\ 0 & 1 & 0 & 0 \\ -sin(\theta) & 0 & cos(\theta) & 0 \\ 0 & 0 & 0 & 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) & -sin(\theta) & 0 & 0 \\ sin(\theta) & cos(\theta) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{pmatrix} \cdot \begin{pmatrix}x\\y\\z\\w\end{pmatrix}$$
The rotation matrix is:
$$\begin{split}R &= R_x \cdot R_y \cdot R_z \\ &= \begin{pmatrix}cos(\theta_y) \cdot cos(\theta_z) & -cos(\theta_y) \cdot sin(\theta_z) & sin(\theta_y) & 0 \\ sin(\theta_x) \cdot \sin(\theta_y) \cdot cos(\theta_z) + cos(\theta_x) \cdot sin(\theta_z) & -sin(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + cos(\theta_x) \cdot cos (\theta_z) & -sin(\theta_x) \cdot cos(\theta_y) & 0\\ -cos(\theta_x) \cdot sin(\theta_y) \cdot cos(\theta_z) + sin(\theta_x) \cdot sin(\theta_z) & cos(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + sin(\theta_x) + cos(\theta_z) & cos(\theta_x) \cdot cos(\theta_y) & 0 \\ 0 & 0 & 0 & 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) & -cos(\theta_y) \cdot sin(\theta_z) & sin(\theta_y) & t_x \\ sin(\theta_x) \cdot \sin(\theta_y) \cdot cos(\theta_z) + cos(\theta_x) \cdot sin(\theta_z) & -sin(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + cos(\theta_x) \cdot cos (\theta_z) & -sin(\theta_x) \cdot cos(\theta_y) & t_y\\ -cos(\theta_x) \cdot sin(\theta_y) \cdot cos(\theta_z) + sin(\theta_x) \cdot sin(\theta_z) & cos(\theta_x) \cdot sin(\theta_y) \cdot sin(\theta_z) + sin(\theta_x) + cos(\theta_z) & cos(\theta_x) \cdot cos(\theta_y) & t_z \\ 0 & 0 & 0 & 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}$$
